-
Notifications
You must be signed in to change notification settings - Fork 0
/
logit_bayes.R
64 lines (50 loc) · 1.91 KB
/
logit_bayes.R
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
library("MASS")
data("OME")
?OME # background on the data
head(OME)
any(is.na(OME)) # check for missing values
dat = subset(OME, OME != "N/A") # manually remove OME missing values identified with "N/A"
dat$OME = factor(dat$OME)
str(dat)
plot(dat$Age, dat$Correct / dat$Trials )
plot(dat$OME, dat$Correct / dat$Trials )
plot(dat$Loud, dat$Correct / dat$Trials )
plot(dat$Noise, dat$Correct / dat$Trials )
mod_glm = glm(Correct/Trials ~ Age + OME + Loud + Noise, data=dat,
weights=Trials, family="binomial")
summary(mod_glm)
plot(residuals(mod_glm, type="deviance"))
plot(fitted(mod_glm), dat$Correct/dat$Trials)
X = model.matrix(mod_glm)[,-1] # -1 removes the column of 1s for the intercept
head(X)
mod_string = " model {
for (i in 1:length(y)) {
y[i] ~ dbin(phi[i], n[i])
logit(phi[i]) = b0 + b[1]*Age[i] + b[2]*OMElow[i] + b[3]*Loud[i] + b[4]*Noiseincoherent[i]
}
b0 ~ dnorm(0.0, 1.0/5.0^2)
for (j in 1:4) {
b[j] ~ dnorm(0.0, 1.0/4.0^2)
}
} "
data_jags = as.list(as.data.frame(X))
data_jags$y = dat$Correct # this will not work if there are missing values in dat (because they would be ignored by model.matrix). Always make sure that the data are accurately pre-processed for JAGS.
data_jags$n = dat$Trials
str(data_jags) # make sure that all variables have the same number of observations (712).
params <- c('b0', 'b')
library(rjags)
mod <- jags.model(textConnection(mod_string), data = data_jags, n.chains = 3)
update(mod, 1000)
mod_sim <- coda.samples(mod, variable.names = params, n.iter = 5e3)
mod_csim <- do.call(rbind, mod_sim)
raftery.diag(mod_csim)
summary(mod_sim)
#predict new observation
mp <- colMeans(mod_csim)
xb <- mp[['b0']] + mp[['b[1]']]*60 + mp[['b[2]']]*0 + mp[['b[3]']]*50 + mp[['b[4]']]*0
1.0 / (1+exp(-xb))
#see model accuracy
pm_Xb = mp[['b0']] + X %*% mp[1:4]
phat = 1.0 / (1.0 + exp(-pm_Xb))
tab0.7 = table(phat > 0.7, (dat$Correct / dat$Trials) > 0.7)
sum(diag(tab0.7)) / sum(tab0.7)