TABLE OF CONTENTS

CMetropolisHastings/MH_Frail [ Functions ]

NAME

    MH_Frail --- MH step for frailties

FUNCTION

Update frailty estimates by Metropolis-Hastings. See also mh.frail.

SYNOPSIS

1924 void MH_Frail(curveP hazard, curveP frailty, regressionP regression)

INPUTS

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

SOURCE

1928 {
1929     double acc = 0;
1930     double baselik, candlik;
1931     int j;
1932     // update the frailties one at a time
1933     for(int i=0; i < regression->m; i++){
1934         double u = frailty->X[i];
1935         double y = frailty->Y[i];
1936         double v = frailty->tun[0];
1937         double uj,yj,candj;
1938         // generate candidate from gamma distribution
1939         double cand = rgamma( pow(u,2)/v, v/u);
1940         double pcu = 1; double puc = 1;
1941         // if the candidate is invalid (e.g. beyond boundary knots), fail
1942         if(isnan(cand) || ( frailty->hasSpline &&
1943             (cand > frailty->SplineKnots[frailty->nj + frailty->SplineOrd -1]
1944              | cand < frailty->SplineKnots[0]) ) ){
1945             continue;
1946         }else{
1947             // find another frailty to move by the same amount
1948             j = i;
1949             while(j == i){
1950                 j = (int) floor(runif(0,(double) frailty->nx));
1951             }
1952             //Rprintf("i: %d,  j: %d\n",i,j);
1953             uj = frailty->X[j];
1954             yj = frailty->Y[j];
1955             candj = u + uj - cand;
1956             // if the moved element is invalid, fail
1957             if( frailty->hasSpline &&
1958                 (candj > frailty->SplineKnots[frailty->nj + frailty->SplineOrd -1]
1959                  | candj < frailty->SplineKnots[0]) ) continue;
1960             // base likelihood
1961             baselik = LikelihoodFrailty(i, hazard, frailty, regression)
1962                 + LikelihoodFrailty(j, hazard, frailty, regression);
1963             // update the curve with the candidate, without touching the basis
1964             frailty->X[i] = cand;
1965             frailty->Y[i] = EvalCurveAtOnePoint(frailty, cand);
1966             frailty->X[j] = candj;
1967             frailty->Y[j] = EvalCurveAtOnePoint(frailty, candj);
1968             // candidate likelihood
1969             candlik = LikelihoodFrailty(i, hazard, frailty, regression)
1970                 + LikelihoodFrailty(j, hazard, frailty, regression);
1971             // transition probabilities
1972             puc = dgamma(cand, pow(u,2)/v, v/u, 0);
1973             pcu = dgamma(u, pow(cand,2)/v, v/cand, 0);
1974         }
1975         int thisacc = AcceptReject(baselik, candlik, pcu/puc);
1976         if(thisacc==0) { //Did not accept, so undo the damage
1977             frailty->X[i] = u;
1978             frailty->Y[i] = y;
1979             frailty->X[j] = uj;
1980             frailty->Y[j] = yj;
1981         }else{ // accepted, so update the curve fully
1982             UpdateCurveX(frailty, cand, i);
1983             UpdateCurveX(frailty, candj, j);
1984             for(int k=regression->Jicum[i]; k<regression->Jicum[i+1]; k++){
1985                 regression->frailrep[k] = frailty->X[i];
1986                 regression->frailelp[k] = frailty->X[i]*regression->elp[k];
1987             }
1988             for(int k=regression->Jicum[j]; k<regression->Jicum[j+1]; k++){
1989                 regression->frailrep[k] = frailty->X[j];
1990                 regression->frailelp[k] = frailty->X[j]*regression->elp[k];
1991             }
1992         }
1993         // update acceptance rate
1994         acc += (double) thisacc;
1995     }
1996     acc /= regression->m;
1997     frailty->Accept[0] = acc;
1998 }