EFPBSA-Beta-Prior/UpdateRPrw.R at main · Dydx1989/EFPBSA-Beta-Prior · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
UpdateRPrw <-
function(survObj, priorPara, mcmcPara, ini){

	n 				<- survObj$n
	p				<- survObj$p
	x				<- survObj$x

	J				<- priorPara$J
	ind.r			<- priorPara$ind.r
	ind.d			<- priorPara$ind.d
	ind.r_d			<- priorPara$ind.r_d
	numBeta			<- mcmcPara$numBeta
	beta.prop.me	<- mcmcPara$beta.prop.me
	beta.prop.var	<- mcmcPara$beta.prop.var

	xbeta			<- ini$xbeta
	be.ini			<- ini$beta.ini
	h				<- ini$h
	sd.be			<- ini$sd.be


	updatej	<- sample(1:p, numBeta)
	accept 	<- rep(0, p)


	for(j in updatej){

		be.prop <- be.ini

		xbeta[xbeta > 700] <- 700
		exp.xbeta		<- exp(xbeta)
		exp.xbeta.mat	<- matrix(rep(exp.xbeta, J), n, J)

		first.sum 	<- colSums(exp.xbeta.mat*ind.r_d)

		h.mat					<- matrix(rep(h, n), n, J, byrow = T)
		h.exp.xbeta.mat 	<- - h.mat * exp.xbeta.mat
		h.exp.xbeta.mat[h.exp.xbeta.mat > -10^(-7)] <- -10^(-7)
		second.sum 	<- colSums(log(1 - exp(h.exp.xbeta.mat))*ind.d)

		loglh.ini <- sum(-h*first.sum + second.sum)

		be.prop[j] <- rnorm(1, mean = beta.prop.me[j], sd = sqrt(beta.prop.var))

		xbeta.prop			<- xbeta - x[,j]*be.ini[j] + x[,j]*be.prop[j]
		xbeta.prop[xbeta.prop > 700] <- 700
		exp.xbeta.prop		<- exp(xbeta.prop)
		exp.xbeta.mat.prop	<- matrix(rep(exp.xbeta.prop, J), n, J)

		first.sum.prop 	<- colSums(exp.xbeta.mat.prop*ind.r_d)

		h.exp.xbeta.mat.prop 	<- - h.mat * exp.xbeta.mat.prop
		h.exp.xbeta.mat.prop[h.exp.xbeta.mat.prop > -10^(-7)] <- -10^(-7)
		second.sum.prop 	<- colSums(log(1 - exp(h.exp.xbeta.mat.prop))*ind.d)

		loglh.prop <- sum(-h*first.sum.prop + second.sum.prop)

		logprior.prop   <- dnorm(be.prop[j] , mean = 0 , sd = sd.be[j], log = TRUE)
		logprior.ini    <- dnorm(be.ini[j]  , mean = 0 , sd = sd.be[j], log = TRUE)
		logprop.prop    <- dnorm(be.prop[j] , mean = beta.prop.me[j], sd = sqrt(beta.prop.var), log = TRUE)
		logprop.ini     <- dnorm(be.ini[j]  , mean = beta.prop.me[j], sd = sqrt(beta.prop.var), log = TRUE)

		logR  <- loglh.prop - loglh.ini + logprior.prop - logprior.ini + logprop.ini - logprop.prop

		u = log(runif(1)) < logR

		if(u == 1){
			be.ini[j] <- be.prop[j]
			xbeta	<- xbeta.prop
			}

		accept[j]<-accept[j] + u

		} # end of for loop for j

	list(beta.ini = be.ini, accept = accept, xbeta = xbeta)

	}