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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
# function
# written by Liam J. Revell 2011, 2012, 2013
mcmcLambda<-function(tree,x,ngen=10000,control=list()){
# starting values (for now)
n<-length(tree$tip)
temp<-aggregate(x,list(species=as.factor(names(x))),mean)
xbar<-temp[,2]; names(xbar)<-temp[,1]; xbar<-xbar[tree$tip.label]
sig2<-mean(pic(xbar,tree)^2)
lambda<-1.0; max.lambda<-maxLambda(tree)
a<-mean(xbar)
intV<-mean(aggregate(x,list(species=as.factor(names(x))),var)[,2],na.rm=T)
prop<-c(0.01*sig2,0.02,0.01*sig2,rep(0.01*sig2*max(vcv(tree)),n),0.01*intV)
pr.mean<-c(1000,max.lambda/2,rep(0,n+1),1000)
pr.var<-c(pr.mean[1]^2,max.lambda,rep(1000,n+1),pr.mean[length(pr.mean)]^2)
# populate control list
con=list(sig2=sig2,lambda=lambda,a=a,xbar=xbar,intV=intV,pr.mean=pr.mean,pr.var=pr.var,prop=prop,sample=100)
con[(namc<-names(control))]<-control
con<-con[!sapply(con,is.null)]
# print control parameters to screen
message("Control parameters (set by user or default):"); str(con)
# function returns the log-likelihood
likelihood<-function(C,invC,detC,x,sig2,a,xbar,intV){
z<-xbar-a
logLik<--z%*%invC%*%z/(2*sig2)-nrow(C)*log(2*pi)/2-nrow(C)*log(sig2)/2-detC/2+sum(dnorm(x,xbar[names(x)],sd=sqrt(intV),log=T))
return(logLik)
}
# function returns the log prior probability
log.prior<-function(sig2,lambda,a,xbar,intV){
pp<-dexp(sig2,rate=1/con$pr.mean[1],log=T)+dunif(lambda,min=con$pr.mean[2]-con$pr.var[2]/2,max=con$pr.mean[2]+con$pr.var[2]/2,log=T)+sum(dnorm(c(a,xbar),mean=con$pr.mean[1+1:(n+1)],sd=sqrt(con$pr.var[1+1:(n+1)]),log=T))+dexp(intV,rate=1/con$pr.mean[length(con$pr.mean)],log=T)
return(pp)
}
# compute (starting values for) C
C1<-vcv.phylo(tree) # used for updates of lambda
C<-lambda.transform(con$lambda,C1)
invC<-solve(C)
detC<-determinant(C,logarithm=TRUE)$modulus[1]
# now set starting values for MCMC
sig2<-con$sig2; lambda<-con$lambda; a<-con$a; xbar<-con$xbar; intV<-con$intV
L<-likelihood(C,invC,detC,x,sig2,a,xbar,intV)
Pr<-log.prior(sig2,lambda,a,xbar,intV)
# store
X<-matrix(NA,ngen/con$sample+1,n+6,dimnames=list(NULL,c("gen","sig2","lambda","a",tree$tip.label,"var","logLik")))
X[1,]<-c(0,sig2,lambda,a,xbar,intV,L)
message("Starting MCMC...")
# start MCMC
for(i in 1:ngen){
j<-(i-1)%%(n+4)
if(j==0){
# update sig2
sig2.prime<-sig2+rnorm(n=1,sd=sqrt(con$prop[j+1]))
if(sig2.prime<0) sig2.prime<--sig2.prime
L.prime<-likelihood(C,invC,detC,x,sig2.prime,a,xbar,intV)
Pr.prime<-log.prior(sig2.prime,lambda,a,xbar,intV)
post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L),na.rm=T)
if(post.odds>runif(n=1)){
if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2.prime,lambda,a,xbar,intV,L.prime)
sig2<-sig2.prime
L<-L.prime
Pr<-Pr.prime
} else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar,intV,L)
} else if(j==1){
# update lambda
lambda.prime<-lambda+rnorm(n=1,sd=sqrt(con$prop[j+1]))
while(lambda.prime<0||lambda.prime>max.lambda){
if(lambda.prime<0) lambda.prime<--lambda.prime
if(lambda.prime>max.lambda) lambda.prime<-2*max.lambda-lambda.prime
}
# update C with new lambda
C.prime<-lambda.transform(lambda.prime,C1)
invC.prime<-solve(C.prime)
detC.prime<-determinant(C.prime,logarithm=TRUE)$modulus[1]
L.prime<-likelihood(C.prime,invC.prime,detC.prime,x,sig2,a,xbar,intV)
Pr.prime<-log.prior(sig2,lambda.prime,a,xbar,intV)
post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L),na.rm=T)
if(post.odds>runif(n=1)){
if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda.prime,a,xbar,intV,L.prime)
lambda<-lambda.prime
C<-C.prime; invC<-invC.prime; detC<-detC.prime
L<-L.prime
Pr<-Pr.prime
} else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar,intV,L)
} else if(j==2){
# update a
a.prime<-a+rnorm(n=1,sd=sqrt(con$prop[j+1]))
L.prime<-likelihood(C,invC,detC,x,sig2,a.prime,xbar,intV)
Pr.prime<-log.prior(sig2,lambda,a.prime,xbar,intV)
post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L),na.rm=T)
if(post.odds>runif(n=1)){
if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a.prime,xbar,intV,L.prime)
a<-a.prime
L<-L.prime
Pr<-Pr.prime
} else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar,intV,L)
} else if(j>2&&j<=(n+2)) {
k<-j-2 # update tip mean k
xbar.prime<-xbar
xbar.prime[k]<-xbar[k]+rnorm(n=1,sd=sqrt(con$prop[j+1]))
L.prime<-likelihood(C,invC,detC,x,sig2,a,xbar.prime,intV)
Pr.prime<-log.prior(sig2,lambda,a,xbar.prime,intV)
post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L),na.rm=T)
if(post.odds>runif(n=1)){
if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar.prime,intV,L.prime)
xbar<-xbar.prime
L<-L.prime
Pr<-Pr.prime
} else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar,intV,L)
} else if(j>(n+2)){
# update var
intV.prime<-intV+rnorm(n=1,sd=sqrt(con$prop[j+1]))
if(intV.prime<0) intV.prime<--intV.prime
L.prime<-likelihood(C,invC,detC,x,sig2,a,xbar,intV.prime)
Pr.prime<-log.prior(sig2,lambda,a,xbar,intV.prime)
post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L),na.rm=T)
if(post.odds>runif(n=1)){
if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar,intV.prime,L.prime)
intV<-intV.prime
L<-L.prime
Pr<-Pr.prime
} else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,lambda,a,xbar,intV,L)
}
}
# done MCMC
message("Done MCMC.")
return(X)
}
|
