CalculateEvidence: Calculates the model evidence from MCMC output
In rich-d-wilkinson/precision: Statistical methods for detecting non-binomial sex allocation when developmental mortality operates
Description
Usage
Arguments
Value
Author(s)
Examples
Description
Uses the method of Chib 1995 and Chib and Jeliazkov 2001 to calculate the log model evidence for each model. For the double and multiplicative binomial models, it is necessary to run a further MCMC analysis with p and psi fixed.
Usage
1
Arguments
mcmc.out
Output from the function MCMCWithinGibbs
data
Dataset
hyper
See the MCMCWithinGibbs help page
family
See the MCMCWithinGibbs help page
sd
logical; if TRUE, returns an estimate of the standard deviation of the log evidence. Generally unreliable due to extreme variance.
nbatch
See the MCMCWithinGibbs help page
step.size
See the MCMCWithinGibbs help page
thin
See the MCMCWithinGibbs help page
thinby
See the MCMCWithinGibbs help page
burnin
See the MCMCWithinGibbs help page
Value
Returns the logarithm of the model evidence.
Author(s)
Richard Wilkinson
Examples
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#demo(CalculateEvidence, package="precision")
#Load some data
data(GlegneriSecondary)
meelis.out <- Meelis.test(GlegneriSecondary, TwoSided=TRUE)
james.out = James.test(GlegneriSecondary, TwoSided=TRUE)
###########################################################
#### Define hyperparameters/priors
#### - the four built in secondary datasets have priors for d and lambda built in, which load when you load the data
#### If you are using a different datasets you'll need to run the commands
# hyper<-list()
## prior for the mortality rate is Gamma(alpha, beta)
# hyper$a.m <- 6 # replace the value with something sensible for your data
# hyper$b.m <- 10
## prior for average clutch size, lambda, is Gamma(alpha, beta)
# hyper$alpha.lambda <- 4
# hyper$beta.lambda <- 1
## prior for p is beta(a, b)
hyper$a.p <- 1 ## 20 ## Hyper parameters for p's prior distribution
hyper$b.p <- 1##100
# prior for psi - assumed to be Gaussian
hyper$mu.psi <- 0
hyper$sd.psi <- 1
#plot.prior(hyper=hyper, show=TRUE, family="multbinom")
#plot.prior(hyper=hyper, show=FALSE, filename="Priors.pdf", family="multbinom") # save plot rather than printing it to screen
##########################################
### MCMC parameters
nbatch <- 1*10^3 ## number of MCMC samples to generate - usually we'll require at least 10^5 iterations
thin = FALSE
burnin <- 10^4
thinby <- 3
#######################################################################################################
# BINOMIAL MODEL
#######################################################################################################
## Define a start point for the MCMC chain - important to name the parameter and the column of the NM matrix
b.theta0 <-c(10, 0.1, 0.5)
names(b.theta0) <-c("lambda", "p", "mort")
b.mcmc.out <- MCMCWithinGibbs( theta0=b.theta0, data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="binomial", keepNM=TRUE)
## Thin the chain or not
if(thin){
b.mcmc.out.t <- ThinChain(b.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
b.mcmc.out.t <- b.mcmc.out
}
rm(b.mcmc.out) ## delete as not needed and takes a lot of memory
#plot.posterior(chain=b.mcmc.out.t$chain, hyper=hyper, show=TRUE, family="binomial", theta.true=NULL)
#plot.trace(chain=b.mcmc.out.t$chain, show=TRUE, family="binomial")
#######################################################################################################
# MULTIPLICATIVE BINOMIAL MODEL
#######################################################################################################
## Define the Metropolis-Hastings random walk step size
m.step.size<-c(0.3,0.2) ## 0.3 and 0.2 aree
names(m.step.size) <- c("p.logit", "psi")
m.theta0 <-c(10, 0.1, 0, 0.1)
names(m.theta0) <-c("lambda", "p", "psi", "mort")
m.mcmc.out <- MCMCWithinGibbs( theta0=m.theta0, data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="multbinom", step.size=m.step.size, keepNM=TRUE)
if(thin){
m.mcmc.out.t <- ThinChain(m.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
m.mcmc.out.t <- m.mcmc.out
}
rm(m.mcmc.out) ## to save space
#plot.posterior(chain=m.mcmc.out.t$chain, hyper=hyper, show=TRUE, family="multbinom", theta.true=NULL)
#plot.trace(chain=m.mcmc.out.t$chain, show=TRUE, family="multbinom")
#######################################################################################################
# DOUBLE BINOMIAL MODEL
#######################################################################################################
d.step.size<-c(0.3,0.2) ## 0.3 and 0.2 aree
names(d.step.size) <- c("p.logit", "psi")
d.theta0 <-c(10, 0.1, 0, 0.1)
names(d.theta0) <-c("lambda", "p", "psi", "mort")
d.mcmc.out <- MCMCWithinGibbs( theta0=d.theta0, data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="doublebinom", step.size=d.step.size, keepNM=TRUE)
if(thin){
d.mcmc.out.t <- ThinChain(d.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
d.mcmc.out.t <- d.mcmc.out
}
rm(d.mcmc.out) ## to save space
#plot.posterior(chain=d.mcmc.out.t$chain, hyper=hyper, show=TRUE, family="doublebinom", theta.true=NULL)
#plot.trace(chain=d.mcmc.out.t$chain, show=TRUE, family="doublebinom")
############################## Estimate Model Evidences
b.log.evidence <- CalculateEvidence(mcmc.out=b.mcmc.out.t, data=GlegneriSecondary, hyper=hyper, family="binomial", sd=FALSE)
m.log.evidence <- CalculateEvidence(mcmc.out=m.mcmc.out.t, data=GlegneriSecondary, hyper=hyper, family="multbinom", sd=FALSE, nbatch=nbatch, step.size=m.step.size, thin=thin, thinby=thinby, burnin=burnin)
d.log.evidence <- CalculateEvidence(mcmc.out=d.mcmc.out.t, data=GlegneriSecondary, hyper=hyper, family="doublebinom", sd=FALSE, nbatch=nbatch, step.size=d.step.size, thin=thin, thinby=thinby, burnin=burnin)
## For the multiplicative and double models, we have to run an MCMC with theta fixed, so it will take longer
log.evidence <- c(b.log.evidence$log.evidence, m.log.evidence$log.evidence, d.log.evidence$log.evidence)
BF<-CalcBF(log.evidence)
chib.out <- list(BF=BF$BF, probH0 = BF$probH0 ,
ProbPosPsi = c("multbinom"=sum((m.mcmc.out.t$chain[,"psi"]>0))/length(m.mcmc.out.t$chain[,"psi"]), "doublebinom"=sum((d.mcmc.out.t$chain[,"psi"]>0))/length(d.mcmc.out.t$chain[,"psi"])),
log.BF=log(BF$BF), log.evidence=log.evidence,
R= c("R"=meelis.out$R.av),
meelis = c("U"=meelis.out$U.av, "p"=meelis.out$p.av, "conclusion"=ifelse(meelis.out$p.av<0.05, "RejectH0", "AcceptH0")),
james = c("U"=james.out$U, "p"=james.out$p.val, "conclusion" = ifelse(james.out$p.val<0.05, "RejectH0", "AcceptH0") ) )
print(chib.out)
rich-d-wilkinson/precision documentation built on May 27, 2019, 7:41 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value Author(s) Examples
Description
Uses the method of Chib 1995 and Chib and Jeliazkov 2001 to calculate the log model evidence for each model. For the double and multiplicative binomial models, it is necessary to run a further MCMC analysis with p and psi fixed.
Usage
1 |
Arguments
mcmc.out |
Output from the function |
data |
Dataset |
hyper |
See the |
family |
See the |
sd |
logical; if TRUE, returns an estimate of the standard deviation of the log evidence. Generally unreliable due to extreme variance. |
nbatch |
See the |
step.size |
See the |
thin |
See the |
thinby |
See the |
burnin |
See the |
Value
Returns the logarithm of the model evidence.
Author(s)
Richard Wilkinson
Examples
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 138 139 140 141 142 143 144 145 146 147 148 149 150 151 | #demo(CalculateEvidence, package="precision")
#Load some data
data(GlegneriSecondary)
meelis.out <- Meelis.test(GlegneriSecondary, TwoSided=TRUE)
james.out = James.test(GlegneriSecondary, TwoSided=TRUE)
###########################################################
#### Define hyperparameters/priors
#### - the four built in secondary datasets have priors for d and lambda built in, which load when you load the data
#### If you are using a different datasets you'll need to run the commands
# hyper<-list()
## prior for the mortality rate is Gamma(alpha, beta)
# hyper$a.m <- 6 # replace the value with something sensible for your data
# hyper$b.m <- 10
## prior for average clutch size, lambda, is Gamma(alpha, beta)
# hyper$alpha.lambda <- 4
# hyper$beta.lambda <- 1
## prior for p is beta(a, b)
hyper$a.p <- 1 ## 20 ## Hyper parameters for p's prior distribution
hyper$b.p <- 1##100
# prior for psi - assumed to be Gaussian
hyper$mu.psi <- 0
hyper$sd.psi <- 1
#plot.prior(hyper=hyper, show=TRUE, family="multbinom")
#plot.prior(hyper=hyper, show=FALSE, filename="Priors.pdf", family="multbinom") # save plot rather than printing it to screen
##########################################
### MCMC parameters
nbatch <- 1*10^3 ## number of MCMC samples to generate - usually we'll require at least 10^5 iterations
thin = FALSE
burnin <- 10^4
thinby <- 3
#######################################################################################################
# BINOMIAL MODEL
#######################################################################################################
## Define a start point for the MCMC chain - important to name the parameter and the column of the NM matrix
b.theta0 <-c(10, 0.1, 0.5)
names(b.theta0) <-c("lambda", "p", "mort")
b.mcmc.out <- MCMCWithinGibbs( theta0=b.theta0, data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="binomial", keepNM=TRUE)
## Thin the chain or not
if(thin){
b.mcmc.out.t <- ThinChain(b.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
b.mcmc.out.t <- b.mcmc.out
}
rm(b.mcmc.out) ## delete as not needed and takes a lot of memory
#plot.posterior(chain=b.mcmc.out.t$chain, hyper=hyper, show=TRUE, family="binomial", theta.true=NULL)
#plot.trace(chain=b.mcmc.out.t$chain, show=TRUE, family="binomial")
#######################################################################################################
# MULTIPLICATIVE BINOMIAL MODEL
#######################################################################################################
## Define the Metropolis-Hastings random walk step size
m.step.size<-c(0.3,0.2) ## 0.3 and 0.2 aree
names(m.step.size) <- c("p.logit", "psi")
m.theta0 <-c(10, 0.1, 0, 0.1)
names(m.theta0) <-c("lambda", "p", "psi", "mort")
m.mcmc.out <- MCMCWithinGibbs( theta0=m.theta0, data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="multbinom", step.size=m.step.size, keepNM=TRUE)
if(thin){
m.mcmc.out.t <- ThinChain(m.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
m.mcmc.out.t <- m.mcmc.out
}
rm(m.mcmc.out) ## to save space
#plot.posterior(chain=m.mcmc.out.t$chain, hyper=hyper, show=TRUE, family="multbinom", theta.true=NULL)
#plot.trace(chain=m.mcmc.out.t$chain, show=TRUE, family="multbinom")
#######################################################################################################
# DOUBLE BINOMIAL MODEL
#######################################################################################################
d.step.size<-c(0.3,0.2) ## 0.3 and 0.2 aree
names(d.step.size) <- c("p.logit", "psi")
d.theta0 <-c(10, 0.1, 0, 0.1)
names(d.theta0) <-c("lambda", "p", "psi", "mort")
d.mcmc.out <- MCMCWithinGibbs( theta0=d.theta0, data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="doublebinom", step.size=d.step.size, keepNM=TRUE)
if(thin){
d.mcmc.out.t <- ThinChain(d.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
d.mcmc.out.t <- d.mcmc.out
}
rm(d.mcmc.out) ## to save space
#plot.posterior(chain=d.mcmc.out.t$chain, hyper=hyper, show=TRUE, family="doublebinom", theta.true=NULL)
#plot.trace(chain=d.mcmc.out.t$chain, show=TRUE, family="doublebinom")
############################## Estimate Model Evidences
b.log.evidence <- CalculateEvidence(mcmc.out=b.mcmc.out.t, data=GlegneriSecondary, hyper=hyper, family="binomial", sd=FALSE)
m.log.evidence <- CalculateEvidence(mcmc.out=m.mcmc.out.t, data=GlegneriSecondary, hyper=hyper, family="multbinom", sd=FALSE, nbatch=nbatch, step.size=m.step.size, thin=thin, thinby=thinby, burnin=burnin)
d.log.evidence <- CalculateEvidence(mcmc.out=d.mcmc.out.t, data=GlegneriSecondary, hyper=hyper, family="doublebinom", sd=FALSE, nbatch=nbatch, step.size=d.step.size, thin=thin, thinby=thinby, burnin=burnin)
## For the multiplicative and double models, we have to run an MCMC with theta fixed, so it will take longer
log.evidence <- c(b.log.evidence$log.evidence, m.log.evidence$log.evidence, d.log.evidence$log.evidence)
BF<-CalcBF(log.evidence)
chib.out <- list(BF=BF$BF, probH0 = BF$probH0 ,
ProbPosPsi = c("multbinom"=sum((m.mcmc.out.t$chain[,"psi"]>0))/length(m.mcmc.out.t$chain[,"psi"]), "doublebinom"=sum((d.mcmc.out.t$chain[,"psi"]>0))/length(d.mcmc.out.t$chain[,"psi"])),
log.BF=log(BF$BF), log.evidence=log.evidence,
R= c("R"=meelis.out$R.av),
meelis = c("U"=meelis.out$U.av, "p"=meelis.out$p.av, "conclusion"=ifelse(meelis.out$p.av<0.05, "RejectH0", "AcceptH0")),
james = c("U"=james.out$U, "p"=james.out$p.val, "conclusion" = ifelse(james.out$p.val<0.05, "RejectH0", "AcceptH0") ) )
print(chib.out)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.