kenya_dat2<-kenya_dat[,c("contracept","urban_stat","x1993","x1998","x2003","x2008","x2014","more_than1_union","husband_home","kid_desire","femwork","hs_or_more","cheb","sonsathome","dausathome","sonsdied","dausdied","literacy","mig_rurl_urb","mig_urb_rurl","geo_ke1989_2014","fridge_stat","phone","tv_stat","radio_stat","newspaper","car_stat","bike_stat","motorcycle_stat")]
imp<-mice(data=kenya_dat2,seed = 1234)
## 
##  iter imp variable
##   1   1  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   1   2  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   1   3  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   1   4  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   1   5  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   2   1  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   2   2  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   2   3  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   2   4  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   2   5  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   3   1  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   3   2  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   3   3  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   3   4  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   3   5  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   4   1  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   4   2  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   4   3  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   4   4  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   4   5  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   5   1  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   5   2  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   5   3  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   5   4  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
##   5   5  more_than1_union  husband_home  kid_desire  femwork  hs_or_more  literacy  fridge_stat  phone  tv_stat  radio_stat  newspaper  car_stat  bike_stat  motorcycle_stat
kenya_dat3<-complete(imp)
modernity.pc<-prcomp(~fridge_stat+phone+tv_stat+radio_stat+newspaper+car_stat+bike_stat+motorcycle_stat,data=kenya_dat3,center=T,scale=T,retx=T)
summary(modernity.pc)
## Importance of components:
##                          PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.510 1.0876 1.0044 0.9553 0.86248 0.84431 0.80126
## Proportion of Variance 0.285 0.1479 0.1261 0.1141 0.09298 0.08911 0.08025
## Cumulative Proportion  0.285 0.4329 0.5590 0.6731 0.76605 0.85516 0.93541
##                            PC8
## Standard deviation     0.71885
## Proportion of Variance 0.06459
## Cumulative Proportion  1.00000
scores<-data.frame(modernity.pc$x)
scores$name<-rownames(modernity.pc$x)
kenya_dat3$name<-rownames(kenya_dat3)
kenya_dat3<-merge(kenya_dat3, scores, by.x="name", by.y="name", all.x=F)
tail(names(kenya_dat3), 20)
##  [1] "literacy"        "mig_rurl_urb"    "mig_urb_rurl"   
##  [4] "geo_ke1989_2014" "fridge_stat"     "phone"          
##  [7] "tv_stat"         "radio_stat"      "newspaper"      
## [10] "car_stat"        "bike_stat"       "motorcycle_stat"
## [13] "PC1"             "PC2"             "PC3"            
## [16] "PC4"             "PC5"             "PC6"            
## [19] "PC7"             "PC8"
modernity.pc$rotation
##                        PC1        PC2        PC3         PC4         PC5
## fridge_stat     0.47826164 -0.2769190  0.1822003 -0.02739979  0.05198816
## phone           0.39696863 -0.2816384  0.1443698 -0.13292768  0.19396320
## tv_stat         0.45239770  0.1126950 -0.1944750  0.25526395  0.04266343
## radio_stat      0.26241427  0.5251392 -0.4085084 -0.02816912  0.60751642
## newspaper       0.36600403  0.1347650 -0.4372899  0.10193297 -0.73134246
## car_stat        0.43252876 -0.2095035  0.2781477 -0.16168380  0.03515345
## bike_stat       0.10149554  0.5558230  0.3370797 -0.68587313 -0.22156582
## motorcycle_stat 0.09389016  0.4295717  0.5992074  0.63927488 -0.05991801
##                         PC6         PC7          PC8
## fridge_stat      0.20196988 -0.16926356 -0.767130520
## phone           -0.80311165 -0.04511198  0.199845420
## tv_stat          0.25539445 -0.67031430  0.404089348
## radio_stat      -0.03261830  0.30496486 -0.156688500
## newspaper       -0.17451178  0.27199307 -0.083488254
## car_stat         0.44629856  0.53876711  0.418128196
## bike_stat       -0.01267849 -0.21783722 -0.003096312
## motorcycle_stat -0.13465530  0.12363916 -0.043840874
#urban only model


fit1b<-glmer(contracept~urban_stat+(1|geo_ke1989_2014), family=binomial, kenya_dat3,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
fit11b<-glmer(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+(1|geo_ke1989_2014), family=binomial, kenya_dat3,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00420155 (tol =
## 0.001, component 1)
fit2b<-glmer(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+cheb+sonsathome+dausathome+sonsdied+dausdied+(1|geo_ke1989_2014), family=binomial, kenya_dat3,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00523655 (tol =
## 0.001, component 1)
fit3b<-glmer(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+cheb+sonsathome+dausathome+sonsdied+dausdied+PC1+PC2+PC3+(1|geo_ke1989_2014), family=binomial, kenya_dat3,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00518432 (tol =
## 0.001, component 1)
fit4b<-glmer(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+cheb+sonsathome+dausathome+sonsdied+dausdied+PC1+PC2+PC3+literacy+mig_rurl_urb+mig_urb_rurl+(1|geo_ke1989_2014), family=binomial, kenya_dat3,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00629405 (tol =
## 0.001, component 1)
stargazer(fit1b,fit11b,fit2b,fit3b,fit4b,type="html")
Dependent variable:
contracept
(1) (2) (3) (4) (5)
urban_stat 0.533*** 0.318*** 0.440*** 0.372*** 0.349***
(0.020) (0.022) (0.023) (0.024) (0.025)
x1993 0.232*** 0.369*** 0.353*** 0.346***
(0.046) (0.047) (0.047) (0.048)
x1998 0.301*** 0.500*** 0.496*** 0.496***
(0.045) (0.046) (0.046) (0.047)
x2003 0.230*** 0.430*** 0.429*** 0.430***
(0.045) (0.046) (0.046) (0.047)
x2008 0.514*** 0.719*** 0.685*** 0.677***
(0.044) (0.045) (0.045) (0.046)
x2014 0.975*** 1.155*** 1.189*** 1.204***
(0.038) (0.039) (0.040) (0.040)
more_than1_union -0.262*** -0.175*** -0.137*** -0.108***
(0.038) (0.038) (0.038) (0.039)
husband_home 0.154*** 0.115*** 0.081*** 0.105***
(0.021) (0.022) (0.022) (0.022)
kid_desire -0.931*** -0.572*** -0.560*** -0.490***
(0.019) (0.022) (0.022) (0.023)
femwork 0.688*** 0.645*** 0.624*** 0.580***
(0.019) (0.020) (0.020) (0.020)
hs_or_more 0.334*** 0.433*** 0.287*** 0.244***
(0.023) (0.024) (0.025) (0.026)
cheb 0.017** 0.028*** 0.069***
(0.009) (0.009) (0.009)
sonsathome 0.212*** 0.217*** 0.223***
(0.011) (0.011) (0.012)
dausathome 0.161*** 0.167*** 0.162***
(0.011) (0.011) (0.011)
sonsdied -0.215*** -0.211*** -0.223***
(0.025) (0.025) (0.025)
dausdied -0.188*** -0.185*** -0.187***
(0.027) (0.027) (0.027)
PC1 0.109*** 0.082***
(0.007) (0.007)
PC2 0.167*** 0.115***
(0.009) (0.009)
PC3 -0.078*** -0.048***
(0.009) (0.009)
literacy 0.906***
(0.031)
mig_rurl_urb 0.152***
(0.047)
mig_urb_rurl 0.263***
(0.050)
Constant -1.453*** -2.048*** -2.854*** -2.829*** -3.711***
(0.277) (0.333) (0.402) (0.365) (0.326)
Observations 70,289 70,289 70,289 70,289 70,289
Log Likelihood -40,142.900 -36,601.590 -35,813.480 -35,496.350 -35,023.030
Akaike Inf. Crit. 80,291.800 73,229.180 71,662.960 71,034.700 70,094.050
Bayesian Inf. Crit. 80,319.280 73,348.260 71,827.850 71,227.070 70,313.900
Note: p<0.1; p<0.05; p<0.01