TABLE OF CONTENTS
MetropolisHastings/mh.bdm [ Functions ]
NAME
mh.bdm --- RJMCMC for Birth-death-move steps
FUNCTION
This function handles the reversible-jump MCMC steps for adding, removing, or moving knots in the adaptive knot selection procedure. It can work with either the hazard or frailty curve.
SYNOPSIS
2896 mh.bdm <- function(which, hazard, frailty, regression)
INPUTS
which string, can be either "hazard" or "frailty"
hazard RCurve for hazard
frailty RCurve for frailty
regression RRegression structure
OUTPUTS
curve updated RCurve for either hazard or frailty
SOURCE
2899 { 2900 # get all the needed components 2901 if(which == "hazard") curve <- hazard 2902 if(which == "frailty") curve <- frailty 2903 if(!curve$hasspline) return(curve); 2904 knots <- curve$spline.knots 2905 ord <- curve$spline.ord 2906 params <- curve$spline.par 2907 nknots <- curve$spline.nknots 2908 candknots <- curve$spline.candknots 2909 ncandknots <- curve$spline.ncandknots 2910 occ <- attr(candknots, "occupied") 2911 occind <- which(occ == 1) 2912 # evaluate the prior probability of k knots, k+1 knots, and k-1 knots 2913 pk <- nknotsPrior(nknots, curve) 2914 pkp1 <- if(nknots < curve$spline.maxoccknots) nknotsPrior(nknots + 1, curve) else 0 2915 pkm1 <- if(nknots > 1) nknotsPrior(nknots - 1, curve) else 0 2916 # compute probability of birth, death and move steps 2917 pb <- curve$spline.bdmconst * min(1, pkp1 / pk) # P(birth) 2918 pd <- curve$spline.bdmconst * min(1, pkm1 / pk) # P(death) 2919 pm <- max(0, 1 - pb - pd) # P(move) 2920 u <- runif(1) 2921 if(curve$name == "hazard") 2922 baselik <- mklik.spline.haz(curve$spline.par, hazard, frailty, regression) 2923 if(curve$name == "frailty") 2924 baselik <- mklik.spline.frail(curve$spline.par, hazard, frailty, regression) 2925 if(u < pd){ 2926 # Death step 2927 j <- sample(1:nknots, 1) + ord # index of dying knot 2928 x <- knots[j] 2929 #cat("Remove knot ", j, " at ", knots[j], "(occ", occ[occind[j - ord]], ")\n") 2930 knots2 <- knots 2931 params2 <- params 2932 knots2 <- knots2[ - j] # remove the knot 2933 # update parameters symetrically to birth step 2934 params2 <- params2[ - (j - 1)] 2935 if(ord > 2) for(j2 in (j - ord + 1):(j - 2)){ 2936 r2 <- (x - knots2[j2]) / (knots2[j2 + ord - 1] - knots2[j2]) 2937 inner <- 1 / r2 *exp(params2[j2]) - (1 - r2) / r2 * exp(params2[j2 - 1]) 2938 if(inner > 0) params2[j2] <- log(inner) else params2[j2] <- curve$spline.min 2939 } 2940 attr(knots2, "boundary") <- attr(knots, "boundary") 2941 attr(knots2, "order") <- attr(knots, "order") 2942 attr(knots2, "index") <- 1:(nknots - 1) - ord 2943 # create a temporary curve that is identical to the original, except for 2944 # the spline.knots and spline.params components 2945 temp <- curve 2946 temp$spline.knots <- knots2 2947 temp$spline.nknots <- nknots - 1 2948 occ2 <- occ; occ2[occind[j - ord]] <- 0 2949 attr(temp$spline.candknots, "occupied") <- occ2 2950 # Compute the spline basis of the temp curve 2951 temp <- makesplinebasis(temp) 2952 # Jacobian of the transformation 2953 J <- exp(params2[j - 1]) / (exp(params[j - 1]) - exp(params[j - 2])) 2954 if(ord > 2) for(j2 in (j - ord + 2):(j - 2)) { 2955 r2 <- (x - knots2[j2]) / (knots2[j2 + ord - 1] - knots2[j2]) 2956 J <- J * exp(params2[j2]) / (r2 * exp(params[j2])) 2957 } 2958 # candidate likelihood 2959 if(curve$name == "hazard") { 2960 temp <- updatespline(temp, params2) 2961 candlik <- mklik.spline.haz(temp$spline.par, temp, frailty, regression) 2962 } 2963 if(curve$name == "frailty"){ 2964 # for the frailty, adjust the mean to make sure it is 1 2965 newmean <- exp(params2)%*%temp$spline.basisexp - 2966 exp(params2[j - 2]) * temp$spline.basisexp[j - 2] 2967 newinner <- -newmean / temp$spline.basisexp[j - 2] 2968 if(newinner < 0){ 2969 candlik<- -Inf 2970 }else{ 2971 params2[j - 2] <- log(newinner) 2972 temp <- updatespline(temp, params2) 2973 candlik <- mklik.spline.frail(temp$spline.par, hazard, temp, regression) 2974 } 2975 } 2976 # compute acceptance ratio, and accept-reject 2977 ratio <- sqrt(2 * pi * curve$spline.priorvar) * abs(J) 2978 acc <- acceptreject(baselik, candlik, ratio) 2979 #cat("Lik: ", baselik, candlik, J, ratio, acc, "\n") 2980 if(any(knots2[(ord + 1):(length(knots2) - ord)] != candknots[occ2 == 1])) browser() 2981 if(acc) return(temp) else return(curve) 2982 } 2983 2984 if(u > pd & u < pd + pb){ 2985 # Birth of a knot 2986 params <- curve$spline.par 2987 # choose an unoccupied candidate knot location at random 2988 birthind <- occind[1] 2989 while(occ[birthind] > 0) birthind <- floor(runif(1, 1, ncandknots + 1)) + ord 2990 # get the position and interval of candidate knot 2991 x <- candknots[birthind] 2992 i <- findInterval(x, knots) 2993 #cat("Birth at ", x, " after knot ", i, "\n"); 2994 # Compute new knots and parameters 2995 knots2 <- c(knots[1:i], x, knots[(i + 1):length(knots)]) 2996 attr(knots2, "boundary") <- attr(knots, "boundary") 2997 attr(knots2, "order") <- attr(knots, "order") 2998 attr(knots2, "index") <- 1:(nknots + 1) - ord 2999 params2 <- c(params[1:(i - ord + 1)], 0, params[(i - ord + 2):length(params)]) 3000 # Even though it is possible to add a knot non-destructively, we need to generate 3001 # another random number to maintain dimension matching and detailed balance 3002 for(i2 in (i - ord + 2):i){ 3003 r2 <- (x - knots[i2]) / (knots[i2 + ord - 1] - knots[i2]) 3004 if(i2 == i) r2 <- runif(1) 3005 params2[i2] <- log(r2 * exp(params[i2]) + (1 - r2) * exp(params[i2 - 1])) 3006 } 3007 # Create a temp curve with new knots and params 3008 temp <- curve 3009 temp$spline.knots <- knots2 3010 temp$spline.nknots <- nknots + 1 3011 occ2 <- occ; occ2[birthind] <- 1 3012 attr(temp$spline.candknots, "occupied") <- occ2 3013 temp <- makesplinebasis(temp) 3014 # Jacobian of the transformation 3015 J <- (exp(params[i]) - exp(params[i - 1])) / exp(params2[i]) 3016 if(ord > 2) for(i2 in (i - ord + 2):(i - 1)) { 3017 r2 <- (x - knots[i2]) / (knots[i2 + ord - 1] - knots[i2]) 3018 J <- J * r2 * exp(params[i2]) / exp(params2[i2]) 3019 } 3020 if(curve$name == "hazard") { 3021 temp <- updatespline(temp, params2) 3022 candlik <- mklik.spline.haz(temp$spline.par, temp, frailty, regression) 3023 } 3024 if(curve$name == "frailty"){ 3025 # adjust frailty mean 3026 newmean <- exp(params2)%*%temp$spline.basisexp - exp(params2[i]) * 3027 temp$spline.basisexp[i] 3028 newinner <- -newmean / temp$spline.basisexp[i] 3029 if(newinner < 0){ 3030 candlik<- -Inf 3031 }else{ 3032 params2[i] <- log(newinner) 3033 temp <- updatespline(temp, params2) 3034 candlik <- mklik.spline.frail(temp$spline.par, hazard, temp, regression) 3035 } 3036 } 3037 # acceptance ratio 3038 ratio <- 1 / sqrt(2 * pi * curve$spline.priorvar) * abs(J) 3039 acc <- acceptreject(baselik, candlik, ratio) 3040 #cat("Lik: ", baselik, candlik, J, ratio, acc, "\n") 3041 if(any(knots2[(ord + 1):(length(knots2) - ord)] != candknots[occ2 == 1])) browser() 3042 if(acc) return(temp) else return(curve) 3043 3044 } 3045 3046 if(u > pd + pb) { 3047 # Move a knot 3048 # choose a random knot to move 3049 moveind <- floor(runif(1, 1, nknots + 1)) 3050 # compute the range of candidate knot indices to which the knot can be moved 3051 leftknotind <- if(moveind == 1) ord + 1 else occind[moveind - 1] + 1 3052 rightknotind <- if(moveind == nknots) ncandknots + ord else occind[moveind + 1] - 1 3053 newknotind <- floor(runif(1, leftknotind, rightknotind + 1)) 3054 oldknotind <- occind[moveind] 3055 #cat("Moveind: ", moveind, "\n"); 3056 #cat("Old / New: ", candknots[oldknotind], candknots[newknotind], "\n"); 3057 if(newknotind == oldknotind) return(curve) 3058 # create new knots and parameters 3059 knots2 <- knots 3060 params2 <- params 3061 knots2[moveind + ord] <- candknots[newknotind] 3062 # update occupied knot indices 3063 occ2 <- occ 3064 occ2[newknotind] <- 1 3065 occ2[occind[moveind]] <- 0 3066 if(sum(occ == 1) < nknots) browser() 3067 if(any(diff(knots2) < 0)) browser() 3068 if(any(knots2[(ord + 1):(length(knots2) - ord)] != candknots[occ2 == 1])) browser() 3069 temp <- curve 3070 temp$spline.knots <- knots2 3071 attr(temp$spline.candknots, "occupied") <- occ2 3072 # update the spline basis 3073 temp <- makesplinebasis(temp) 3074 if(curve$name == "hazard") { 3075 candlik <- mklik.spline.haz(params2, temp, frailty, regression) 3076 } 3077 if(curve$name == "frailty"){ 3078 # adjust frailty mean 3079 newmean <- exp(params2)%*%temp$spline.basisexp - exp(params2[moveind + ord]) * 3080 temp$spline.basisexp[moveind + ord] 3081 newinner <- -newmean / temp$spline.basisexp[moveind + ord] 3082 if(newinner < 0){ 3083 candlik<- -Inf 3084 }else{ 3085 params2[moveind + ord] <- log(newinner) 3086 temp <- updatespline(temp, params2) 3087 candlik <- mklik.spline.frail(params2, hazard, temp, regression) 3088 } 3089 } 3090 acc <- acceptreject(baselik, candlik, 1) 3091 #cat("Lik: ", baselik, candlik, acc, "\n") 3092 if(acc) return(temp) else return(curve) 3093 } 3094 3095 }