#Question 1: Standard Method Crude OR
X.E <- c(rep(1,147), rep(0,423), rep(1,58), rep(0,441))
Y.D <- c(rep(1,570), rep(0,499))
dput( list(N=1069, X.E=X.E, Y.D=Y.D), "data.txt")
mod <- glm(Y.D ~ X.E, family = binomial)
summary(mod)
##
## Call:
## glm(formula = Y.D ~ X.E, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5891 -1.1598 0.8156 1.1952 1.1952
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.04167 0.06806 -0.612 0.54
## X.E 0.97166 0.16934 5.738 9.58e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1477.2 on 1068 degrees of freedom
## Residual deviance: 1441.6 on 1067 degrees of freedom
## AIC: 1445.6
##
## Number of Fisher Scoring iterations: 4
#Question1: Bayes Approach Crude OR
library(rjags)
## Warning: package 'rjags' was built under R version 3.4.3
## Loading required package: coda
## Warning: package 'coda' was built under R version 3.4.3
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
library(coda)
model.1 <-
"model {
###data model
for (i in 1:N){
Y.D[i] ~ dbin(p[i], 1)
logit(p[i]) <- beta_0 + beta_1*X.E[i]
}
OR <- exp(beta_1)
pos.prob <- step(beta_1)
### prior
beta_1 ~ dnorm(0,0.0001)
beta_0 ~dnorm(0,0.0001)
}"
data.1 <- list(N =1069, X.E=X.E, Y.D=Y.D)
model_odds <-jags.model(textConnection(model.1), data=data.1)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 1069
## Unobserved stochastic nodes: 2
## Total graph size: 2152
##
## Initializing model
test_odds <-coda.samples(model_odds, c("OR", "beta_0", "beta_1"), n.iter= 10000)
summary(test_odds)
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## OR 2.68960 0.46534 0.0046534 0.0067037
## beta_0 -0.04239 0.06828 0.0006828 0.0009903
## beta_1 0.97473 0.17103 0.0017103 0.0024649
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## OR 1.9044 2.36028 2.64883 2.9764 3.70960
## beta_0 -0.1786 -0.08736 -0.04209 0.0037 0.09198
## beta_1 0.6441 0.85878 0.97412 1.0907 1.31092
graphics.off()
par("mar")
## [1] 5.1 4.1 4.1 2.1
par(mar=c(1,1,1,1))
plot(test_odds)
#Question2 Adjusted OR
X.E <- c( rep(1, 30), rep(0, 207), rep(1, 117), rep(0,216), rep(1, 33), rep(0, 327), rep(1, 25), rep(0, 114))
X.smoke <- c( rep(0, 237), rep(1, 333),
rep(0, 360), rep(1, 139))
Y.D <- c( rep(1, 570), rep(0, 499))
table(Y.D, X.E, X.smoke)
## , , X.smoke = 0
##
## X.E
## Y.D 0 1
## 0 327 33
## 1 207 30
##
## , , X.smoke = 1
##
## X.E
## Y.D 0 1
## 0 114 25
## 1 216 117
#Traditional
mod <- glm(Y.D ~ X.E+X.smoke, family=binomial)
summary(mod)
##
## Call:
## glm(formula = Y.D ~ X.E + X.smoke, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7796 -0.9776 0.6779 0.8998 1.3913
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.49008 0.08636 -5.675 1.39e-08 ***
## X.E 0.65856 0.17814 3.697 0.000218 ***
## X.smoke 1.18524 0.13406 8.841 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1477.2 on 1068 degrees of freedom
## Residual deviance: 1360.1 on 1066 degrees of freedom
## AIC: 1366.1
##
## Number of Fisher Scoring iterations: 4
#Bayes Approach
model.1 <-
"model {
###data model
for (i in 1:N){
Y.D[i] ~ dbin(p[i], 1)
logit(p[i]) <- beta_0 + beta_1*X.E[i]+beta_2*X.smoke[i]
}
OR <- exp(beta_1)
pos.prob <- step(beta_1)
### prior
beta_0 ~dnorm(0,0.0001)
beta_1 ~ dnorm(0,0.0001)
beta_2 ~dnorm(0,0.0001)
}"
data.1 <- list(N =1069, X.E=X.E, Y.D=Y.D, X.smoke=X.smoke)
model_odds <-jags.model(textConnection(model.1), data=data.1)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 1069
## Unobserved stochastic nodes: 3
## Total graph size: 3228
##
## Initializing model
#specify parameters for analysis
test_odds <-coda.samples(model_odds, c("OR", "beta_0", "beta_1", "beta_2"), n.iter= 10000)
summary(test_odds)
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## OR 1.9655 0.35360 0.0035360 0.005106
## beta_0 -0.4911 0.08629 0.0008629 0.001664
## beta_1 0.6599 0.17812 0.0017812 0.002585
## beta_2 1.1910 0.13393 0.0013393 0.002409
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## OR 1.3594 1.7179 1.9357 2.1776 2.7511
## beta_0 -0.6663 -0.5489 -0.4900 -0.4332 -0.3242
## beta_1 0.3071 0.5411 0.6605 0.7782 1.0120
## beta_2 0.9337 1.1015 1.1910 1.2805 1.4613
graphics.off()
par("mar")
## [1] 5.1 4.1 4.1 2.1
par(mar=c(1,1,1,1))
plot(test_odds)
traceplot(test_odds)
densplot(test_odds)
autocorr.plot(test_odds)
#Question3: Interaction
model.1 <-
"model{
for(i in 1:N){
Y.D[i] ~ dbin(p[i], 1)
logit(p[i]) <- beta_0+beta_1*X.E[i]+beta_2*X.smoke[i]+beta_3*X.E[i]*X.smoke[i]
}
OR.notcurrent <- exp(beta_1)
OR.current <- exp(beta_1+beta_3)
ROR <- exp(beta_3)
### prior
beta_0 ~ dnorm(0,0.0001)
beta_1 ~ dnorm(0,0.0001)
beta_2 ~ dnorm(0,0.0001)
beta_3 ~ dnorm(0,0.0001)
}"
data.1 <- list(N=1069, X.E=X.E, Y.D=Y.D, X.smoke=X.smoke)
model_odds <- jags.model(textConnection(model.1), data=data.1)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 1069
## Unobserved stochastic nodes: 4
## Total graph size: 3235
##
## Initializing model
test_odds<-coda.samples(model_odds, c("OR.notcurrent","OR.current","ROR","beta_1","beta_2","beta_3"),n.iter=10000)
summary(test_odds)
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## OR.current 2.5701 0.6660 0.006660 0.010636
## OR.notcurrent 1.4987 0.4130 0.004130 0.010967
## ROR 1.8471 0.7242 0.007242 0.019738
## beta_1 0.3677 0.2723 0.002723 0.007355
## beta_2 1.1066 0.1455 0.001455 0.003359
## beta_3 0.5439 0.3712 0.003712 0.010157
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## OR.current 1.5190 2.0941 2.4778 2.9523 4.0973
## OR.notcurrent 0.8435 1.2076 1.4489 1.7318 2.4370
## ROR 0.8397 1.3391 1.7103 2.2062 3.6363
## beta_1 -0.1702 0.1887 0.3708 0.5491 0.8908
## beta_2 0.8263 1.0079 1.1063 1.2017 1.4018
## beta_3 -0.1747 0.2920 0.5367 0.7913 1.2910
graphics.off()
par("mar")
## [1] 5.1 4.1 4.1 2.1
par(mar=c(1,1,1,1))
plot(test_odds)
traceplot(test_odds)
densplot(test_odds)
autocorr.plot(test_odds)
