TABLE OF CONTENTS

MetropolisHastings/mh.frail [ Functions ]

NAME

    mh.frail --- MH for frailties

FUNCTION

Update the frailty estimates by Metropolis-Hastings

SYNOPSIS

2611 mh.frail <- function(hazard, frailty, regression)

INPUTS

    hazard     RCurve for hazard
    frailty    RCurve for frailty
    regression RRegression structure

OUTPUTS

    frailty    Rcurve with updated frailty values

SOURCE

2614 {
2615     acc <- rep(0, regression$m)
2616     # update each of the frailties separately
2617     for(i in 1:regression$m){
2618         u <- frailty$x[i]
2619         v <- frailty$tun
2620         # generate candidates from gamma distribution
2621         cand <- rgamma(1, shape = u^2 / v, scale = v / u)
2622         # check if the candidate is out of bounds or NaN
2623         if(is.nan(cand) || (frailty$hasspline && 
2624           (cand > attr(frailty$spline.knots, "b")[2] | 
2625             cand < attr(frailty$spline.knots, "b")[1]))) next;
2626         # choose another fraity to compensate for the change in the mean
2627         j <- i
2628         while(j == i) j <- floor(runif(1, 1, length(frailty$x) + 1))
2629         # adjust the second candidate to make sure the mean remains 1
2630         candj <- u + frailty$x[j] - cand
2631         # check that candj is in bounds as well.
2632         if(is.nan(candj) || (frailty$hasspline && 
2633           (candj > attr(frailty$spline.knots, "b")[2] | 
2634             candj < attr(frailty$spline.knots, "b")[1]))) next;
2635 
2636         # base likelihood
2637         baselik <- mklik.frail(i, hazard, frailty, regression) +
2638             mklik.frail(j, hazard, frailty, regression)
2639         temp <- frailty
2640         temp$x[i] <- cand
2641         temp$x[j] <- candj
2642         temp <- updatecurvex(temp, i)
2643         temp <- updatecurvex(temp, j)
2644         # candidate likelihoood
2645         candlik <- mklik.frail(i, hazard, temp, regression) +
2646             mklik.frail(j, hazard, temp, regression)
2647 
2648         # transition u->cand
2649         puc <- suppressWarnings(dgamma(cand, shape = u^2 / v, scale = v / u)) 
2650         # transition cand->u
2651         pcu <- suppressWarnings(dgamma(u, shape = cand^2 / v, scale = v / cand)) 
2652 
2653         acc[i] <- acceptreject(baselik, candlik, pcu / puc)
2654         if(acc[i]) frailty <- temp
2655     }
2656     frailty$accept <- mean(acc)
2657     return(frailty)
2658 }