modernity.pc<-prcomp(~fridge_stat+phone+tv_stat+radio_stat+newspaper+car_stat+bike_stat+motorcycle_stat,data=kenya_dat,center=T,scale=T,retx=T)
summary(modernity.pc)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.5161 1.0924 1.0033 0.9544 0.86270 0.84464 0.78710
## Proportion of Variance 0.2873 0.1492 0.1258 0.1139 0.09303 0.08918 0.07744
## Cumulative Proportion  0.2873 0.4365 0.5623 0.6762 0.76920 0.85838 0.93582
##                            PC8
## Standard deviation     0.71654
## Proportion of Variance 0.06418
## Cumulative Proportion  1.00000
scores<-data.frame(modernity.pc$x)
scores$name<-rownames(modernity.pc$x)
kenya_dat$name<-rownames(kenya_dat)
kenya_dat<-merge(kenya_dat, scores, by.x="name", by.y="name", all.x=F)
tail(names(kenya_dat), 20)
##  [1] "childhood_residence" "previous_residence"  "mig_rurl_urb"       
##  [4] "mig_urb_rurl"        "urban_stat"          "urban_rural"        
##  [7] "x1989"               "x1993"               "x1998"              
## [10] "x2003"               "x2008"               "x2014"              
## [13] "PC1"                 "PC2"                 "PC3"                
## [16] "PC4"                 "PC5"                 "PC6"                
## [19] "PC7"                 "PC8"
modernity.pc$rotation
##                        PC1        PC2        PC3         PC4         PC5
## fridge_stat     0.47460264 -0.2691482  0.1987914 -0.02015657 -0.04553629
## phone           0.39651207 -0.2565401  0.1326930 -0.15441722 -0.17085210
## tv_stat         0.45385087  0.1217493 -0.2114413  0.25020012 -0.01972992
## radio_stat      0.25873886  0.5288544 -0.4026169 -0.03789473 -0.61358756
## newspaper       0.37026348  0.1170035 -0.4344985  0.10501272  0.71898313
## car_stat        0.43712481 -0.2217283  0.3035192 -0.13448930 -0.08862383
## bike_stat       0.09696691  0.5545907  0.3309355 -0.68782122  0.25136876
## motorcycle_stat 0.08383537  0.4436356  0.5879079  0.63992175  0.06228851
##                         PC6        PC7           PC8
## fridge_stat     -0.25597735 -0.1370325  0.7589541440
## phone            0.82054713 -0.1232549 -0.1335418724
## tv_stat         -0.26810582 -0.6684396 -0.3909060288
## radio_stat       0.04064905  0.2933739  0.1600643652
## newspaper        0.15820678  0.3177504  0.0804183633
## car_stat        -0.36930733  0.5347155 -0.4683849571
## bike_stat       -0.03250887 -0.1901748  0.0008693396
## motorcycle_stat  0.15865311  0.1004899  0.0432979708
round(cor(kenya_dat[,c("fridge_stat","phone","tv_stat","radio_stat","newspaper","car_stat","bike_stat","motorcycle_stat", "PC1", "PC2","PC3")], method = "spearman"),3)
##                 fridge_stat  phone tv_stat radio_stat newspaper car_stat
## fridge_stat           1.000  0.359   0.363      0.083     0.229    0.454
## phone                 0.359  1.000   0.236      0.094     0.200    0.312
## tv_stat               0.363  0.236   1.000      0.271     0.331    0.273
## radio_stat            0.083  0.094   0.271      1.000     0.207    0.090
## newspaper             0.229  0.200   0.331      0.207     1.000    0.192
## car_stat              0.454  0.312   0.273      0.090     0.192    1.000
## bike_stat             0.020  0.023   0.036      0.147     0.043    0.065
## motorcycle_stat       0.032  0.003   0.091      0.067     0.010    0.045
## PC1                   0.430  0.334   0.749      0.620     0.571    0.389
## PC2                  -0.252 -0.227   0.216      0.619     0.219   -0.201
## PC3                   0.173  0.112  -0.268     -0.575    -0.452    0.260
##                 bike_stat motorcycle_stat    PC1    PC2    PC3
## fridge_stat         0.020           0.032  0.430 -0.252  0.173
## phone               0.023           0.003  0.334 -0.227  0.112
## tv_stat             0.036           0.091  0.749  0.216 -0.268
## radio_stat          0.147           0.067  0.620  0.619 -0.575
## newspaper           0.043           0.010  0.571  0.219 -0.452
## car_stat            0.065           0.045  0.389 -0.201  0.260
## bike_stat           1.000           0.103  0.269  0.601  0.260
## motorcycle_stat     0.103           1.000  0.183  0.363  0.394
## PC1                 0.269           0.183  1.000  0.499 -0.372
## PC2                 0.601           0.363  0.499  1.000 -0.300
## PC3                 0.260           0.394 -0.372 -0.300  1.000
kenya_dat$pc1q<-cut(kenya_dat$PC1, breaks = quantile(kenya_dat$PC1,probs = c(0,.25,.5,.75,1) , na.rm=T))
kenya_dat$pc2q<-cut(kenya_dat$PC2, breaks = quantile(kenya_dat$PC2,probs = c(0,.25,.5,.75,1) , na.rm=T))
kenya_dat$pc3q<-cut(kenya_dat$PC3, breaks = quantile(kenya_dat$PC3,probs = c(0,.25,.5,.75,1) , na.rm=T))
#urban only model
options(survey.lonely.psu = "adjust")
des1<-svydesign(ids=~caseid,strata=~strata,weights=~perweight,data=kenya_dat[!is.na(kenya_dat$strata),],nest=T)
fit1<-svyglm(contracept~urban_stat,design=des1,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
#urban + year
options(survey.lonely.psu = "adjust")
fit11<-svyglm(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more,design=des1,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
#urban + demographic model
options(survey.lonely.psu = "adjust")
fit2<-svyglm(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more,design=des1,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
#urban + demographic + children model
options(survey.lonely.psu = "adjust")
fit3<-svyglm(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+cheb+sonsathome+dausathome+sonsdied+dausdied,design=des1,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
#urban + demographic + children + modernity model
options(survey.lonely.psu = "adjust")
fit4<-svyglm(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+cheb+sonsathome+dausathome+sonsdied+dausdied+pc1q,design=des1,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
#urban + demographic + children + modernity + migratory model
options(survey.lonely.psu = "adjust")
fit5<-svyglm(contracept~urban_stat+x1993+x1998+x2003+x2008+x2014+more_than1_union+husband_home+kid_desire+femwork+hs_or_more+cheb+sonsathome+dausathome+sonsdied+dausdied+pc1q+literacy+mig_rurl_urb+mig_urb_rurl,design=des1,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
stargazer(fit1,fit11,fit2,fit3,fit4,fit5,type="html")
Dependent variable:
contracept
(1) (2) (3) (4) (5) (6)
urban_stat 0.355*** 0.363*** 0.363*** 0.308*** 0.106** 0.049
(0.027) (0.047) (0.047) (0.047) (0.052) (0.061)
x1993
x1998 -0.957*** -0.957*** -0.927*** -0.980***
(0.049) (0.049) (0.050) (0.057)
x2003 -0.875*** -0.875*** -0.838*** -0.837*** -0.843***
(0.047) (0.047) (0.048) (0.052) (0.057)
x2008 -0.566*** -0.566*** -0.554*** -0.586*** -0.601***
(0.056) (0.056) (0.057) (0.061) (0.062)
x2014
more_than1_union -0.414*** -0.414*** -0.260*** -0.248*** -0.184*
(0.079) (0.079) (0.082) (0.093) (0.104)
husband_home 0.154*** 0.154*** 0.142*** 0.116** 0.122**
(0.045) (0.045) (0.046) (0.051) (0.059)
kid_desire -0.448*** -0.448*** -0.613*** -0.454*** -0.383***
(0.039) (0.039) (0.045) (0.050) (0.057)
femwork 0.405*** 0.405*** 0.447*** 0.309*** 0.348***
(0.041) (0.041) (0.041) (0.046) (0.052)
hs_or_more 0.601*** 0.601*** 0.499*** 0.327*** 0.217***
(0.051) (0.051) (0.052) (0.056) (0.064)
cheb -0.182*** -0.159*** -0.160***
(0.019) (0.021) (0.025)
sonsathome 0.189*** 0.226*** 0.245***
(0.024) (0.027) (0.032)
dausathome 0.144*** 0.167*** 0.167***
(0.024) (0.027) (0.031)
sonsdied -0.157*** -0.150*** -0.100
(0.050) (0.055) (0.064)
dausdied -0.121** -0.111* -0.052
(0.053) (0.059) (0.068)
pc1q(-1.02,-0.476] 0.404*** 0.229**
(0.097) (0.112)
pc1q(-0.476,0.293] 0.635*** 0.481***
(0.106) (0.123)
pc1q(0.293,7.54] 0.776*** 0.605***
(0.104) (0.121)
literacy 0.500***
(0.074)
mig_rurl_urb 0.141
(0.133)
mig_urb_rurl 0.118
(0.092)
Constant -0.834*** -0.360*** -0.360*** -0.056 -0.493*** -0.841***
(0.014) (0.058) (0.058) (0.075) (0.126) (0.154)
Observations 53,867 21,636 21,636 21,636 16,660 13,391
Log Likelihood -33,352.900 -13,069.870 -13,069.870 -12,862.230 -10,826.330 -8,928.553
Akaike Inf. Crit. 66,709.800 26,159.740 26,159.740 25,754.460 21,688.660 17,897.110
Note: p<0.1; p<0.05; p<0.01
#urban only model


fit1b<-glmer(contracept~urban_stat+(1|geo_ke1989_2014), family=binomial, kenya_dat,
                   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_dat,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
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_dat,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
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+pc1q+(1|geo_ke1989_2014), family=binomial, kenya_dat,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
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+pc1q+literacy+mig_rurl_urb+mig_urb_rurl+(1|geo_ke1989_2014), family=binomial, kenya_dat,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
stargazer(fit1b,fit11b,fit2b,fit3b,fit4b,type="html")
Dependent variable:
contracept
(1) (2) (3) (4) (5)
urban_stat 0.446*** 0.408*** 0.382*** 0.187*** 0.154***
(0.022) (0.038) (0.039) (0.043) (0.049)
x1998 -0.886*** -0.861*** -0.943***
(0.044) (0.044) (0.050)
x2003 -0.880*** -0.852*** -0.870*** -0.860***
(0.044) (0.044) (0.048) (0.051)
x2008 -0.556*** -0.541*** -0.564*** -0.572***
(0.042) (0.042) (0.046) (0.048)
more_than1_union -0.332*** -0.198*** -0.199*** -0.160**
(0.065) (0.066) (0.074) (0.080)
husband_home 0.203*** 0.187*** 0.157*** 0.180***
(0.037) (0.038) (0.042) (0.047)
kid_desire -0.414*** -0.520*** -0.353*** -0.306***
(0.031) (0.037) (0.041) (0.045)
femwork 0.426*** 0.454*** 0.315*** 0.348***
(0.033) (0.033) (0.037) (0.041)
hs_or_more 0.565*** 0.500*** 0.336*** 0.228***
(0.039) (0.040) (0.044) (0.048)
cheb -0.153*** -0.131*** -0.134***
(0.016) (0.018) (0.020)
sonsathome 0.198*** 0.223*** 0.239***
(0.020) (0.022) (0.025)
dausathome 0.132*** 0.155*** 0.157***
(0.020) (0.022) (0.024)
sonsdied -0.139*** -0.136*** -0.104**
(0.042) (0.047) (0.052)
dausdied -0.071 -0.065 -0.007
(0.045) (0.050) (0.056)
pc1q(-1.02,-0.476] 0.385*** 0.192**
(0.082) (0.092)
pc1q(-0.476,0.293] 0.624*** 0.437***
(0.088) (0.100)
pc1q(0.293,7.54] 0.701*** 0.489***
(0.087) (0.098)
literacy 0.564***
(0.063)
mig_rurl_urb 0.114
(0.083)
mig_urb_rurl 0.110
(0.083)
Constant -1.283*** -0.788** -0.618* -0.887*** -1.225***
(0.306) (0.356) (0.361) (0.278) (0.259)
Observations 53,867 21,636 21,636 16,660 13,391
Log Likelihood -31,872.690 -12,509.650 -12,372.790 -10,220.730 -8,358.986
Akaike Inf. Crit. 63,751.380 25,041.300 24,777.580 20,479.470 16,759.970
Bayesian Inf. Crit. 63,778.060 25,129.100 24,905.300 20,626.160 16,917.520
Note: p<0.1; p<0.05; p<0.01
summary(fit4b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## contracept ~ urban_stat + x1993 + x1998 + x2003 + x2008 + x2014 +  
##     more_than1_union + husband_home + kid_desire + femwork +  
##     hs_or_more + cheb + sonsathome + dausathome + sonsdied +  
##     dausdied + pc1q + literacy + mig_rurl_urb + mig_urb_rurl +  
##     (1 | geo_ke1989_2014)
##    Data: kenya_dat
## Control: 
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##  16760.0  16917.5  -8359.0  16718.0    13370 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4427 -0.8345 -0.4540  0.9201  5.5943 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  geo_ke1989_2014 (Intercept) 0.3967   0.6299  
## Number of obs: 13391, groups:  geo_ke1989_2014, 8
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -1.225161   0.258577  -4.738 2.16e-06 ***
## urban_stat          0.153565   0.049270   3.117  0.00183 ** 
## x2003              -0.859932   0.051104 -16.827  < 2e-16 ***
## x2008              -0.572104   0.048094 -11.895  < 2e-16 ***
## more_than1_union   -0.159573   0.080202  -1.990  0.04663 *  
## husband_home        0.180268   0.046616   3.867  0.00011 ***
## kid_desire         -0.305662   0.044565  -6.859 6.94e-12 ***
## femwork             0.347842   0.041230   8.437  < 2e-16 ***
## hs_or_more          0.228066   0.048346   4.717 2.39e-06 ***
## cheb               -0.134485   0.020099  -6.691 2.21e-11 ***
## sonsathome          0.238964   0.024942   9.581  < 2e-16 ***
## dausathome          0.157482   0.024328   6.473 9.60e-11 ***
## sonsdied           -0.103641   0.051981  -1.994  0.04617 *  
## dausdied           -0.007021   0.055516  -0.126  0.89936    
## pc1q(-1.02,-0.476]  0.192309   0.091615   2.099  0.03581 *  
## pc1q(-0.476,0.293]  0.436677   0.099732   4.378 1.20e-05 ***
## pc1q(0.293,7.54]    0.488530   0.098261   4.972 6.63e-07 ***
## literacy            0.563616   0.063012   8.945  < 2e-16 ***
## mig_rurl_urb        0.113935   0.082751   1.377  0.16856    
## mig_urb_rurl        0.110170   0.082890   1.329  0.18381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients