HW_4_Final
Muntasir
October 25, 2018
Loading Packages
library(readr)
library(car)
library(haven)
library(tidyr)
library(labelled)
library(ggplot2)
library(foreign)
library(lme4)
library(lmerTest)
library(sjstats)
library(MuMIn)
library(MASS)
library(arm)
library(lattice)
library(rstanarm)
library(bayesplot)
library(INLA)
library(brms)
library(pscl)
library(rstan)
library(brinla)Loading NHIS Data
load("P:/dem7473/HW Dem 7473/homework4.Rdata")Constructing regression models with 3 predictors
The Outcome Variable is mortality which is treated as a dichotomous variable with two possible outcomes: dead or alive. The three predictors are alcohol consumption, education and employment.
Bayesian Regression using STAN/MCMC
Logistic Regression Model
fit.glm<- glm(I(d.event)~alcohol+moredu+employment, data=women4, family = "binomial")
summary(fit.glm)$coef## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.6789828 0.2366013 -15.5492952 1.608663e-54
## alcoholHeavyDrinker -0.4333907 0.3809312 -1.1377139 2.552400e-01
## alcoholLifetimeAbstainer -0.1954936 0.2267801 -0.8620403 3.886653e-01
## alcoholLightDrinker -0.6326526 0.2102251 -3.0094055 2.617595e-03
## alcoholModerateDrinker -0.2088512 0.2935776 -0.7114004 4.768362e-01
## moreduHS 0.5332494 0.2069516 2.5766868 9.975229e-03
## moreduLHS 0.7288453 0.2586158 2.8182548 4.828548e-03
## moreduSomeCollege 0.2114438 0.2110215 1.0020013 3.163429e-01
## employmentUnemployed 0.1079063 0.3176641 0.3396867 7.340925e-01confint(fit.glm)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -4.1612284 -3.2324278
## alcoholHeavyDrinker -1.2404448 0.2709273
## alcoholLifetimeAbstainer -0.6356763 0.2565504
## alcoholLightDrinker -1.0358423 -0.2090407
## alcoholModerateDrinker -0.8020889 0.3558680
## moreduHS 0.1327596 0.9464979
## moreduLHS 0.2147069 1.2324914
## moreduSomeCollege -0.1984892 0.6314490
## employmentUnemployed -0.5735084 0.6836115Bayesian Model
options(mc.cores = parallel::detectCores())
#No Priors
fit.glm.00<-brm(I(d.event)|weights(mortwt)~alcohol+moredu+employment, data=women4, family=bernoulli(),
warmup=1000, #burnin for 500 interations for each chain = 1000 burnin
iter=5000, chains=2, #2*5000 = 10000 - 1000 burnin = 9000 total iterations
cores=2,seed = 1220 #number of computer cores, 1 per chain is good.
)## Compiling the C++ model## Start sampling#200 secNow I print the summaries of the model fit by likelihood and the Bayesian model:
summary(fit.glm) #Logistic Regression ##
## Call:
## glm(formula = I(d.event) ~ alcohol + moredu + employment, family = "binomial",
## data = women4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.3366 -0.2342 -0.2126 -0.1813 2.9411
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.6790 0.2366 -15.549 < 2e-16 ***
## alcoholHeavyDrinker -0.4334 0.3809 -1.138 0.25524
## alcoholLifetimeAbstainer -0.1955 0.2268 -0.862 0.38867
## alcoholLightDrinker -0.6327 0.2102 -3.009 0.00262 **
## alcoholModerateDrinker -0.2089 0.2936 -0.711 0.47684
## moreduHS 0.5332 0.2070 2.577 0.00998 **
## moreduLHS 0.7288 0.2586 2.818 0.00483 **
## moreduSomeCollege 0.2114 0.2110 1.002 0.31634
## employmentUnemployed 0.1079 0.3177 0.340 0.73409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1815.1 on 8230 degrees of freedom
## Residual deviance: 1787.8 on 8222 degrees of freedom
## AIC: 1805.8
##
## Number of Fisher Scoring iterations: 7print(summary(fit.glm.00), digits=3) #Bayesian Logistic Regression, no priors## Family: bernoulli
## Links: mu = logit
## Formula: I(d.event) | weights(mortwt) ~ alcohol + moredu + employment
## Data: women4 (Number of observations: 8231)
## Samples: 2 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept -3.758 0.004 -3.766 -3.751 5143
## alcoholHeavyDrinker -0.800 0.007 -0.815 -0.786 6731
## alcoholLifetimeAbstainer -0.189 0.004 -0.196 -0.181 6774
## alcoholLightDrinker -0.684 0.004 -0.691 -0.677 6143
## alcoholModerateDrinker -0.066 0.005 -0.075 -0.057 6758
## moreduHS 0.683 0.003 0.676 0.690 6196
## moreduLHS 0.835 0.005 0.826 0.844 5913
## moreduSomeCollege 0.319 0.004 0.312 0.326 6754
## employmentUnemployed 0.052 0.006 0.040 0.063 8480
## Rhat
## Intercept 1.000
## alcoholHeavyDrinker 1.000
## alcoholLifetimeAbstainer 1.000
## alcoholLightDrinker 1.000
## alcoholModerateDrinker 1.000
## moreduHS 1.000
## moreduLHS 1.000
## moreduSomeCollege 1.000
## employmentUnemployed 1.000
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).We can not see much dissimilarity between the two models, however, Bayesian shows more significant effects.
Specifying new priors for the analysis
prior1<-c(set_prior("normal(0,1)", class="b"))#prior for the beta's
# set_prior("normal(0,1)", class="Intercept"), #prior for the intercept
# set_prior("inv_gamma(.5,.5)", class="sigma"))#prior for the residual std. deviation
fit.glm.1<-brm(I(d.event)|weights(mortwt)~alcohol+moredu+employment, data=women4, family=bernoulli(),
prior = prior1,
warmup=1000, #burnin for 1000 interations for each chain = 1000 burnin
iter=5000, chains=2, #2*5000 = 10000 - 2000 burnin = 8000 total iterations
cores=2,seed = 1220, #number of computer cores, 1 per chain is good.
save_model="fit_lm_gamma_stan.txt")## Compiling the C++ model## Start sampling#150 secprint(summary(fit.glm.1), digits=3) #Bayesian Logistic Regression, with priors## Family: bernoulli
## Links: mu = logit
## Formula: I(d.event) | weights(mortwt) ~ alcohol + moredu + employment
## Data: women4 (Number of observations: 8231)
## Samples: 2 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept -3.759 0.004 -3.766 -3.751 3961
## alcoholHeavyDrinker -0.800 0.007 -0.814 -0.786 5360
## alcoholLifetimeAbstainer -0.189 0.004 -0.197 -0.181 4754
## alcoholLightDrinker -0.684 0.004 -0.691 -0.677 4274
## alcoholModerateDrinker -0.066 0.005 -0.075 -0.057 5298
## moreduHS 0.683 0.004 0.677 0.690 5505
## moreduLHS 0.835 0.005 0.826 0.844 5702
## moreduSomeCollege 0.319 0.004 0.312 0.326 5763
## employmentUnemployed 0.052 0.006 0.040 0.063 7005
## Rhat
## Intercept 1.000
## alcoholHeavyDrinker 1.000
## alcoholLifetimeAbstainer 1.000
## alcoholLightDrinker 1.000
## alcoholModerateDrinker 1.000
## moreduHS 1.000
## moreduLHS 1.000
## moreduSomeCollege 1.000
## employmentUnemployed 1.000
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).out00<-as.data.frame(fit.glm.00) #storing in dataframe
out1<-as.data.frame(fit.glm.1)
mysum<-function(x){
mu=mean(x,na.rm=T)
sdev=sd(x,na.rm=T)
lci=quantile(x, p=.025,na.rm=T) #lower credible interval
uci=quantile(x, p=.975,na.rm=T) #upper credible interval
list(mu=mu, sd=sdev, lci=lci, uci=uci)
}
sapply(out00, mysum )## b_Intercept b_alcoholHeavyDrinker b_alcoholLifetimeAbstainer
## mu -3.758469 -0.8001722 -0.188973
## sd 0.003978966 0.007388848 0.003840475
## lci -3.766152 -0.8149961 -0.1964875
## uci -3.750771 -0.7856767 -0.1812482
## b_alcoholLightDrinker b_alcoholModerateDrinker b_moreduHS b_moreduLHS
## mu -0.6838897 -0.06597456 0.6831812 0.8348106
## sd 0.003532439 0.004650322 0.003497126 0.004679034
## lci -0.6907853 -0.07511936 0.6762286 0.8255839
## uci -0.6769052 -0.05682343 0.6899329 0.8439341
## b_moreduSomeCollege b_employmentUnemployed lp__
## mu 0.3187993 0.05157707 -3043283
## sd 0.003600947 0.005764005 2.125012
## lci 0.3117603 0.04024672 -3043288
## uci 0.3258599 0.0630192 -3043280sapply(out1, mysum) ## b_Intercept b_alcoholHeavyDrinker b_alcoholLifetimeAbstainer
## mu -3.758559 -0.8000793 -0.1889428
## sd 0.004035721 0.007264764 0.003949477
## lci -3.76637 -0.8143867 -0.1967932
## uci -3.750633 -0.7857311 -0.1810738
## b_alcoholLightDrinker b_alcoholModerateDrinker b_moreduHS b_moreduLHS
## mu -0.6838818 -0.06595046 0.6832386 0.8348698
## sd 0.003585804 0.004640781 0.00350162 0.004646418
## lci -0.6910339 -0.07505439 0.6765019 0.8255492
## uci -0.6770319 -0.05690847 0.6900371 0.8437927
## b_moreduSomeCollege b_employmentUnemployed lp__
## mu 0.3188504 0.05157563 -3043291
## sd 0.003519751 0.005771402 2.160354
## lci 0.311953 0.040122 -3043297
## uci 0.3256713 0.06267143 -3043288The parameters show same values upto 3 decimals places, which means that the analysis is absolutely robust.
Bayesian P values
bayespval<-function(x){
if(mean(x)<0) {
pval=mean(x<0)
}
else{
pval=mean(x>0)
}
list(pval=pval)
}sapply(out1, bayespval)## $b_Intercept.pval
## [1] 1
##
## $b_alcoholHeavyDrinker.pval
## [1] 1
##
## $b_alcoholLifetimeAbstainer.pval
## [1] 1
##
## $b_alcoholLightDrinker.pval
## [1] 1
##
## $b_alcoholModerateDrinker.pval
## [1] 1
##
## $b_moreduHS.pval
## [1] 1
##
## $b_moreduLHS.pval
## [1] 1
##
## $b_moreduSomeCollege.pval
## [1] 1
##
## $b_employmentUnemployed.pval
## [1] 1
##
## $lp__.pval
## [1] 1Comparing Bayesian Models
WAIC1<-WAIC(fit.glm.00)
WAIC2<-WAIC(fit.glm.1)We can evaluate the model fit to the sensitivity of the prior from the normal model:
WAIC(fit.glm.00, fit.glm.1, compare = T)## WAIC SE
## fit.glm.00 6138906.91 370966.96
## fit.glm.1 6138860.19 370964.26
## fit.glm.00 - fit.glm.1 46.72 210.09It looks like the model fit is not that affected by the choice of prior in this case.
Additional Predictor:
Race and Ethnicity is included in the 3 predictors above as a fourth variable.
Logistic Regression Model
fit.glm.a<- glm(I(d.event)~alcohol+moredu+employment+race_eth, data=women4, family = "binomial")
summary(fit.glm)$coef## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.6789828 0.2366013 -15.5492952 1.608663e-54
## alcoholHeavyDrinker -0.4333907 0.3809312 -1.1377139 2.552400e-01
## alcoholLifetimeAbstainer -0.1954936 0.2267801 -0.8620403 3.886653e-01
## alcoholLightDrinker -0.6326526 0.2102251 -3.0094055 2.617595e-03
## alcoholModerateDrinker -0.2088512 0.2935776 -0.7114004 4.768362e-01
## moreduHS 0.5332494 0.2069516 2.5766868 9.975229e-03
## moreduLHS 0.7288453 0.2586158 2.8182548 4.828548e-03
## moreduSomeCollege 0.2114438 0.2110215 1.0020013 3.163429e-01
## employmentUnemployed 0.1079063 0.3176641 0.3396867 7.340925e-01confint(fit.glm)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -4.1612284 -3.2324278
## alcoholHeavyDrinker -1.2404448 0.2709273
## alcoholLifetimeAbstainer -0.6356763 0.2565504
## alcoholLightDrinker -1.0358423 -0.2090407
## alcoholModerateDrinker -0.8020889 0.3558680
## moreduHS 0.1327596 0.9464979
## moreduLHS 0.2147069 1.2324914
## moreduSomeCollege -0.1984892 0.6314490
## employmentUnemployed -0.5735084 0.6836115Bayesian Model
options(mc.cores = parallel::detectCores())
#No Priors
fit.glm.a00<-brm(I(d.event)|weights(mortwt)~alcohol+moredu+employment+race_eth, data=women4, family=bernoulli(),
warmup=1000, #burnin for 500 interations for each chain = 1000 burnin
iter=5000, chains=2, #2*5000 = 10000 - 1000 burnin = 9000 total iterations
thin = 1,
cores=2,seed = 1220 #number of computer cores, 1 per chain is good.
)## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling#105 SecPrinting the summaries of the model fit by likelihood and the Bayesian model:
summary(fit.glm.a) #Logistic Regression ##
## Call:
## glm(formula = I(d.event) ~ alcohol + moredu + employment + race_eth,
## family = "binomial", data = women4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4113 -0.2358 -0.2027 -0.1769 3.1334
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.28776 0.34339 -12.486 < 2e-16 ***
## alcoholHeavyDrinker -0.43998 0.38136 -1.154 0.248612
## alcoholLifetimeAbstainer -0.16542 0.22978 -0.720 0.471587
## alcoholLightDrinker -0.61399 0.21053 -2.916 0.003541 **
## alcoholModerateDrinker -0.21060 0.29382 -0.717 0.473508
## moreduHS 0.53887 0.20764 2.595 0.009454 **
## moreduLHS 0.88467 0.26693 3.314 0.000919 ***
## moreduSomeCollege 0.20395 0.21200 0.962 0.336038
## employmentUnemployed 0.03865 0.32022 0.121 0.903917
## race_ethNHBlack 0.93724 0.27636 3.391 0.000696 ***
## race_ethNHWhite 0.57868 0.26038 2.222 0.026254 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1815.1 on 8230 degrees of freedom
## Residual deviance: 1775.2 on 8220 degrees of freedom
## AIC: 1797.2
##
## Number of Fisher Scoring iterations: 7print(summary(fit.glm.a00), digits=3) #Bayesian Logistic Regression, no priors## Family: bernoulli
## Links: mu = logit
## Formula: I(d.event) | weights(mortwt) ~ alcohol + moredu + employment + race_eth
## Data: women4 (Number of observations: 8231)
## Samples: 2 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept -4.453 0.007 -4.467 -4.440 5834
## alcoholHeavyDrinker -0.808 0.007 -0.823 -0.794 7239
## alcoholLifetimeAbstainer -0.177 0.004 -0.184 -0.169 7688
## alcoholLightDrinker -0.675 0.004 -0.682 -0.669 7311
## alcoholModerateDrinker -0.072 0.005 -0.081 -0.063 7512
## moreduHS 0.682 0.004 0.675 0.689 7214
## moreduLHS 0.949 0.005 0.939 0.958 6718
## moreduSomeCollege 0.310 0.004 0.303 0.317 7232
## employmentUnemployed 0.010 0.006 -0.001 0.021 8756
## race_ethNHBlack 0.972 0.006 0.960 0.984 5904
## race_ethNHWhite 0.682 0.006 0.671 0.693 5860
## Rhat
## Intercept 1.000
## alcoholHeavyDrinker 1.000
## alcoholLifetimeAbstainer 1.000
## alcoholLightDrinker 1.000
## alcoholModerateDrinker 1.000
## moreduHS 1.000
## moreduLHS 1.000
## moreduSomeCollege 1.000
## employmentUnemployed 1.000
## race_ethNHBlack 1.000
## race_ethNHWhite 1.000
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).We see great similarity between the two models. Although the Bayesian model shows more significant effects.
Specifying new priors for the analysis
prior1<-c(set_prior("normal(0,1)", class="b"))#prior for the beta's
# set_prior("normal(0,1)", class="Intercept"), #prior for the intercept
# set_prior("inv_gamma(.5,.5)", class="sigma"))#prior for the residual std. deviation
fit.glm.a1<-brm(I(d.event)|weights(mortwt)~alcohol+moredu+employment+race_eth, data=women4, family=bernoulli(),
prior = prior1,
warmup=1000, #burnin for 1000 interations for each chain = 1000 burnin
iter=5000, chains=2, #2*5000 = 10000 - 2000 burnin = 8000 total iterations
thin = 1,
cores=2,seed = 1220, #number of computer cores, 1 per chain is good.
save_model="fit_lm_gamma_stan.txt")## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling#124 secprint(summary(fit.glm.a1), digits=3) #Bayesian Logistic Regression, with priors## Family: bernoulli
## Links: mu = logit
## Formula: I(d.event) | weights(mortwt) ~ alcohol + moredu + employment + race_eth
## Data: women4 (Number of observations: 8231)
## Samples: 2 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept -4.453 0.007 -4.467 -4.440 4936
## alcoholHeavyDrinker -0.808 0.007 -0.822 -0.794 7267
## alcoholLifetimeAbstainer -0.177 0.004 -0.184 -0.169 6609
## alcoholLightDrinker -0.675 0.004 -0.682 -0.668 6411
## alcoholModerateDrinker -0.072 0.005 -0.082 -0.063 7267
## moreduHS 0.682 0.003 0.675 0.689 7582
## moreduLHS 0.949 0.005 0.939 0.958 7033
## moreduSomeCollege 0.310 0.004 0.303 0.318 7954
## employmentUnemployed 0.010 0.006 -0.002 0.021 7946
## race_ethNHBlack 0.972 0.006 0.960 0.983 4910
## race_ethNHWhite 0.682 0.006 0.670 0.693 4977
## Rhat
## Intercept 1.000
## alcoholHeavyDrinker 1.000
## alcoholLifetimeAbstainer 1.000
## alcoholLightDrinker 1.000
## alcoholModerateDrinker 1.000
## moreduHS 1.000
## moreduLHS 1.000
## moreduSomeCollege 1.000
## employmentUnemployed 1.000
## race_ethNHBlack 1.000
## race_ethNHWhite 1.000
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).out000<-as.data.frame(fit.glm.a00)
out1a<-as.data.frame(fit.glm.a1)
mysum<-function(x){
mu=mean(x,na.rm=T)
sdev=sd(x,na.rm=T)
lci=quantile(x, p=.025,na.rm=T)
uci=quantile(x, p=.975,na.rm=T)
list(mu=mu, sd=sdev, lci=lci, uci=uci)
}
sapply(out000, mysum )## b_Intercept b_alcoholHeavyDrinker b_alcoholLifetimeAbstainer
## mu -4.453386 -0.8084124 -0.1766402
## sd 0.006787199 0.007300905 0.003968754
## lci -4.466814 -0.8229084 -0.184409
## uci -4.440114 -0.7939159 -0.1687983
## b_alcoholLightDrinker b_alcoholModerateDrinker b_moreduHS b_moreduLHS
## mu -0.6754719 -0.0723674 0.682053 0.9487533
## sd 0.003526414 0.004658845 0.003540484 0.00472207
## lci -0.6823777 -0.08128406 0.6749762 0.9393064
## uci -0.6685948 -0.06337779 0.6889084 0.9579816
## b_moreduSomeCollege b_employmentUnemployed b_race_ethNHBlack
## mu 0.3104788 0.009704734 0.971949
## sd 0.00361255 0.005779338 0.006073533
## lci 0.3033339 -0.001437887 0.9600891
## uci 0.317461 0.0211974 0.9838692
## b_race_ethNHWhite lp__
## mu 0.6817498 -3028293
## sd 0.005664759 2.33051
## lci 0.6708317 -3028298
## uci 0.6928859 -3028289sapply(out1a, mysum) ## b_Intercept b_alcoholHeavyDrinker b_alcoholLifetimeAbstainer
## mu -4.453294 -0.8083941 -0.1766347
## sd 0.006862014 0.007295046 0.003934714
## lci -4.466803 -0.822445 -0.1842223
## uci -4.439731 -0.7941895 -0.1686969
## b_alcoholLightDrinker b_alcoholModerateDrinker b_moreduHS b_moreduLHS
## mu -0.6754254 -0.0723267 0.6820526 0.9486821
## sd 0.003551389 0.004670658 0.003494446 0.004797601
## lci -0.6824541 -0.08163878 0.6751413 0.9393967
## uci -0.6684915 -0.06335021 0.6889235 0.9582757
## b_moreduSomeCollege b_employmentUnemployed b_race_ethNHBlack
## mu 0.3104538 0.009666257 0.9717578
## sd 0.003590197 0.005852722 0.006017397
## lci 0.3033954 -0.002038752 0.9599352
## uci 0.3175201 0.02109389 0.9834735
## b_race_ethNHWhite lp__
## mu 0.6816801 -3028304
## sd 0.00569307 2.294779
## lci 0.6704062 -3028309
## uci 0.6926277 -3028300Comparing Bayesian Models
WAIC1<-WAIC(fit.glm.a00)
WAIC2<-WAIC(fit.glm.a1)We can also evaluate the model fit to the sensitivity of the prior from the normal model:
WAIC(fit.glm.a00, fit.glm.a1, compare = T)## WAIC SE
## fit.glm.a00 6118793.0 369748.03
## fit.glm.a1 6119044.2 369766.96
## fit.glm.a00 - fit.glm.a1 -251.2 179.43It looks like the bayesian model fit is affected by the choice of priors and the inclusion of the additional predictor.
Bayesian Regression using INLA
Logistic regression in INLA
95% Bayesian Credible Interval
inla.setOption("enable.inla.argument.weights", TRUE)
mod3<-inla(I(d.event)~alcohol+moredu+employment+race_eth, weights = women4$mortwt,
data=women4,
family = "binomial", Ntrials = 1,
control.compute = list(waic=T),
num.threads = 2)
summary(mod3)##
## Call:
## c("inla(formula = I(d.event) ~ alcohol + moredu + employment + race_eth, ", " family = \"binomial\", data = women4, weights = women4$mortwt, ", " Ntrials = 1, control.compute = list(waic = T), num.threads = 2)" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.7194 4.6926 1.6401 9.0521
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -4.4530 0.0068 -4.4665 -4.4530 -4.4396
## alcoholHeavyDrinker -0.8084 0.0073 -0.8226 -0.8083 -0.7941
## alcoholLifetimeAbstainer -0.1766 0.0039 -0.1844 -0.1766 -0.1689
## alcoholLightDrinker -0.6755 0.0035 -0.6824 -0.6755 -0.6685
## alcoholModerateDrinker -0.0727 0.0047 -0.0818 -0.0727 -0.0635
## moreduHS 0.6818 0.0035 0.6749 0.6818 0.6887
## moreduLHS 0.9486 0.0048 0.9392 0.9486 0.9580
## moreduSomeCollege 0.3103 0.0036 0.3032 0.3103 0.3174
## employmentUnemployed 0.0097 0.0058 -0.0018 0.0097 0.0211
## race_ethNHBlack 0.9719 0.0061 0.9600 0.9719 0.9838
## race_ethNHWhite 0.6816 0.0057 0.6704 0.6815 0.6927
## mode kld
## (Intercept) -4.4530 0
## alcoholHeavyDrinker -0.8083 0
## alcoholLifetimeAbstainer -0.1766 0
## alcoholLightDrinker -0.6755 0
## alcoholModerateDrinker -0.0727 0
## moreduHS 0.6818 0
## moreduLHS 0.9486 0
## moreduSomeCollege 0.3103 0
## employmentUnemployed 0.0097 0
## race_ethNHBlack 0.9719 0
## race_ethNHWhite 0.6815 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 14.87(0.00)
## Number of equivalent replicates : 553.47
##
## Watanabe-Akaike information criterion (WAIC) ...: 1659027.44
## Effective number of parameters .................: 54657.82
##
## Marginal log-Likelihood: -3020856.50
## Posterior marginals for linear predictor and fitted values computedbri.fixed.plot(mod3)
90% Bayesian Credible Interval
inla.setOption("enable.inla.argument.weights", TRUE)
mod3.1<-inla(I(d.event)~alcohol+moredu+employment+race_eth, weights = women4$mortwt,
data=women4,
family = "binomial", Ntrials = 1,
control.compute = list(waic=T),
num.threads = 2, quantiles = c(0.05, 0.5, 0.95))
summary(mod3.1)##
## Call:
## c("inla(formula = I(d.event) ~ alcohol + moredu + employment + race_eth, ", " family = \"binomial\", data = women4, quantiles = c(0.05, 0.5, ", " 0.95), weights = women4$mortwt, Ntrials = 1, control.compute = list(waic = T), ", " num.threads = 2)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.7764 5.4418 1.9191 10.1373
##
## Fixed effects:
## mean sd 0.05quant 0.5quant 0.95quant
## (Intercept) -4.4530 0.0068 -4.4643 -4.4530 -4.4418
## alcoholHeavyDrinker -0.8084 0.0073 -0.8203 -0.8083 -0.7964
## alcoholLifetimeAbstainer -0.1766 0.0039 -0.1831 -0.1766 -0.1702
## alcoholLightDrinker -0.6755 0.0035 -0.6813 -0.6755 -0.6697
## alcoholModerateDrinker -0.0727 0.0047 -0.0804 -0.0727 -0.0650
## moreduHS 0.6818 0.0035 0.6760 0.6818 0.6876
## moreduLHS 0.9486 0.0048 0.9407 0.9486 0.9565
## moreduSomeCollege 0.3103 0.0036 0.3044 0.3103 0.3163
## employmentUnemployed 0.0097 0.0058 0.0001 0.0097 0.0193
## race_ethNHBlack 0.9719 0.0061 0.9619 0.9719 0.9818
## race_ethNHWhite 0.6816 0.0057 0.6722 0.6815 0.6909
## mode kld
## (Intercept) -4.4530 0
## alcoholHeavyDrinker -0.8083 0
## alcoholLifetimeAbstainer -0.1766 0
## alcoholLightDrinker -0.6755 0
## alcoholModerateDrinker -0.0727 0
## moreduHS 0.6818 0
## moreduLHS 0.9486 0
## moreduSomeCollege 0.3103 0
## employmentUnemployed 0.0097 0
## race_ethNHBlack 0.9719 0
## race_ethNHWhite 0.6815 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 14.87(0.00)
## Number of equivalent replicates : 553.47
##
## Watanabe-Akaike information criterion (WAIC) ...: 1659027.44
## Effective number of parameters .................: 54657.82
##
## Marginal log-Likelihood: -3020856.50
## Posterior marginals for linear predictor and fitted values computedbri.fixed.plot(mod3.1)
The WAIC values is: 1659027.44
Basic Random Intercept Model in INLA
95% Bayesian Credible Interval
inla.setOption("enable.inla.argument.weights", TRUE)
mod4<-inla(I(d.event)~alcohol+moredu+employment+race_eth+f(cohort2, model='iid'), weights = women4$mortwt,
data=women4,
family = "binomial", Ntrials = 1, #Ntrials=1 for bernoulli, logistic regression
verbose = TRUE,
control.compute = list(waic=T),
num.threads = 2)
summary(mod4)##
## Call:
## c("inla(formula = I(d.event) ~ alcohol + moredu + employment + race_eth + ", " f(cohort2, model = \"iid\"), family = \"binomial\", data = women4, ", " weights = women4$mortwt, Ntrials = 1, verbose = TRUE, control.compute = list(waic = T), ", " num.threads = 2)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.6381 18.4471 1.5759 22.6611
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -3.4645 0.5042 -4.4677 -3.4645 -2.4626
## alcoholHeavyDrinker -0.3901 0.0075 -0.4048 -0.3901 -0.3755
## alcoholLifetimeAbstainer -0.1020 0.0041 -0.1101 -0.1020 -0.0939
## alcoholLightDrinker -0.2822 0.0038 -0.2896 -0.2822 -0.2748
## alcoholModerateDrinker 0.3042 0.0049 0.2945 0.3042 0.3138
## moreduHS 0.4362 0.0036 0.4290 0.4362 0.4433
## moreduLHS 0.4547 0.0051 0.4448 0.4547 0.4646
## moreduSomeCollege 0.2008 0.0037 0.1935 0.2008 0.2081
## employmentUnemployed 0.3526 0.0060 0.3408 0.3526 0.3642
## race_ethNHBlack 0.8860 0.0061 0.8739 0.8860 0.8980
## race_ethNHWhite 0.0972 0.0058 0.0859 0.0972 0.1086
## mode kld
## (Intercept) -3.4645 0
## alcoholHeavyDrinker -0.3901 0
## alcoholLifetimeAbstainer -0.1020 0
## alcoholLightDrinker -0.2822 0
## alcoholModerateDrinker 0.3042 0
## moreduHS 0.4362 0
## moreduLHS 0.4547 0
## moreduSomeCollege 0.2008 0
## employmentUnemployed 0.3526 0
## race_ethNHBlack 0.8859 0
## race_ethNHWhite 0.0972 0
##
## Random effects:
## Name Model
## cohort2 IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for cohort2 0.3538 0.1301 0.1449 0.3381 0.6503 0.3046
##
## Expected number of effective parameters(std dev): 26.60(1e-04)
## Number of equivalent replicates : 309.44
##
## Watanabe-Akaike information criterion (WAIC) ...: 1647972.72
## Effective number of parameters .................: 99529.63
##
## Marginal log-Likelihood: -2659215.81
## Posterior marginals for linear predictor and fitted values computedbri.fixed.plot(mod4)
bri.hyperpar.plot(mod4)
90% Bayesian Credible Interval
mod4.1<-inla(I(d.event)~alcohol+moredu+employment+race_eth+f(cohort2, model='iid'), weights = women4$mortwt,
data=women4,
family = "binomial", Ntrials = 1, #Ntrials=1 for bernoulli, logistic regression
verbose = TRUE,
control.compute = list(waic=T),
num.threads = 2, quantiles = c(0.05, 0.5, 0.95))
summary(mod4.1)##
## Call:
## c("inla(formula = I(d.event) ~ alcohol + moredu + employment + race_eth + ", " f(cohort2, model = \"iid\"), family = \"binomial\", data = women4, ", " quantiles = c(0.05, 0.5, 0.95), weights = women4$mortwt, ", " Ntrials = 1, verbose = TRUE, control.compute = list(waic = T), ", " num.threads = 2)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 3.0292 18.2834 1.7494 23.0620
##
## Fixed effects:
## mean sd 0.05quant 0.5quant 0.95quant
## (Intercept) -3.4645 0.5042 -4.2890 -3.4645 -2.6406
## alcoholHeavyDrinker -0.3901 0.0075 -0.4024 -0.3901 -0.3778
## alcoholLifetimeAbstainer -0.1020 0.0041 -0.1088 -0.1020 -0.0952
## alcoholLightDrinker -0.2822 0.0038 -0.2884 -0.2822 -0.2760
## alcoholModerateDrinker 0.3042 0.0049 0.2960 0.3042 0.3123
## moreduHS 0.4362 0.0036 0.4302 0.4362 0.4421
## moreduLHS 0.4547 0.0051 0.4464 0.4547 0.4630
## moreduSomeCollege 0.2008 0.0037 0.1947 0.2008 0.2069
## employmentUnemployed 0.3526 0.0060 0.3427 0.3526 0.3624
## race_ethNHBlack 0.8860 0.0061 0.8759 0.8860 0.8960
## race_ethNHWhite 0.0972 0.0058 0.0877 0.0972 0.1067
## mode kld
## (Intercept) -3.4645 0
## alcoholHeavyDrinker -0.3901 0
## alcoholLifetimeAbstainer -0.1020 0
## alcoholLightDrinker -0.2822 0
## alcoholModerateDrinker 0.3042 0
## moreduHS 0.4362 0
## moreduLHS 0.4547 0
## moreduSomeCollege 0.2008 0
## employmentUnemployed 0.3526 0
## race_ethNHBlack 0.8859 0
## race_ethNHWhite 0.0972 0
##
## Random effects:
## Name Model
## cohort2 IID model
##
## Model hyperparameters:
## mean sd 0.05quant 0.5quant 0.95quant mode
## Precision for cohort2 0.3538 0.1301 0.1679 0.3381 0.5943 0.3046
##
## Expected number of effective parameters(std dev): 26.60(1e-04)
## Number of equivalent replicates : 309.44
##
## Watanabe-Akaike information criterion (WAIC) ...: 1647972.75
## Effective number of parameters .................: 99529.64
##
## Marginal log-Likelihood: -2659215.81
## Posterior marginals for linear predictor and fitted values computedbri.fixed.plot(mod4.1)
bri.hyperpar.plot(mod4.1)
The WAIC values is: 1647972.75
Random Slopes Model
coh<-table(women4$cohort2)
women4$coho<-rep(1:length(unique(women4$cohort2)), coh)
women4$coho2<- women4$coho #copy it
n<- length(unique(women4$cohort2)) Random Slopes - uncorrelated with random intercepts
95% Bayesian Credible Interval
mod6<-inla(I(d.event)~alcohol+moredu+employment+race_eth+f(coho, model="iid")+f( coho2, employment, model="iid"), weights = women4$mortwt,
data=women4,
family = "binomial", Ntrials = 1, #Ntrials=1 for bernoulli, logistic regression
control.compute = list(waic=T),
num.threads = 2)## Warning in `[<-.factor`(`*tmp*`, is.na(www), value = 0): invalid factor
## level, NA generatedsummary(mod6)##
## Call:
## c("inla(formula = I(d.event) ~ alcohol + moredu + employment + race_eth + ", " f(coho, model = \"iid\") + f(coho2, employment, model = \"iid\"), ", " family = \"binomial\", data = women4, weights = women4$mortwt, ", " Ntrials = 1, control.compute = list(waic = T), num.threads = 2)" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 3.0292 111.6244 1.6389 116.2925
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -5.2765 2.2339 -9.7805 -5.2554 -0.8947
## alcoholHeavyDrinker -0.8091 0.0073 -0.8235 -0.8091 -0.7948
## alcoholLifetimeAbstainer -0.1568 0.0040 -0.1645 -0.1568 -0.1490
## alcoholLightDrinker -0.6686 0.0036 -0.6755 -0.6686 -0.6616
## alcoholModerateDrinker -0.0947 0.0047 -0.1039 -0.0947 -0.0855
## moreduHS 0.6748 0.0035 0.6679 0.6748 0.6818
## moreduLHS 0.9451 0.0048 0.9357 0.9451 0.9546
## moreduSomeCollege 0.2937 0.0036 0.2866 0.2937 0.3008
## employmentUnemployed -6.1292 1.7172 -9.8474 -6.0225 -3.0167
## race_ethNHBlack 0.9829 0.0061 0.9710 0.9829 0.9948
## race_ethNHWhite 0.6867 0.0057 0.6755 0.6867 0.6979
## mode kld
## (Intercept) -5.2171 1e-04
## alcoholHeavyDrinker -0.8091 0e+00
## alcoholLifetimeAbstainer -0.1568 0e+00
## alcoholLightDrinker -0.6686 0e+00
## alcoholModerateDrinker -0.0947 0e+00
## moreduHS 0.6748 0e+00
## moreduLHS 0.9451 0e+00
## moreduSomeCollege 0.2937 0e+00
## employmentUnemployed -5.8308 2e-04
## race_ethNHBlack 0.9829 0e+00
## race_ethNHWhite 0.6867 0e+00
##
## Random effects:
## Name Model
## coho IID model
## coho2 IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for coho 0.0383 0.0173 0.0144 0.0351 0.0810 0.0292
## Precision for coho2 0.0420 0.0205 0.0143 0.0381 0.0932 0.0308
##
## Expected number of effective parameters(std dev): 36.28(0.0499)
## Number of equivalent replicates : 226.86
##
## Watanabe-Akaike information criterion (WAIC) ...: 1774644.25
## Effective number of parameters .................: 116592.17
##
## Marginal log-Likelihood: -2963945.43
## Posterior marginals for linear predictor and fitted values computedbri.fixed.plot(mod6)
bri.hyperpar.plot(mod6)
90% Bayesian Credible Interval
mod6.1<-inla(I(d.event)~alcohol+moredu+employment+race_eth+f(coho, model="iid")+f( coho2, employment, model="iid"), weights = women4$mortwt,
data=women4,
family = "binomial", Ntrials = 1, #Ntrials=1 for bernoulli, logistic regression
control.compute = list(waic=T),
num.threads = 2, quantiles = c(0.05, 0.5, 0.95))## Warning in `[<-.factor`(`*tmp*`, is.na(www), value = 0): invalid factor
## level, NA generatedsummary(mod6.1)##
## Call:
## c("inla(formula = I(d.event) ~ alcohol + moredu + employment + race_eth + ", " f(coho, model = \"iid\") + f(coho2, employment, model = \"iid\"), ", " family = \"binomial\", data = women4, quantiles = c(0.05, 0.5, ", " 0.95), weights = women4$mortwt, Ntrials = 1, control.compute = list(waic = T), ", " num.threads = 2)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 3.1982 113.7000 1.5066 118.4048
##
## Fixed effects:
## mean sd 0.05quant 0.5quant 0.95quant
## (Intercept) -5.2388 2.2322 -8.9235 -5.2223 -1.6166
## alcoholHeavyDrinker -0.8091 0.0073 -0.8211 -0.8091 -0.7971
## alcoholLifetimeAbstainer -0.1568 0.0040 -0.1633 -0.1568 -0.1502
## alcoholLightDrinker -0.6686 0.0036 -0.6744 -0.6686 -0.6627
## alcoholModerateDrinker -0.0947 0.0047 -0.1024 -0.0947 -0.0870
## moreduHS 0.6748 0.0035 0.6690 0.6748 0.6806
## moreduLHS 0.9451 0.0048 0.9372 0.9451 0.9531
## moreduSomeCollege 0.2937 0.0036 0.2878 0.2937 0.2997
## employmentUnemployed -6.1091 1.7101 -9.0719 -6.0058 -3.5045
## race_ethNHBlack 0.9829 0.0061 0.9729 0.9829 0.9929
## race_ethNHWhite 0.6867 0.0057 0.6773 0.6867 0.6961
## mode kld
## (Intercept) -5.1925 1e-04
## alcoholHeavyDrinker -0.8091 0e+00
## alcoholLifetimeAbstainer -0.1568 0e+00
## alcoholLightDrinker -0.6686 0e+00
## alcoholModerateDrinker -0.0947 0e+00
## moreduHS 0.6748 0e+00
## moreduLHS 0.9451 0e+00
## moreduSomeCollege 0.2937 0e+00
## employmentUnemployed -5.8199 1e-04
## race_ethNHBlack 0.9829 0e+00
## race_ethNHWhite 0.6867 0e+00
##
## Random effects:
## Name Model
## coho IID model
## coho2 IID model
##
## Model hyperparameters:
## mean sd 0.05quant 0.5quant 0.95quant mode
## Precision for coho 0.0383 0.0173 0.0167 0.0351 0.0711 0.0292
## Precision for coho2 0.0420 0.0205 0.0168 0.0381 0.0806 0.0308
##
## Expected number of effective parameters(std dev): 36.28(0.0499)
## Number of equivalent replicates : 226.86
##
## Watanabe-Akaike information criterion (WAIC) ...: 1767467.09
## Effective number of parameters .................: 113003.57
##
## Marginal log-Likelihood: -2963945.43
## Posterior marginals for linear predictor and fitted values computedbri.fixed.plot(mod6.1)
bri.hyperpar.plot(mod6.1)
The WAIC values is: 1788195.28
Compaing the logisitic, random intercept and random slopes model using INLA approach, the WAIC value shows that the Random intercept model is the best among the three. Which suggest that, more complicated models are not always the better model.