/
s2_get_dir_est.R
103 lines (77 loc) · 2.77 KB
/
s2_get_dir_est.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
rm(list=ls())
library(data.table)
library(survey)
library(sp)
library(scales)
library(RColorBrewer)
library(gridExtra)
library(rgeos)
library(maptools)
library(latticeExtra)
library(xtable)
library(Hmisc)
options(survey.lonely.psu = "adjust")
load("all_chil_6m.RData")
##################################
# create storage tables
##################################
state_dir_est_temp <- NULL
#######
# DHS
#######
for (year in c(2003, 2008, 2013, 2018)){
b_6m_vt <- sort(unique(all_chil_dt_final[svy_year == year, ]$b_6m))
for (i in 1:length(b_6m_vt)){
b_6m_i <- b_6m_vt[i]
svydata <- all_chil_dt_final[b_6m == b_6m_i & svy_year == year, ]
# design object
svydes <- svydesign(data = svydata,
id = ~cluster+hh,
weights = ~sampwgt,
strata = ~strata,
nest = T)
# state-level mcv estimates
mcv_state_dt <- as.data.table(svyby(formula = ~mcv, by = ~StateName, design = svydes, FUN = svymean))
mcv_state_dt$b_6m <- unique(svydata$b_6m)
mcv_state_dt$t_6m <- unique(svydata$t_6m)
mcv_state_dt$svy_year <- unique(svydata$svy_year)
mcv_state_dt$svy_type <- unique(svydata$svy_type)
state_dir_est_temp <- rbind(state_dir_est_temp, mcv_state_dt)
print(b_6m_i)
}
}
################
# MICS & SMART
################
for (year in c(2007, 2011, 2014, 2015, 2016)){
b_6m_vt <- sort(unique(all_chil_dt_final[svy_year == year, ]$b_6m))
for (i in 1:length(b_6m_vt)){
b_6m_i <- b_6m_vt[i]
svydata <- all_chil_dt_final[b_6m == b_6m_i & svy_year == year, ]
# design object
svydes <- svydesign(data = svydata,
id = ~cluster+hh,
weights = ~sampwgt,
strata = ~StateName,
nest = T)
# state-level mcv estimates
mcv_state_dt <- as.data.table(svyby(formula = ~mcv, by = ~StateName, design = svydes, FUN = svymean))
mcv_state_dt$b_6m <- unique(svydata$b_6m)
mcv_state_dt$t_6m <- unique(svydata$t_6m)
mcv_state_dt$svy_year <- unique(svydata$svy_year)
mcv_state_dt$svy_type <- unique(svydata$svy_type)
state_dir_est_temp <- rbind(state_dir_est_temp, mcv_state_dt)
print(b_6m_i)
}
}
##################################
# match to sia status and StateID
##################################
state_svy_b_t_lookup_dt <- unique(all_chil_dt_final[, c("StateName", "StateID", "b_6m", "t_6m", "N_sia")])
state_dir_est_final <- merge(state_dir_est_temp, state_svy_b_t_lookup_dt,
by = c("StateName", "b_6m", "t_6m"), all.x = T)
##################################
# save results
##################################
save(state_dir_est_final,
file = "state_dir_est_6m.RData")