TABLE OF CONTENTS

CMetropolisHastings/MH_SplineFrailty [ Functions ]

NAME

    MH_SplineFrailty --- MH for frailty spline parameters

FUNCTION

Update frailty spline parameters by Metropolis-Hastings. See also mh.frailty.spline.

SYNOPSIS

2150 void MH_SplineFrailty(curveP hazard, curveP frailty, regressionP regression)

INPUTS

    hazard        CCurve for the hazard
    frailty       CCurve for the frailty
    regression    CRegression structure

SOURCE

2154 {
2155     if(!frailty->hasSpline) return;
2156     double baselik, candlik;
2157     double sumacc=0;
2158     // base likelihood
2159     baselik = LikelihoodSplineFrailty(hazard,frailty,regression); 
2160     double * cand = (double *) calloc( frailty->nj, sizeof(double));
2161     // allocate storage for parameters of old spline
2162     double * oldPar = (double *) malloc( (frailty->nj) * sizeof(double));
2163     double * oldEPar = (double *) malloc( frailty->nj * sizeof(double));
2164     double * oldY = (double *) malloc( frailty->nx * sizeof(double));
2165     double * oldSplineY = (double *) malloc( frailty->nx * sizeof(double));
2166     dcopyWrapper(frailty->nj, frailty->SplinePar, oldPar);
2167     dcopyWrapper(frailty->nj, frailty->SplineEPar, oldEPar);
2168     dcopyWrapper(frailty->nx, frailty->Y, oldY);
2169     dcopyWrapper(frailty->nx, frailty->SplineY, oldSplineY);
2170     double oldSplineEParSum = frailty->SplineEParSum;
2171     double oldSplineFvar = frailty->SplineFvar;
2172     int ord2 = frailty->SplineOrd / 2;
2173 
2174     // update the frailty spline parameters one by one
2175     for(int j=0; j<frailty->nj; j++){
2176         dcopyWrapper(frailty->nj, oldPar, cand);
2177         // choose which parameter will compensate for j
2178         // to ensure the frailty mean remains 1
2179         int k = j;
2180         while(j == k | k<ord2-1 | k>frailty->nj-ord2-1) 
2181             k = (int) floor(runif(0,(double) frailty->nj));
2182         
2183         // Generate candidate parameter at j
2184         cand[j] = frailty->SplinePar[j]+frailty->SplineTun[0]*
2185             rnorm(0,frailty->SplineCandSD[j]);
2186         // Try to compute value at k to compensate for change at j
2187         double newmean = frailty->SplineBasisExp[j] * (exp(cand[j])-exp(oldPar[j]));
2188         double candk = log(oldEPar[k] - newmean/frailty->SplineBasisExp[k]);
2189         if(isnan(candk)) continue;
2190         cand[k] = candk;
2191         // Compute candidate likelihood
2192         UpdateSplinePar(frailty,cand,-1); 
2193         candlik = LikelihoodSplineFrailty(hazard,frailty,regression);
2194         int acc = AcceptReject(baselik, candlik, 1);
2195         if(acc){ // accepted, save and continue
2196             baselik = candlik;
2197             sumacc++;
2198             dcopyWrapper(frailty->nj, frailty->SplinePar, oldPar);
2199             dcopyWrapper(frailty->nj, frailty->SplineEPar, oldEPar);
2200             dcopyWrapper(frailty->nx, frailty->Y, oldY);
2201             dcopyWrapper(frailty->nx, frailty->SplineY, oldSplineY);
2202             oldSplineEParSum = frailty->SplineEParSum;
2203             oldSplineFvar = frailty->SplineFvar;
2204         }else{ // rejected, restore old values and continue
2205             dcopyWrapper(frailty->nj, oldPar, frailty->SplinePar);
2206             dcopyWrapper(frailty->nj, oldEPar, frailty->SplineEPar);
2207             dcopyWrapper(frailty->nx, oldY, frailty->Y);
2208             dcopyWrapper(frailty->nx, oldSplineY, frailty->SplineY);
2209             frailty->SplineEParSum = oldSplineEParSum ;
2210             frailty->SplineFvar = oldSplineFvar ;
2211         }
2212     }
2213     // acceptance rate
2214     frailty->SplineAccept[0] = sumacc / ((double) frailty->nj);
2215     // free memory
2216     free(cand);
2217     free(oldPar);
2218     free(oldEPar);
2219     free(oldY);
2220     free(oldSplineY);
2221 }