Solid cooking fuel, kitchen arrangement and child haemoglobin

# Solid cooking fuel, kitchen arrangement and child haemoglobin
# Data: NFHS-5 (2019-21), children aged 6-59 months
 
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(survey)

Loading required package: grid
Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

Loading required package: survival

Attaching package: 'survey'

The following object is masked from 'package:graphics':

    dotchart

library(broom) 
 
setwd("~/Desktop/Statistics Research")
nfhs <- read_csv("nfhs5_final_analysis.csv")

Rows: 181378 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): case_id, mother_edu, wealth, region, state
dbl (10): psu, strata, weight, child_hb, solid_fuel, no_sep_kitchen, child_a...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(nfhs)

Rows: 181,378
Columns: 15
$ case_id          <chr> "0100101399-02_1", "0100101357-02_1", "0100101395-02_…
$ psu              <dbl> 113, 113, 113, 113, 113, 130, 130, 106, 106, 106, 131…
$ strata           <dbl> 121, 121, 121, 121, 121, 121, 121, 121, 121, 121, 121…
$ weight           <dbl> 0.193444, 0.193444, 0.193444, 0.193444, 0.193444, 0.2…
$ child_hb         <dbl> 8.3, 7.8, 9.2, 7.2, 7.2, 12.6, 10.2, 6.6, 12.9, 7.1, …
$ solid_fuel       <dbl> 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,…
$ no_sep_kitchen   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ child_age_months <dbl> 51, 12, 25, 13, 39, 57, 11, 36, 6, 12, 52, 31, 42, 32…
$ child_male       <dbl> 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0,…
$ mother_hb        <dbl> 13.8, 7.9, 8.4, 10.2, 10.2, 8.6, 10.9, 11.1, 7.7, 11.…
$ mother_edu       <chr> "Secondary", "Primary", "Secondary", "Secondary", "Se…
$ wealth           <chr> "Middle", "Poorest", "Poorest", "Middle", "Poorer", "…
$ rural            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ region           <chr> "North", "North", "North", "North", "North", "North",…
$ state            <chr> "Jammu & Kashmir", "Jammu & Kashmir", "Jammu & Kashmi…

# For each category variable, the first level is the base group.
 
nfhs <- nfhs %>%
  mutate(
    mother_edu = factor(mother_edu, levels = c("No Education", "Primary", "Secondary", "Higher")),
    wealth     = factor(wealth, levels = c("Poorest", "Poorer", "Middle", "Richer", "Richest")),
    state      = fct_relevel(state, "Uttar Pradesh")
  )
 
 
nrow(nfhs)                 # number of children

[1] 181378

n_distinct(nfhs$psu)       # number of clusters (PSUs)

[1] 29192

n_distinct(nfhs$state)     # number of states and UTs

[1] 36

sum(is.na(nfhs))           # missing values (should be 0)

[1] 0

# Table 1: averages and shares for the whole sample
 
nfhs %>%
  summarise(
    avg_child_hb       = weighted.mean(child_hb, weight),
    avg_mother_hb      = weighted.mean(mother_hb, weight),
    avg_child_age      = weighted.mean(child_age_months, weight),
    share_solid_fuel   = weighted.mean(solid_fuel, weight),
    share_no_kitchen   = weighted.mean(no_sep_kitchen, weight),
    share_male         = weighted.mean(child_male, weight),
    share_rural        = weighted.mean(rural, weight),
    share_anaemic      = weighted.mean(child_hb < 11, weight)
  )

# A tibble: 1 × 8
  avg_child_hb avg_mother_hb avg_child_age share_solid_fuel share_no_kitchen
         <dbl>         <dbl>         <dbl>            <dbl>            <dbl>
1         10.2          11.4          32.7            0.497            0.272
# ℹ 3 more variables: share_male <dbl>, share_rural <dbl>, share_anaemic <dbl>

# Shares of mother's education and wealth
nfhs %>%
  group_by(mother_edu) %>%
  summarise(total_weight = sum(weight)) %>%
  mutate(share = total_weight / sum(total_weight))

# A tibble: 4 × 3
  mother_edu   total_weight share
  <fct>               <dbl> <dbl>
1 No Education       37841. 0.214
2 Primary            22004. 0.124
3 Secondary          90494. 0.511
4 Higher             26644. 0.151

nfhs %>%
  group_by(wealth) %>%
  summarise(total_weight = sum(weight)) %>%
  mutate(share = total_weight / sum(total_weight))

# A tibble: 5 × 3
  wealth  total_weight share
  <fct>          <dbl> <dbl>
1 Poorest       43437. 0.245
2 Poorer        38699. 0.219
3 Middle        34962. 0.198
4 Richer        32736. 0.185
5 Richest       27149. 0.153

# Table 2: average haemoglobin for each fuel and kitchen group
group_means <- nfhs %>%
  group_by(solid_fuel, no_sep_kitchen) %>%
  summarise(
    mean_hb  = weighted.mean(child_hb, weight),
    children = n(),
    .groups  = "drop"
  )
 
group_means

# A tibble: 4 × 4
  solid_fuel no_sep_kitchen mean_hb children
       <dbl>          <dbl>   <dbl>    <int>
1          0              0    10.4    64563
2          0              1    10.2    17531
3          1              0    10.2    70078
4          1              1    10.1    29206

# Weights make the results represent India.
 
options(survey.lonely.psu = "adjust")
 
nfhs_design <- svydesign(
  ids     = ~psu,
  strata  = ~strata,
  weights = ~weight,
  data    = nfhs,
  nest    = TRUE
)
 
# solid_fuel * no_sep_kitchen adds both variables AND their interaction.
 
# Model 1: key variables only
model1 <- svyglm(child_hb ~ solid_fuel * no_sep_kitchen,
                 design = nfhs_design)
 
# Model 2: full model
model2 <- svyglm(child_hb ~ solid_fuel * no_sep_kitchen +
                   child_age_months + child_male + mother_hb +
                   mother_edu + wealth + rural + state,
                 design = nfhs_design)
 
# Model 3: full model without mother's haemoglobin
model3 <- svyglm(child_hb ~ solid_fuel * no_sep_kitchen +
                   child_age_months + child_male +
                   mother_edu + wealth + rural + state,
                 design = nfhs_design)
 
summary(model1)

Call:
svyglm(formula = child_hb ~ solid_fuel * no_sep_kitchen, design = nfhs_design)

Survey design:
svydesign(ids = ~psu, strata = ~strata, weights = ~weight, data = nfhs, 
    nest = TRUE)

Coefficients:
                          Estimate Std. Error  t value Pr(>|t|)    
(Intercept)               10.35044    0.01001 1034.204  < 2e-16 ***
solid_fuel                -0.15325    0.01276  -12.007  < 2e-16 ***
no_sep_kitchen            -0.18926    0.01889  -10.018  < 2e-16 ***
solid_fuel:no_sep_kitchen  0.09104    0.02339    3.893 9.92e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 2.144477)

Number of Fisher Scoring iterations: 2

summary(model2)

Call:
svyglm(formula = child_hb ~ solid_fuel * no_sep_kitchen + child_age_months + 
    child_male + mother_hb + mother_edu + wealth + rural + state, 
    design = nfhs_design)

Survey design:
svydesign(ids = ~psu, strata = ~strata, weights = ~weight, data = nfhs, 
    nest = TRUE)

Coefficients:
                                            Estimate Std. Error t value
(Intercept)                                7.5272719  0.0491153 153.257
solid_fuel                                -0.0015454  0.0145386  -0.106
no_sep_kitchen                            -0.0560542  0.0187565  -2.989
child_age_months                           0.0238646  0.0003071  77.707
child_male                                -0.0193933  0.0093464  -2.075
mother_hb                                  0.1544616  0.0033346  46.321
mother_eduPrimary                          0.0487163  0.0163677   2.976
mother_eduSecondary                        0.1397601  0.0131261  10.647
mother_eduHigher                           0.2568627  0.0190380  13.492
wealthPoorer                               0.0613118  0.0133103   4.606
wealthMiddle                               0.1051029  0.0167318   6.282
wealthRicher                               0.1461251  0.0205528   7.110
wealthRichest                              0.2013727  0.0260342   7.735
rural                                      0.0076062  0.0173079   0.439
stateAndaman & Nicobar Islands             1.0962437  0.1030383  10.639
stateAndhra Pradesh                        0.0131602  0.0406275   0.324
stateArunachal Pradesh                     0.2773082  0.0396204   6.999
stateAssam                                 0.0978661  0.0258261   3.789
stateBihar                                 0.0697595  0.0210457   3.315
stateChandigarh                            0.1908083  0.1147149   1.663
stateChhattisgarh                          0.0250314  0.0334761   0.748
stateDadra & Nagar Haveli and Daman & Diu -0.2718680  0.0694763  -3.913
stateDelhi                                -0.3206215  0.0511797  -6.265
stateGoa                                   0.3036305  0.0934082   3.251
stateGujarat                              -0.4613912  0.0283339 -16.284
stateHaryana                              -0.2465872  0.0296438  -8.318
stateHimachal Pradesh                      0.2401958  0.0449396   5.345
stateJammu & Kashmir                      -0.3346756  0.0463082  -7.227
stateJharkhand                             0.1262345  0.0269535   4.683
stateKarnataka                            -0.0956000  0.0319726  -2.990
stateKerala                                0.6344728  0.0452591  14.019
stateLadakh                               -2.1608826  0.1724348 -12.532
stateLakshadweep                           0.5420168  0.1156889   4.685
stateMadhya Pradesh                       -0.2431574  0.0251871  -9.654
stateMaharashtra                          -0.2083738  0.0319642  -6.519
stateManipur                               0.6722089  0.0359445  18.701
stateMeghalaya                             0.9554525  0.0479735  19.916
stateMizoram                               0.4996603  0.0550791   9.072
stateNagaland                              0.6870914  0.0424772  16.176
stateOdisha                                0.1902481  0.0248724   7.649
statePuducherry                           -0.1575223  0.1179758  -1.335
statePunjab                               -0.3723440  0.0407837  -9.130
stateRajasthan                            -0.1182181  0.0245149  -4.822
stateSikkim                                0.3253355  0.1165731   2.791
stateTamil Nadu                            0.1776185  0.0323605   5.489
stateTelangana                            -0.2547429  0.0372539  -6.838
stateTripura                               0.2309524  0.0428403   5.391
stateUttarakhand                           0.1098088  0.0539244   2.036
stateWest Bengal                           0.1486782  0.0283863   5.238
solid_fuel:no_sep_kitchen                  0.0159152  0.0227349   0.700
                                          Pr(>|t|)    
(Intercept)                                < 2e-16 ***
solid_fuel                                0.915349    
no_sep_kitchen                            0.002806 ** 
child_age_months                           < 2e-16 ***
child_male                                0.038001 *  
mother_hb                                  < 2e-16 ***
mother_eduPrimary                         0.002919 ** 
mother_eduSecondary                        < 2e-16 ***
mother_eduHigher                           < 2e-16 ***
wealthPoorer                              4.12e-06 ***
wealthMiddle                              3.40e-10 ***
wealthRicher                              1.19e-12 ***
wealthRichest                             1.07e-14 ***
rural                                     0.660328    
stateAndaman & Nicobar Islands             < 2e-16 ***
stateAndhra Pradesh                       0.745998    
stateArunachal Pradesh                    2.64e-12 ***
stateAssam                                0.000151 ***
stateBihar                                0.000919 ***
stateChandigarh                           0.096259 .  
stateChhattisgarh                         0.454623    
stateDadra & Nagar Haveli and Daman & Diu 9.13e-05 ***
stateDelhi                                3.79e-10 ***
stateGoa                                  0.001153 ** 
stateGujarat                               < 2e-16 ***
stateHaryana                               < 2e-16 ***
stateHimachal Pradesh                     9.12e-08 ***
stateJammu & Kashmir                      5.07e-13 ***
stateJharkhand                            2.84e-06 ***
stateKarnataka                            0.002792 ** 
stateKerala                                < 2e-16 ***
stateLadakh                                < 2e-16 ***
stateLakshadweep                          2.81e-06 ***
stateMadhya Pradesh                        < 2e-16 ***
stateMaharashtra                          7.21e-11 ***
stateManipur                               < 2e-16 ***
stateMeghalaya                             < 2e-16 ***
stateMizoram                               < 2e-16 ***
stateNagaland                              < 2e-16 ***
stateOdisha                               2.10e-14 ***
statePuducherry                           0.181820    
statePunjab                                < 2e-16 ***
stateRajasthan                            1.43e-06 ***
stateSikkim                               0.005261 ** 
stateTamil Nadu                           4.08e-08 ***
stateTelangana                            8.20e-12 ***
stateTripura                              7.07e-08 ***
stateUttarakhand                          0.041726 *  
stateWest Bengal                          1.64e-07 ***
solid_fuel:no_sep_kitchen                 0.483912    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.893536)

Number of Fisher Scoring iterations: 2

summary(model3)

Call:
svyglm(formula = child_hb ~ solid_fuel * no_sep_kitchen + child_age_months + 
    child_male + mother_edu + wealth + rural + state, design = nfhs_design)

Survey design:
svydesign(ids = ~psu, strata = ~strata, weights = ~weight, data = nfhs, 
    nest = TRUE)

Coefficients:
                                            Estimate Std. Error t value
(Intercept)                                9.2910394  0.0306043 303.586
solid_fuel                                -0.0097490  0.0148946  -0.655
no_sep_kitchen                            -0.0640576  0.0191110  -3.352
child_age_months                           0.0238340  0.0003073  77.559
child_male                                -0.0150726  0.0094492  -1.595
mother_eduPrimary                          0.0569422  0.0166393   3.422
mother_eduSecondary                        0.1515601  0.0133590  11.345
mother_eduHigher                           0.2891668  0.0194777  14.846
wealthPoorer                               0.0717795  0.0136419   5.262
wealthMiddle                               0.1249779  0.0171317   7.295
wealthRicher                               0.1728360  0.0208583   8.286
wealthRichest                              0.2419157  0.0263681   9.175
rural                                     -0.0014313  0.0178941  -0.080
stateAndaman & Nicobar Islands             1.0788089  0.1062372  10.155
stateAndhra Pradesh                       -0.0216895  0.0417120  -0.520
stateArunachal Pradesh                     0.3700009  0.0409644   9.032
stateAssam                                 0.0636287  0.0265099   2.400
stateBihar                                 0.0299299  0.0218698   1.369
stateChandigarh                            0.1324284  0.1136197   1.166
stateChhattisgarh                         -0.0055339  0.0344067  -0.161
stateDadra & Nagar Haveli and Daman & Diu -0.2909983  0.0750120  -3.879
stateDelhi                                -0.3508770  0.0517096  -6.786
stateGoa                                   0.3364944  0.0968261   3.475
stateGujarat                              -0.5300333  0.0293305 -18.071
stateHaryana                              -0.3279927  0.0308858 -10.620
stateHimachal Pradesh                      0.2433460  0.0474767   5.126
stateJammu & Kashmir                      -0.3791055  0.0469984  -8.066
stateJharkhand                             0.0705443  0.0277346   2.544
stateKarnataka                            -0.0767790  0.0324848  -2.364
stateKerala                                0.7082643  0.0469707  15.079
stateLadakh                               -2.4948854  0.1871678 -13.330
stateLakshadweep                           0.6500383  0.1192402   5.452
stateMadhya Pradesh                       -0.2604206  0.0258048 -10.092
stateMaharashtra                          -0.2258888  0.0332223  -6.799
stateManipur                               0.7831658  0.0379085  20.659
stateMeghalaya                             0.9422649  0.0494844  19.042
stateMizoram                               0.5917729  0.0559121  10.584
stateNagaland                              0.8185539  0.0437651  18.703
stateOdisha                                0.1554481  0.0257667   6.033
statePuducherry                           -0.1708422  0.1191902  -1.433
statePunjab                               -0.4266680  0.0416332 -10.248
stateRajasthan                            -0.1259528  0.0254364  -4.952
stateSikkim                                0.3745751  0.1177845   3.180
stateTamil Nadu                            0.1754380  0.0329438   5.325
stateTelangana                            -0.2855896  0.0376443  -7.587
stateTripura                               0.2017855  0.0440669   4.579
stateUttarakhand                           0.1473208  0.0575186   2.561
stateWest Bengal                           0.0866062  0.0290666   2.980
solid_fuel:no_sep_kitchen                  0.0198895  0.0231979   0.857
                                          Pr(>|t|)    
(Intercept)                                < 2e-16 ***
solid_fuel                                0.512774    
no_sep_kitchen                            0.000804 ***
child_age_months                           < 2e-16 ***
child_male                                0.110700    
mother_eduPrimary                         0.000622 ***
mother_eduSecondary                        < 2e-16 ***
mother_eduHigher                           < 2e-16 ***
wealthPoorer                              1.44e-07 ***
wealthMiddle                              3.07e-13 ***
wealthRicher                               < 2e-16 ***
wealthRichest                              < 2e-16 ***
rural                                     0.936248    
stateAndaman & Nicobar Islands             < 2e-16 ***
stateAndhra Pradesh                       0.603081    
stateArunachal Pradesh                     < 2e-16 ***
stateAssam                                0.016394 *  
stateBihar                                0.171152    
stateChandigarh                           0.243811    
stateChhattisgarh                         0.872222    
stateDadra & Nagar Haveli and Daman & Diu 0.000105 ***
stateDelhi                                1.18e-11 ***
stateGoa                                  0.000511 ***
stateGujarat                               < 2e-16 ***
stateHaryana                               < 2e-16 ***
stateHimachal Pradesh                     2.99e-07 ***
stateJammu & Kashmir                      7.55e-16 ***
stateJharkhand                            0.010979 *  
stateKarnataka                            0.018109 *  
stateKerala                                < 2e-16 ***
stateLadakh                                < 2e-16 ***
stateLakshadweep                          5.04e-08 ***
stateMadhya Pradesh                        < 2e-16 ***
stateMaharashtra                          1.07e-11 ***
stateManipur                               < 2e-16 ***
stateMeghalaya                             < 2e-16 ***
stateMizoram                               < 2e-16 ***
stateNagaland                              < 2e-16 ***
stateOdisha                               1.63e-09 ***
statePuducherry                           0.151767    
statePunjab                                < 2e-16 ***
stateRajasthan                            7.40e-07 ***
stateSikkim                               0.001474 ** 
stateTamil Nadu                           1.02e-07 ***
stateTelangana                            3.39e-14 ***
stateTripura                              4.69e-06 ***
stateUttarakhand                          0.010434 *  
stateWest Bengal                          0.002889 ** 
solid_fuel:no_sep_kitchen                 0.391241    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.948082)

Number of Fisher Scoring iterations: 2

# 95% confidence intervals for the full model
confint(model2)

                                                 2.5 %       97.5 %
(Intercept)                                7.431003337  7.623540433
solid_fuel                                -0.030041842  0.026951079
no_sep_kitchen                            -0.092818010 -0.019290439
child_age_months                           0.023262654  0.024466564
child_male                                -0.037712741 -0.001073838
mother_hb                                  0.147925701  0.160997546
mother_eduPrimary                          0.016634825  0.080797871
mother_eduSecondary                        0.114032146  0.165488028
mother_eduHigher                           0.219547242  0.294178059
wealthPoorer                               0.035222848  0.087400848
wealthMiddle                               0.072307571  0.137898201
wealthRicher                               0.105840552  0.186409607
wealthRichest                              0.150344289  0.252401088
rural                                     -0.026318178  0.041530582
stateAndaman & Nicobar Islands             0.894283155  1.298204303
stateAndhra Pradesh                       -0.066471811  0.092792257
stateArunachal Pradesh                     0.199650047  0.354966408
stateAssam                                 0.047245473  0.148486680
stateBihar                                 0.028508700  0.111010281
stateChandigarh                           -0.034038962  0.415655584
stateChhattisgarh                         -0.040583532  0.090646422
stateDadra & Nagar Haveli and Daman & Diu -0.408045162 -0.135690746
stateDelhi                                -0.420936388 -0.220306586
stateGoa                                   0.120545508  0.486715525
stateGujarat                              -0.516927172 -0.405855301
stateHaryana                              -0.304690615 -0.188483754
stateHimachal Pradesh                      0.152111872  0.328279815
stateJammu & Kashmir                      -0.425442160 -0.243909004
stateJharkhand                             0.073404263  0.179064772
stateKarnataka                            -0.158267913 -0.032931996
stateKerala                                0.545762459  0.723183144
stateLadakh                               -2.498864077 -1.822901093
stateLakshadweep                           0.315260453  0.768773203
stateMadhya Pradesh                       -0.292525445 -0.193789305
stateMaharashtra                          -0.271025318 -0.145722341
stateManipur                               0.601755802  0.742662048
stateMeghalaya                             0.861421941  1.049483112
stateMizoram                               0.391702331  0.607618361
stateNagaland                              0.603833932  0.770348950
stateOdisha                                0.141496793  0.238999381
statePuducherry                           -0.388761233  0.073716643
statePunjab                               -0.452282190 -0.292405786
stateRajasthan                            -0.166268729 -0.070167561
stateSikkim                                0.096846027  0.553825030
stateTamil Nadu                            0.114190307  0.241046762
stateTelangana                            -0.327762539 -0.181723297
stateTripura                               0.146983052  0.314921826
stateUttarakhand                           0.004114024  0.215503596
stateWest Bengal                           0.093039470  0.204316985
solid_fuel:no_sep_kitchen                 -0.028646361  0.060476753

#Effect of solid fuel with no separate kitchen 
# This is b1 + b3 (solid fuel + interaction).
 
svycontrast(model1, c("solid_fuel" = 1, "solid_fuel:no_sep_kitchen" = 1))

          contrast     SE
contrast -0.062212 0.0205

svycontrast(model2, c("solid_fuel" = 1, "solid_fuel:no_sep_kitchen" = 1))

         contrast     SE
contrast  0.01437 0.0204

svycontrast(model3, c("solid_fuel" = 1, "solid_fuel:no_sep_kitchen" = 1))

         contrast     SE
contrast  0.01014 0.0208

# Step A: Model 1 + child's age, sex and mother's haemoglobin
modelA <- svyglm(child_hb ~ solid_fuel * no_sep_kitchen +
                   child_age_months + child_male + mother_hb,
                 design = nfhs_design)
 
# Step B: Step A + mother's education, wealth and rural
modelB <- svyglm(child_hb ~ solid_fuel * no_sep_kitchen +
                   child_age_months + child_male + mother_hb +
                   mother_edu + wealth + rural,
                 design = nfhs_design)
 
# Put the three key coefficients from every step in one table
key_terms <- c("solid_fuel", "no_sep_kitchen", "solid_fuel:no_sep_kitchen")
 
steps <- bind_rows(
  tidy(model1, conf.int = TRUE) %>% mutate(step = "1. No controls"),
  tidy(modelA, conf.int = TRUE) %>% mutate(step = "2. + age, sex, mother's Hb"),
  tidy(modelB, conf.int = TRUE) %>% mutate(step = "3. + education, wealth, rural"),
  tidy(model2, conf.int = TRUE) %>% mutate(step = "4. + state (Model 2)")
) %>%
  filter(term %in% key_terms)
 
steps %>%
  select(step, term, estimate, std.error, p.value)

# A tibble: 12 × 5
   step                          term                estimate std.error  p.value
   <chr>                         <chr>                  <dbl>     <dbl>    <dbl>
 1 1. No controls                solid_fuel          -0.153      0.0128 3.99e-33
 2 1. No controls                no_sep_kitchen      -0.189      0.0189 1.39e-23
 3 1. No controls                solid_fuel:no_sep_…  0.0910     0.0234 9.92e- 5
 4 2. + age, sex, mother's Hb    solid_fuel          -0.0957     0.0124 9.79e-15
 5 2. + age, sex, mother's Hb    no_sep_kitchen      -0.162      0.0182 5.54e-19
 6 2. + age, sex, mother's Hb    solid_fuel:no_sep_…  0.0632     0.0225 5.05e- 3
 7 3. + education, wealth, rural solid_fuel           0.00510    0.0146 7.27e- 1
 8 3. + education, wealth, rural no_sep_kitchen      -0.0849     0.0186 5.28e- 6
 9 3. + education, wealth, rural solid_fuel:no_sep_…  0.0198     0.0227 3.84e- 1
10 4. + state (Model 2)          solid_fuel          -0.00155    0.0145 9.15e- 1
11 4. + state (Model 2)          no_sep_kitchen      -0.0561     0.0188 2.81e- 3
12 4. + state (Model 2)          solid_fuel:no_sep_…  0.0159     0.0227 4.84e- 1

# Are the state dummies jointly significant?
regTermTest(model2, ~state)

Wald test for state
 in svyglm(formula = child_hb ~ solid_fuel * no_sep_kitchen + child_age_months + 
    child_male + mother_hb + mother_edu + wealth + rural + state, 
    design = nfhs_design)
F =  83.25544  on  35  and  26434  df: p= < 2.22e-16 

# R-squared for each model (from the same model with lm and weights)
summary(lm(child_hb ~ solid_fuel * no_sep_kitchen,
           data = nfhs, weights = weight))$r.squared

[1] 0.004288322

summary(lm(child_hb ~ solid_fuel * no_sep_kitchen +
             child_age_months + child_male + mother_hb +
             mother_edu + wealth + rural + state,
           data = nfhs, weights = weight))$r.squared

[1] 0.120804

summary(lm(child_hb ~ solid_fuel * no_sep_kitchen +
             child_age_months + child_male +
             mother_edu + wealth + rural + state,
           data = nfhs, weights = weight))$r.squared

[1] 0.0954775

# Chart 1 - average haemoglobin by group -------
 
chart1_data <- group_means %>%
  mutate(
    fuel    = if_else(solid_fuel == 1, "Solid fuel", "Clean fuel"),
    kitchen = if_else(no_sep_kitchen == 1, "No separate kitchen", "Separate kitchen"),
    kitchen = factor(kitchen, levels = c("Separate kitchen", "No separate kitchen"))
  )
 
chart1 <- ggplot(chart1_data, aes(x = kitchen, y = mean_hb, colour = fuel, group = fuel)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_text(aes(label = format(round(mean_hb, 2), nsmall = 2)), vjust = -1, show.legend = FALSE) +
  scale_colour_manual(values = c("Clean fuel" = "#2f6690", "Solid fuel" = "#b3541e")) +
  scale_y_continuous(limits = c(10.05, 10.40)) +
  labs(
    x = NULL,
    y = "Average haemoglobin (g/dL)",
    colour = NULL,
    caption = "Survey-weighted averages, NFHS-5 (2019-21). The y-axis does not start at zero."
  ) +
  theme_minimal() +
  theme(legend.position = "top")
 
chart1

ggsave("chart1_group_means.png", chart1, width = 7, height = 4.5, dpi = 300, bg = "white")

chart2_data <- steps %>%
  mutate(
    term = recode(term,
                  "solid_fuel" = "Solid fuel",
                  "no_sep_kitchen" = "No separate kitchen",
                  "solid_fuel:no_sep_kitchen" = "Interaction"),
    term = factor(term, levels = c("Solid fuel", "No separate kitchen", "Interaction")),
    step = fct_rev(step)
  )
 
chart2 <- ggplot(chart2_data, aes(x = estimate, y = step)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "black") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0.2, colour = "pink") +
  geom_point(size = 2.5, colour = "#9b1c31") +
  facet_wrap(~term) +
  labs(
    x = "Change in child haemoglobin (g/dL), with 95% confidence interval",
    y = NULL,
    caption = "Survey-weighted OLS, NFHS-5 (2019-21), N = 181,378"
  ) +
  theme_minimal()
 
chart2

ggsave("chart2_coefficients.png", chart2, width = 9, height = 4.2, dpi = 300, bg = "white")