Both Homeworks of Multilevel Models 1 & 2
Urie Lomelí
April 8, 2018
The data set I have consists on individual characteristics from birth certificates of the six states of the Mexico-US Border in 2011. Birth certificates collected information of the mother (age, marital status, number of preganancies, if received prenatal care, education and institutional affiliation) and the newborn (sex, weight, place of birth, state and municipality of birth), I also include three variables at the municipality level, total homicides, male and female homicides of the period 2006-2010. The number of complete cases of this dataset (births2011) was 166 200 observations.
First we are going to use weight at birth as our outcome variables to test if the accumulation of homicides had an effect on the weight of newborns in these municipalities. Following what we saw in class, we need to see if there is any variation on the higher level units, the municipalities, by doing an ANOVA.
names(births2011)## [1] "age_m" "marital_s" "number_p" "prenatal_care"
## [5] "afiliation" "educ_m" "sex_birth" "wgt_birth"
## [9] "low_wgt" "place_birth" "state_birth" "mpio_birth"
## [13] "hom_t" "hom_m" "hom_f"births2011$wgt<-as.numeric(scale(births2011$wgt_birth, center = T, scale = T))
fit.anov<-lm(wgt~as.factor(mpio_birth), births2011)
anova(fit.anov)## Analysis of Variance Table
##
## Response: wgt
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(mpio_birth) 60 2594 43.231 43.901 < 2.2e-16 ***
## Residuals 166139 163605 0.985
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We see a significant variation in weight at birth across municipalities, so we can proceed to fit a hierarchical model.
Random intercept model
#install.packages("lmerTest")
#install.packages("arm")
library(lme4)## Loading required package: Matrixlibrary(lmerTest)## Warning: package 'lmerTest' was built under R version 3.4.4##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(arm)## Warning: package 'arm' was built under R version 3.4.4## Loading required package: MASS##
## arm (Version 1.9-3, built: 2016-11-21)## Working directory is C:/Users/gpe637/Desktop/UTSA/Spring 2018/Statistics IIfitb<-lmer(wgt~(1|mpio_birth), data=births2011)
arm::display(fitb, detail=T)## lme4::lmer(formula = wgt ~ (1 | mpio_birth), data = births2011)
## coef.est coef.se t value
## (Intercept) -0.03 0.04 -0.90
##
## Error terms:
## Groups Name Std.Dev.
## mpio_birth (Intercept) 0.27
## Residual 0.99
## ---
## number of obs: 166200, groups: mpio_birth, 61
## AIC = 469315, DIC = 469299.4
## deviance = 469304.2# Why it isn't 1.00 in Residual?
births2011$predwgt<-sd(births2011$wgt)*(fitted(fitb))+mean(births2011$wgt)
head(ranef(fitb))## $mpio_birth
## (Intercept)
## ACUÑA 0.1508974275
## AHUMADA -0.0115004756
## ALLENDE 0.0522772505
## ASCENSION 0.1845706163
## BALLEZA 0.0968194961
## BATOPILAS 0.0145441015
## BOCOYNA 0.0126240363
## BUENAVENTURA -0.0654630092
## CAJEME 0.0720117512
## CAMARGO 0.1096879766
## CHIHUAHUA -0.0969877125
## CIUDAD MADERO -0.6143897040
## CUATRO CIENEGAS 0.1781235106
## CUAUHTEMOC -0.0715907629
## DELICIAS 0.0296005360
## ENSENADA 0.1010102733
## FRANCISCO I. MADERO 0.1096115512
## GALEANA 0.2492437044
## GENERAL CEPEDA -0.1570634689
## GOMEZ FARIAS 0.0086705666
## GUACHOCHI 0.0081993279
## GUADALUPE -0.1362075830
## GUADALUPE Y CALVO -0.1294417997
## GUAZAPARES -0.0720198610
## GUERRERO -0.1957554567
## HERMOSILLO 0.0899898782
## HIDALGO DEL PARRAL 0.0965636598
## JANOS 0.2063999991
## JIMENEZ 0.1134783055
## JUAREZ -0.0084680164
## MATAMOROS 0.0232194580
## MEXICALI 0.1475078634
## MONCLOVA 0.0493076183
## MONTERREY -1.3572165253
## MUZQUIZ -0.0280908402
## NAMIQUIPA -0.0434876575
## NAVOJOA 0.2069986343
## NOGALES -0.0008083953
## NUEVO CASAS GRANDES 0.0203026263
## NUEVO LAREDO 0.1433812228
## OJINAGA 0.1195027981
## PARRAS 0.1153977682
## PIEDRAS NEGRAS -0.0508501914
## PLAYAS DE ROSARITO 0.2192883753
## RAMOS ARIZPE -0.1535974889
## RIVA PALACIO 0.3553335009
## SABINAS 0.0588453694
## SALTILLO -0.1391167168
## SAN BUENAVENTURA 0.0879938882
## SAN JUAN DE SABINAS 0.2142460114
## SAN LUIS RIO COLORADO 0.3391402641
## SAN NICOLAS DE LOS GARZA -0.0535602523
## SAN PEDRO 0.0966802899
## SAN PEDRO GARZA GARCIA -0.3240144926
## SAUCILLO -0.0499582779
## TAMPICO -0.4020023195
## TECATE 0.1833988798
## TIJUANA 0.1483976202
## TORREON 0.0179565777
## URIQUE 0.0194027871
## URUACHI -0.2890345147mpiomeans<-aggregate(cbind(wgt, predwgt)~mpio_birth, births2011, mean)
head(mpiomeans, n =15)## mpio_birth wgt predwgt
## 1 ACUÑA 0.118270834 0.117478903
## 2 AHUMADA -0.048016895 -0.044919000
## 3 ALLENDE 0.019503324 0.018858726
## 4 ASCENSION 0.176749463 0.151152092
## 5 BALLEZA 0.148162156 0.063400972
## 6 BATOPILAS -0.015236504 -0.018874423
## 7 BOCOYNA -0.020627037 -0.020794488
## 8 BUENAVENTURA -0.142546271 -0.098881534
## 9 CAJEME 0.061008552 0.038593227
## 10 CAMARGO 0.077885048 0.076269452
## 11 CHIHUAHUA -0.130497980 -0.130406237
## 12 CIUDAD MADERO -0.908594167 -0.647808228
## 13 CUATRO CIENEGAS 0.158339019 0.144704986
## 14 CUAUHTEMOC -0.105205298 -0.105009287
## 15 DELICIAS -0.003718891 -0.003817988plot(predwgt~wgt, mpiomeans)
We include individual level predictors: age of the mother, marital status, if she received prenatal care and education; and test if these variables are explaining birth weight.
# Model with individual level predictors:
fitb2<-lmer(wgt~age_m+marital_s+prenatal_care+educ_m+(1|mpio_birth), data = births2011)
arm::display(fitb2, detail=T)## lme4::lmer(formula = wgt ~ age_m + marital_s + prenatal_care +
## educ_m + (1 | mpio_birth), data = births2011)
## coef.est coef.se t value
## (Intercept) -0.41 0.04 -9.62
## age_m 0.01 0.00 19.71
## marital_sDIVORCIADA -0.13 0.05 -2.89
## marital_sSEPARADA -0.08 0.04 -2.11
## marital_sSOLTERA -0.05 0.01 -6.66
## marital_sUNIÓN LIBRE -0.01 0.01 -1.53
## marital_sVIUDA 0.09 0.06 1.51
## prenatal_careSI 0.19 0.02 12.15
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA -0.01 0.01 -1.11
## educ_mNINGUNA -0.04 0.02 -1.94
## educ_mPOSGRADO -0.19 0.04 -4.81
## educ_mPRIMARIA COMPLETA 0.01 0.01 0.74
## educ_mPRIMARIA INCOMPLETA -0.02 0.01 -1.60
## educ_mPROFESIONAL -0.07 0.01 -7.28
## educ_mSECUNDARIA COMPLETA 0.01 0.01 1.46
## educ_mSECUNDARIA INCOMPLETA -0.01 0.01 -1.22
##
## Error terms:
## Groups Name Std.Dev.
## mpio_birth (Intercept) 0.27
## Residual 0.99
## ---
## number of obs: 166200, groups: mpio_birth, 61
## AIC = 468670, DIC = 468412.1
## deviance = 468523.0# Likedihood ratio test to see if the individual level variables are explaining anything about out ourcome. Spoiler alert: Yes.
anova(fitb, fitb2)## refitting model(s) with ML (instead of REML)## Data: births2011
## Models:
## object: wgt ~ (1 | mpio_birth)
## ..1: wgt ~ age_m + marital_s + prenatal_care + educ_m + (1 | mpio_birth)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 3 469310 469340 -234652 469304
## ..1 18 468559 468739 -234262 468523 781.12 15 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# There is variation in average birth weight across municipalities in the Northen Border States
rand(fitb2)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## mpio_birth 2345 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The sd at the municipality level is 0.27, that means a variance of 0.079, while for the residuals we get an sd of 0.99, a variance of 0.9801. The coefficients are small and only a few can be interpreted. For instance, children of women who went to prenatal care are born with higher birth weight than children of whomen who didn’t went to prenatal care.
#Lets go for the intra-class correlation coefficient
ICC<-VarCorr(fitb)$mpio_birth[1]/( VarCorr(fitb)$mpio_birth[1]+attr(VarCorr(fitb), "sc"))
#About a 6.6% of difference between municipalities.
ICC*100## [1] 6.615799plot(NULL,ylim=c(-2, 2), xlim=c(0,1),
ylab="Intercept", xlab="")
title(main="Fixed Effect Model")
for (i in 1: length(fixcoef)[1]){
abline(a=fixcoef[i], b=0, lwd=1.5, col=4)
}
# Random intercepts
rancoefs1<-ranef(fitb2)$mpio_birth+fixef(fitb2)[1]
rancoefs1<-rancoefs1*sd(births2011$wgt)+mean(births2011$wgt)
summary(rancoefs1)## (Intercept)
## Min. :-1.7889
## 1st Qu.:-0.4666
## Median :-0.3650
## Mean :-0.4069
## 3rd Qu.:-0.2889
## Max. :-0.0414plot(NULL, ylim=c(-2, 2), xlim=c(0,1),
ylab="Intercept", xlab="Prenatal care")
title(main="Random Intercept Models")
for (i in 1: length(rancoefs1[,1])){
abline(a=rancoefs1[i,1], b=.19, lwd=1., col="red")
}
Now, we include a group-level predictor, homicides in the municipality.
fitbhl<-lmer(wgt~age_m+marital_s+prenatal_care+educ_m+homz+(1|mpio_birth), data = births2011)
display(fitbhl, detail=T)## lme4::lmer(formula = wgt ~ age_m + marital_s + prenatal_care +
## educ_m + homz + (1 | mpio_birth), data = births2011)
## coef.est coef.se t value
## (Intercept) -0.44 0.07 -6.61
## age_m 0.01 0.00 19.71
## marital_sDIVORCIADA -0.13 0.05 -2.89
## marital_sSEPARADA -0.08 0.04 -2.11
## marital_sSOLTERA -0.05 0.01 -6.66
## marital_sUNIÓN LIBRE -0.01 0.01 -1.53
## marital_sVIUDA 0.09 0.06 1.51
## prenatal_careSI 0.19 0.02 12.15
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA -0.01 0.01 -1.11
## educ_mNINGUNA -0.04 0.02 -1.94
## educ_mPOSGRADO -0.19 0.04 -4.81
## educ_mPRIMARIA COMPLETA 0.01 0.01 0.74
## educ_mPRIMARIA INCOMPLETA -0.02 0.01 -1.60
## educ_mPROFESIONAL -0.07 0.01 -7.28
## educ_mSECUNDARIA COMPLETA 0.01 0.01 1.46
## educ_mSECUNDARIA INCOMPLETA -0.01 0.01 -1.22
## homz -0.05 0.08 -0.58
##
## Error terms:
## Groups Name Std.Dev.
## mpio_birth (Intercept) 0.27
## Residual 0.99
## ---
## number of obs: 166200, groups: mpio_birth, 61
## AIC = 468675, DIC = 468408.5
## deviance = 468522.7rand(fitbhl)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## mpio_birth 2341 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Lets go for the intra-class correlation coefficient
ICC2<-VarCorr(fitbhl)$mpio_birth[1]/(VarCorr(fitbhl)$mpio_birth[1]+attr(VarCorr(fitbhl), "sc"))
# About a 6.8% of difference between municipalities.
ICC2*100## [1] 6.802304We now compare both models
anova(fitb2, fitbhl)## refitting model(s) with ML (instead of REML)## Data: births2011
## Models:
## object: wgt ~ age_m + marital_s + prenatal_care + educ_m + (1 | mpio_birth)
## ..1: wgt ~ age_m + marital_s + prenatal_care + educ_m + homz + (1 |
## ..1: mpio_birth)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 18 468559 468739 -234262 468523
## ..1 19 468561 468751 -234261 468523 0.3527 1 0.5526# I think is not a significant effect of adding the homicides by municipality into the model. Random slope model
fitbrsm<-lmer(wgt~age_m+marital_s+prenatal_care+educ_m+homz+(prenatal_care|mpio_birth), data = births2011, REML = F)
display(fitbrsm, detail=T)## lme4::lmer(formula = wgt ~ age_m + marital_s + prenatal_care +
## educ_m + homz + (prenatal_care | mpio_birth), data = births2011,
## REML = F)
## coef.est coef.se t value
## (Intercept) -0.46 0.08 -5.97
## age_m 0.01 0.00 19.70
## marital_sDIVORCIADA -0.13 0.05 -2.88
## marital_sSEPARADA -0.08 0.04 -2.09
## marital_sSOLTERA -0.06 0.01 -6.72
## marital_sUNIÓN LIBRE -0.01 0.01 -1.42
## marital_sVIUDA 0.09 0.06 1.52
## prenatal_careSI 0.22 0.04 5.55
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA -0.01 0.01 -1.10
## educ_mNINGUNA -0.04 0.02 -1.95
## educ_mPOSGRADO -0.19 0.04 -4.80
## educ_mPRIMARIA COMPLETA 0.01 0.01 0.79
## educ_mPRIMARIA INCOMPLETA -0.02 0.01 -1.67
## educ_mPROFESIONAL -0.07 0.01 -7.31
## educ_mSECUNDARIA COMPLETA 0.01 0.01 1.48
## educ_mSECUNDARIA INCOMPLETA -0.01 0.01 -1.27
## homz -0.05 0.08 -0.70
##
## Error terms:
## Groups Name Std.Dev. Corr
## mpio_birth (Intercept) 0.36
## prenatal_careSI 0.18 -0.73
## Residual 0.99
## ---
## number of obs: 166200, groups: mpio_birth, 61
## AIC = 468506, DIC = 468464
## deviance = 468464.0We now compare models
anova(fitbhl, fitbrsm)## refitting model(s) with ML (instead of REML)## Data: births2011
## Models:
## object: wgt ~ age_m + marital_s + prenatal_care + educ_m + homz + (1 |
## object: mpio_birth)
## ..1: wgt ~ age_m + marital_s + prenatal_care + educ_m + homz + (prenatal_care |
## ..1: mpio_birth)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 19 468561 468751 -234261 468523
## ..1 21 468506 468716 -234232 468464 58.704 2 1.789e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# The random slope model fits betterplot(NULL, ylim=c(-3, 1), xlim=c(0,1),ylab="Intercept", xlab="Prenatal care (Yes)")
title (main="Random Slope and Intercept Model - Prenatal care by municipality")
cols<-rainbow(length(unique(births2011$mpio_birth)))
for (i in 1:dim(rancoefs2$mpio_birth)[1]){
abline(a=fixef(fitbrsm)["(Intercept)"]+rancoefs2$mpio_birth[[1]][i],
b=fixef(fitbrsm)["prenatal_careSI"]+rancoefs2$mpio_birth[ , "prenatal_careSI"][i], col=cols[i])
}
Up next, the model with a binary variable: low weight at birth coded as yes (1) or no (0).
# We observed there is a difference by municipality
anova(glm(low_wgt ~ mpio_birth, family = binomial, data = births2011), test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: low_wgt
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 166199 79737
## mpio_birth 60 974.41 166139 78763 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fitlw<-glmer(low_wgt ~ age_m + marital_s + prenatal_care + educ_m + (1 | mpio_birth), family = binomial, data = births2011, control = glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun=2e5)))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.021506 (tol =
## 0.001, component 1)fitlw<-refit(fitlw)arm::display(fitlw, detail=T)## glmer(formula = low_wgt ~ age_m + marital_s + prenatal_care +
## educ_m + (1 | mpio_birth), data = births2011, family = binomial,
## control = glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 2e+05)))
## coef.est coef.se z value
## (Intercept) -2.45 0.12 -21.01
## age_m 0.01 0.00 3.55
## marital_sDIVORCIADA 0.34 0.16 2.15
## marital_sSEPARADA 0.22 0.15 1.48
## marital_sSOLTERA 0.09 0.03 2.60
## marital_sUNIÓN LIBRE 0.06 0.02 2.53
## marital_sVIUDA -0.31 0.26 -1.20
## prenatal_careSI -0.58 0.05 -10.91
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA 0.10 0.04 2.47
## educ_mNINGUNA -0.03 0.08 -0.34
## educ_mPOSGRADO 0.21 0.14 1.46
## educ_mPRIMARIA COMPLETA 0.06 0.04 1.51
## educ_mPRIMARIA INCOMPLETA 0.08 0.06 1.40
## educ_mPROFESIONAL 0.09 0.04 2.33
## educ_mSECUNDARIA COMPLETA 0.04 0.03 1.28
## educ_mSECUNDARIA INCOMPLETA 0.12 0.04 2.66
## Pr(>|z|)
## (Intercept) 0.00
## age_m 0.00
## marital_sDIVORCIADA 0.03
## marital_sSEPARADA 0.14
## marital_sSOLTERA 0.01
## marital_sUNIÓN LIBRE 0.01
## marital_sVIUDA 0.23
## prenatal_careSI 0.00
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA 0.01
## educ_mNINGUNA 0.73
## educ_mPOSGRADO 0.14
## educ_mPRIMARIA COMPLETA 0.13
## educ_mPRIMARIA INCOMPLETA 0.16
## educ_mPROFESIONAL 0.02
## educ_mSECUNDARIA COMPLETA 0.20
## educ_mSECUNDARIA INCOMPLETA 0.01
##
## Error terms:
## Groups Name Std.Dev.
## mpio_birth (Intercept) 0.60
## Residual 1.00
## ---
## number of obs: 166200, groups: mpio_birth, 61
## AIC = 78855.4, DIC = 78448.1
## deviance = 78634.7exp(fixef(fitlw)[-1])## age_m
## 1.0062174
## marital_sDIVORCIADA
## 1.3980227
## marital_sSEPARADA
## 1.2433528
## marital_sSOLTERA
## 1.0906696
## marital_sUNIÓN LIBRE
## 1.0619755
## marital_sVIUDA
## 0.7327644
## prenatal_careSI
## 0.5580320
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA
## 1.1101784
## educ_mNINGUNA
## 0.9720934
## educ_mPOSGRADO
## 1.2315022
## educ_mPRIMARIA COMPLETA
## 1.0588891
## educ_mPRIMARIA INCOMPLETA
## 1.0823524
## educ_mPROFESIONAL
## 1.0896364
## educ_mSECUNDARIA COMPLETA
## 1.0389570
## educ_mSECUNDARIA INCOMPLETA
## 1.1229115exp(confint(fitlw, method="Wald"))## 2.5 % 97.5 %
## .sig01 NA NA
## (Intercept) 0.06877318 0.108598
## age_m 1.00278096 1.009666
## marital_sDIVORCIADA 1.02953521 1.898398
## marital_sSEPARADA 0.93105288 1.660406
## marital_sSOLTERA 1.02170161 1.164293
## marital_sUNIÓN LIBRE 1.01372001 1.112528
## marital_sVIUDA 0.44100177 1.217554
## prenatal_careSI 0.50253403 0.619659
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA 1.02174549 1.206265
## educ_mNINGUNA 0.82547858 1.144749
## educ_mPOSGRADO 0.93193173 1.627370
## educ_mPRIMARIA COMPLETA 0.98285764 1.140802
## educ_mPRIMARIA INCOMPLETA 0.96917856 1.208742
## educ_mPROFESIONAL 1.01386389 1.171072
## educ_mSECUNDARIA COMPLETA 0.97969244 1.101807
## educ_mSECUNDARIA INCOMPLETA 1.03103397 1.222976# Women living in municipalities with more homicides will have higher odds of having a low weight birth, net of the individual level factors
fitlw2<-glmer(low_wgt ~ age_m + marital_s + prenatal_care + educ_m + homz+(1 | mpio_birth), family = binomial, data = births2011, control = glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun=2e5)))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00889835 (tol =
## 0.001, component 1)fitlw2<-refit(fitlw2)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00166969 (tol =
## 0.001, component 1)arm::display(fitlw2, detail=T)## glmer(formula = low_wgt ~ age_m + marital_s + prenatal_care +
## educ_m + homz + (1 | mpio_birth), data = births2011, family = binomial,
## control = glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 2e+05)))
## coef.est coef.se z value
## (Intercept) -2.34 0.15 -15.13
## age_m 0.01 0.00 3.55
## marital_sDIVORCIADA 0.34 0.15 2.18
## marital_sSEPARADA 0.22 0.15 1.46
## marital_sSOLTERA 0.09 0.03 2.60
## marital_sUNIÓN LIBRE 0.06 0.02 2.53
## marital_sVIUDA -0.31 0.26 -1.21
## prenatal_careSI -0.58 0.05 -10.93
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA 0.10 0.04 2.47
## educ_mNINGUNA -0.03 0.08 -0.33
## educ_mPOSGRADO 0.21 0.14 1.46
## educ_mPRIMARIA COMPLETA 0.06 0.04 1.51
## educ_mPRIMARIA INCOMPLETA 0.08 0.06 1.41
## educ_mPROFESIONAL 0.09 0.04 2.34
## educ_mSECUNDARIA COMPLETA 0.04 0.03 1.28
## educ_mSECUNDARIA INCOMPLETA 0.12 0.04 2.66
## homz 0.18 0.18 1.01
## Pr(>|z|)
## (Intercept) 0.00
## age_m 0.00
## marital_sDIVORCIADA 0.03
## marital_sSEPARADA 0.15
## marital_sSOLTERA 0.01
## marital_sUNIÓN LIBRE 0.01
## marital_sVIUDA 0.23
## prenatal_careSI 0.00
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA 0.01
## educ_mNINGUNA 0.74
## educ_mPOSGRADO 0.14
## educ_mPRIMARIA COMPLETA 0.13
## educ_mPRIMARIA INCOMPLETA 0.16
## educ_mPROFESIONAL 0.02
## educ_mSECUNDARIA COMPLETA 0.20
## educ_mSECUNDARIA INCOMPLETA 0.01
## homz 0.31
##
## Error terms:
## Groups Name Std.Dev.
## mpio_birth (Intercept) 0.59
## Residual 1.00
## ---
## number of obs: 166200, groups: mpio_birth, 61
## AIC = 78856.4, DIC = 78449.7
## deviance = 78635.0exp(fixef(fitlw2)[-1])## age_m
## 1.0062113
## marital_sDIVORCIADA
## 1.3979945
## marital_sSEPARADA
## 1.2434094
## marital_sSOLTERA
## 1.0905009
## marital_sUNIÓN LIBRE
## 1.0618471
## marital_sVIUDA
## 0.7326883
## prenatal_careSI
## 0.5579736
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA
## 1.1101859
## educ_mNINGUNA
## 0.9727391
## educ_mPOSGRADO
## 1.2316231
## educ_mPRIMARIA COMPLETA
## 1.0588754
## educ_mPRIMARIA INCOMPLETA
## 1.0826722
## educ_mPROFESIONAL
## 1.0897114
## educ_mSECUNDARIA COMPLETA
## 1.0389760
## educ_mSECUNDARIA INCOMPLETA
## 1.1229748
## homz
## 1.1981729exp(confint(fitlw2, method="Wald"))## 2.5 % 97.5 %
## .sig01 NA NA
## (Intercept) 0.07106548 0.1303273
## age_m 1.00277518 1.0096592
## marital_sDIVORCIADA 1.03453314 1.8891504
## marital_sSEPARADA 0.92758886 1.6667588
## marital_sSOLTERA 1.02154442 1.1641120
## marital_sUNIÓN LIBRE 1.01358996 1.1124018
## marital_sVIUDA 0.44313421 1.2114439
## prenatal_careSI 0.50254322 0.6195180
## educ_mBACHILLERATO O PREPARATORIA INCOMPLETA 1.02179235 1.2062261
## educ_mNINGUNA 0.82580636 1.1458150
## educ_mPOSGRADO 0.93107987 1.6291787
## educ_mPRIMARIA COMPLETA 0.98290564 1.1407169
## educ_mPRIMARIA INCOMPLETA 0.96949000 1.2090678
## educ_mPROFESIONAL 1.01390411 1.1711866
## educ_mSECUNDARIA COMPLETA 0.97968548 1.1018547
## educ_mSECUNDARIA INCOMPLETA 1.03098763 1.2231693
## homz 0.84291746 1.7031541# It seems like there isn't a significative difference in low weight births whit the inclusion of a group level variable: homicides.
anova(fitlw, fitlw2)## Data: births2011
## Models:
## fitlw: low_wgt ~ age_m + marital_s + prenatal_care + educ_m + (1 | mpio_birth)
## fitlw2: low_wgt ~ age_m + marital_s + prenatal_care + educ_m + homz +
## fitlw2: (1 | mpio_birth)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fitlw 17 78855 79026 -39411 78821
## fitlw2 18 78856 79037 -39410 78820 0.9579 1 0.3277