TABLE OF CONTENTS

00main/SplineSurvMainLoop [ Functions ]

NAME

    SplineSurvMainLoop --- main loop in C

FUNCTION

This is the main fitting loop in C, which is called by splinesurv.agdata if it is called with usec=TRUE (default). Everything in the set of CFitting routines is implemented to match the RFitting routines, with the caveat of C's memory management.

See also the call graph linked from 00main.

SYNOPSIS

2734 SEXP SplineSurvMainLoop( SEXP Rhazard, SEXP Rfrailty, SEXP Rregression, 
2735         SEXP Rhistory, SEXP Rstartiter, SEXP Renditer, SEXP Rthin, 
2736         SEXP Rverbose)

INPUTS

    Rhazard       RCurve structure for the hazard
    Rfrailty      RCurve structure for the frailty
    Rregression   RRegression structure for the regression information
    Rhistory      RHistory structure for the history of the parameters
    Rstartiter    integer, for the starting iteration from which the loop should begin
    Renditer      integer, for the iteration until which the loop should run
    Rthin         integer, for the number of iterations that should be discarded each time
    Rverbose      integer, containing the amount of output that should be printed to the
                  screen. This is not as well-behaved as the verbose option within R.

SOURCE

2740 {
2741     
2742     // initialize random number generator.
2743     GetRNGstate();
2744     // Create local curve and regression structures, as CCurve, CRegression 
2745     // and CHistory structures. These consist primarily of pointers to the memory
2746     // allocated by R, but they make the addressing easier.
2747     curveP hazard = (struct curve *) R_alloc(1,sizeof(struct curve));
2748     curveP frailty = (struct curve *) R_alloc(1,sizeof(struct curve));
2749     regressionP regression = (struct regress *) R_alloc(1,sizeof(struct regress));
2750     historyP history = (struct hist *) R_alloc(1,sizeof(struct hist));
2751     hazard->isHazard=1;
2752     frailty->isHazard=0;
2753     PopulateLocalCurve(hazard, Rhazard);
2754     PopulateLocalCurve(frailty, Rfrailty);
2755     PopulateLocalRegression(regression, Rregression);
2756     PopulateLocalHistory(history,Rhistory);
2757     for(int j=0; j<regression->m; j++) 
2758         for(int i=regression->Jicum[j]; i<regression->Jicum[j+1]; i++)
2759             regression->frailrep[i] = frailty->X[j];
2760     // compute the value of the frailty*lp vector
2761     diagmvWrapper(regression->n, regression->frailrep, regression->elp, regression->frailelp);
2762 
2763     // Begin main loop
2764     int iter = asInteger(Rstartiter); // iteration counter
2765     int enditer = asInteger(Renditer);
2766     int thin = asInteger(Rthin);
2767     int verbose = asInteger(Rverbose);
2768     int iterEl = 0; // counts the elapsed iterations between history updates
2769     while(iter < enditer)
2770     {
2771         R_CheckUserInterrupt();
2772         iterEl++;
2773         if(verbose >= 4) Rprintf("%d", iterEl);
2774         
2775         // Metropolis-Hastings steps
2776         MH_Frail(hazard,frailty,regression);
2777 
2778         MH_Regression(hazard,frailty,regression);
2779         
2780         MH_SplineHazard(hazard,frailty,regression);
2781 
2782         MH_SplineFrailty(hazard,frailty,regression);
2783         
2784         MH_ParamHazard(hazard,frailty,regression);
2785 
2786         MH_ParamFrailty(hazard,frailty,regression);
2787 
2788         MH_Weight(hazard, hazard, frailty, regression);
2789         MH_Weight(frailty, hazard, frailty, regression);
2790 
2791         UpdatePostvarCurve(hazard);
2792         UpdatePostvarCurve(frailty);
2793         UpdatePostvarRegression(regression);
2794 
2795         // RJMCMC steps for adaptive knot selection
2796         if(hazard->SplineAdaptive) MH_BDM('h',hazard,frailty,regression);
2797         if(frailty->SplineAdaptive) MH_BDM('f',hazard,frailty,regression);
2798         
2799         // store current state in CHistory
2800         if(iterEl == thin){
2801             iter++;
2802             UpdateHistory(hazard, frailty, regression, history, iter);
2803             iterEl = 0;
2804             if(verbose>=3) Rprintf(" %d ",iter);
2805         }
2806         
2807     }
2808     REAL(getListElement(Rhazard,"spline.nknots"))[0] = (double) hazard->SplineNknots;
2809     REAL(getListElement(Rfrailty,"spline.nknots"))[0] = (double) frailty->SplineNknots;
2810     SEXP returnval = R_NilValue;
2811     PutRNGstate();
2812     return returnval;
2813 }