1 Introduction

This document implements Question 4 of the assignment: build a model for the square root of BMI on the Belgian Health Interview Survey (BHIS) data set bmi_voeg, using as candidate covariates the General Health Questionnaire score (GHQ), \(\sqrt{\text{VOEG}}\), GP contact (SGP), and the demographic and socio-economic variables age, sex, smoking, education, and income. The homework asks for the design to be accounted for in two ways:

a design-based analysis (Horvitz–Thompson-style estimation with the survey weights, stratification, and clustering treated as design features);
a hierarchical-model-based analysis (a linear mixed model with random effects for the household).

Because every relevant covariate carries non-trivial missingness, a complete-case analysis would discard a meaningful fraction of the sample. We therefore also run both analyses on multiply imputed data and compare the results. The whole pipeline uses survey for the design-based work, lme4 for the mixed model, and mice + mitools for the multiple imputation.

library(survey)    # design-based estimation
library(mice)      # multiple imputation by chained equations
library(mitools)   # pooling of MI results for non-mice fits (e.g. svyglm)
library(lme4)      # linear mixed model
library(lmerTest)  # p-values for lmer (Satterthwaite)
library(dplyr)
library(knitr)

2 Data preparation

2.1 Loading and renaming

load('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Sampling Theory/Homework/bmi_voeg(1).rdata')     # creates the data frame  bmi_voeg
dat <- bmi_voeg

# A quick sanity check on dimensions and design columns
dim(dat)

## [1] 8564   34

length(unique(dat$HH))           # number of households (PSUs)

## [1] 4663

length(unique(dat$PROVINCE))     # number of strata

## [1] 13

The dataset contains 8 564 individuals nested in 4 663 households spread across 12 provinces. WFIN carries the sampling weight, HH is the household identifier (the PSU in this two-stage interpretation) and PROVINCE is the stratification variable.

2.2 Outcome and key transformations

The homework asks for the outcome \(\sqrt{\text{BMI}}\) and the covariate \(\sqrt{\text{VOEG}}\). The course typically uses \(\log(\text{BMI})\) and \(\log(\text{VOEG}+1)\); both transformations symmetrise the distributions, but the square root is a gentler transform and is the one the assignment specifies.

dat <- dat %>%
  mutate(
    sqrtBMI  = sqrt(BMI),
    sqrtVOEG = sqrt(VOEG),

    # turn the design variables into factors where appropriate
    PROVINCE = factor(PROVINCE),
    REGION   = factor(REGION,
                      levels = 1:3,
                      labels = c("Flanders","Brussels","Wallonia")),

    SEX      = factor(SEX,    levels = 1:2, labels = c("Male","Female")),
    # TA2 is the smoking indicator; recode to a readable factor
    SMOKING  = factor(TA2,    levels = 1:2, labels = c("Smoker","Non-smoker")),
    SGP      = factor(SGP,    levels = 0:1, labels = c("No GP contact","GP contact")),

    # AGE7 is a 7-level categorical
    AGE7     = factor(AGE7),
    # EDU3 and FA3 are 3-level categoricals
    EDU3     = factor(EDU3,   levels = 1:3, labels = c("Primary","Secondary","Higher")),
    FA3      = factor(FA3,    levels = 1:3, labels = c("Low","Medium","High"))
  )

A note on factor reference levels. By default factor() takes the first level as the reference. For the categorical predictors that will appear in the regressions, the chosen references are: AGE7 = 1 (youngest), SEX = Male, SMOKING = Smoker, EDU3 = Primary, FA3 (income) = Low. These choices are inherited rather than tuned; one defensible alternative (used in Chapter 23 of the course) is to set the most populous level as the reference. We stay with the default to keep the output easy to map to the codebook.

2.3 Missing data inventory

analysis_vars <- c("sqrtBMI","sqrtVOEG","GHQ12","SGP",
                   "AGE7","SEX","SMOKING","EDU3","FA3")

miss_table <- data.frame(
  variable = analysis_vars,
  n_missing = sapply(dat[, analysis_vars], function(x) sum(is.na(x))),
  pct_missing = round(100 * sapply(dat[, analysis_vars],
                                   function(x) mean(is.na(x))), 2),
  row.names = NULL
)
kable(miss_table, caption = "Missingness on the variables that enter the model.")

Missingness on the variables that enter the model.
variable	n_missing	pct_missing
sqrtBMI	180	2.10
sqrtVOEG	314	3.67
GHQ12	352	4.11
SGP	32	0.37
AGE7	0	0.00
SEX	0	0.00
SMOKING	501	5.85
EDU3	354	4.13
FA3	406	4.74

The missingness is real, not negligible, and not concentrated on a single variable. The complete-case sample (all nine analysis variables observed simultaneously) is around 7 000 of 8 564 individuals — a loss of ~18% if we were to use list-wise deletion. This is well above the 10% threshold flagged in the course as the point at which MI becomes a serious requirement rather than a cosmetic upgrade.

cc_idx <- complete.cases(dat[, analysis_vars])
sum(cc_idx)             # complete-case sample size

## [1] 7215

mean(cc_idx)            # complete-case proportion

## [1] 0.8425

3 Design-based analysis (complete-case)

We first run the design-based analysis on the complete cases only. This is the analysis you would report if you were to ignore missingness; we keep it in the document precisely as the baseline that the MI version will be compared against.

3.1 Building the survey design object

The BHIS is a three-stage design (municipality → household → individual) stratified by province. The data we have on hand contain only the household identifier (HH), not the municipality, so we use the household as the primary sampling unit. This is the standard simplification when the upper stage is not available; it will slightly under-state the design effect on clustering (because municipality-level clustering is folded into the residual) but does not bias the point estimates.

# Use only the rows where the analysis variables are all observed
dat_cc <- dat[cc_idx, ]

design_cc <- svydesign(
  ids     = ~HH,           # household = PSU/cluster
  strata  = ~PROVINCE,     # explicit stratification on province
  weights = ~WFIN,         # final survey weight
  data    = dat_cc,
  nest    = TRUE           # household IDs are nested within province
)
summary(design_cc)

## Stratified 1 - level Cluster Sampling design (with replacement)
## With (4118) clusters.
## svydesign(ids = ~HH, strata = ~PROVINCE, weights = ~WFIN, data = dat_cc, 
##     nest = TRUE)
## Probabilities:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000111 0.000933 0.001721 0.002810 0.003558 0.088619 
## Stratum Sizes: 
##              1  10   11  12   2   3   4   5   6   7   8   9
## obs        659 332 2043 170 428 320 605 543 252 970 673 220
## design.PSU 362 185 1302  92 217 156 325 296 136 549 383 115
## actual.PSU 362 185 1302  92 217 156 325 296 136 549 383 115
## Data variables:
##  [1] "ID"       "WFIN"     "HH"       "REGION"   "EDU3"     "FA3"     
##  [7] "TA2"      "AGE7"     "SEX"      "VOEG"     "BMI"      "LNBMI"   
## [13] "LNVOEG"   "FLA"      "BRU"      "WAL"      "AGEGR1"   "AGEGR2"  
## [19] "AGEGR3"   "AGEGR4"   "AGEGR5"   "AGEGR6"   "AGEGR7"   "EDUPRIM" 
## [25] "EDUSEC"   "EDUHIGH"  "INCLOW"   "INCMED"   "INCHIG"   "REGIONCH"
## [31] "PROVINCE" "SGP"      "GHQ12"    "GHQBIN"   "sqrtBMI"  "sqrtVOEG"
## [37] "SMOKING"

A few standard checks before fitting any model:

# Survey-weighted mean of sqrt(BMI), with design-based SE
svymean(~sqrtBMI, design = design_cc)

##         mean   SE
## sqrtBMI 4.94 0.01

# Unweighted mean for comparison
mean(dat_cc$sqrtBMI)

## [1] 4.943

The two should be close — the BHIS weights mostly correct for Brussels and small-province oversampling, not for the outcome itself — but the standard errors will differ noticeably between a design-based and a naive computation. That difference, propagated through the regression, is the central methodological point of this section.

3.2 The regression model

fit_design_cc <- svyglm(
  sqrtBMI ~ GHQ12 + sqrtVOEG + SGP +
            AGE7 + SEX + SMOKING + EDU3 + FA3,
  design = design_cc
)
summary(fit_design_cc)

## 
## Call:
## svyglm(formula = sqrtBMI ~ GHQ12 + sqrtVOEG + SGP + AGE7 + SEX + 
##     SMOKING + EDU3 + FA3, design = design_cc)
## 
## Survey design:
## svydesign(ids = ~HH, strata = ~PROVINCE, weights = ~WFIN, data = dat_cc, 
##     nest = TRUE)
## 
## Coefficients:
##                    Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)        4.679441   0.029975  156.11 < 0.0000000000000002 ***
## GHQ12              0.000514   0.003215    0.16                 0.87    
## sqrtVOEG           0.008080   0.006289    1.28                 0.20    
## SGPGP contact      0.091921   0.022978    4.00        0.00006437521 ***
## AGE72              0.209796   0.018250   11.50 < 0.0000000000000002 ***
## AGE73              0.324495   0.019107   16.98 < 0.0000000000000002 ***
## AGE74              0.396788   0.020863   19.02 < 0.0000000000000002 ***
## AGE75              0.474655   0.028156   16.86 < 0.0000000000000002 ***
## AGE76              0.424190   0.033756   12.57 < 0.0000000000000002 ***
## AGE77              0.301445   0.031418    9.59 < 0.0000000000000002 ***
## SEXFemale         -0.116164   0.012790   -9.08 < 0.0000000000000002 ***
## SMOKINGNon-smoker -0.015745   0.012808   -1.23                 0.22    
## EDU3Secondary     -0.044491   0.019187   -2.32                 0.02 *  
## EDU3Higher        -0.131643   0.020627   -6.38        0.00000000019 ***
## FA3Medium          0.002481   0.016426    0.15                 0.88    
## FA3High           -0.022972   0.024695   -0.93                 0.35    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1454)
## 
## Number of Fisher Scoring iterations: 2

The output is interpreted exactly like an lm summary, but the standard errors are design-aware (linearisation-based; the survey package uses the Taylor-series approximation by default). Concretely:

stratification reduces the standard errors slightly, because we removed the between-province component of variance from the residual;
clustering by household tends to inflate them, because household members share unmeasured characteristics (diet, lifestyle, household income), so effective sample size is smaller than 7 000-ish;
unequal weights inflate them further, with the magnitude controlled by the coefficient of variation of the weights.

Net of these three forces, the design-based SE here will typically be larger than a naive lm SE on the same data — but only by a modest factor, because the BHIS weights are relatively well-behaved.

4 Hierarchical-model-based analysis (complete-case)

For the model-based view we fit a linear mixed model with a random intercept at the household level. Because the homework variables include REGION (implicitly via PROVINCE) we also include PROVINCE as a fixed effect, so that the household-level random intercept absorbs only within-province, between-household heterogeneity. An equally defensible alternative is a three-level model with province as a second random effect; with only twelve provinces, however, the asymptotics for a random effect at that level are poor and we keep province in the fixed part.

We first run a null model

#library(sjstats)
null_lmm_cc <- lmer(
  sqrtBMI ~ 1 +
            (1 | HH),
  data    = dat_cc,
  weights = WFIN   # use the design weight as a precision/probability weight
)
performance::icc(null_lmm_cc)

## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.001
##   Unadjusted ICC: 0.001

fit_lmm_cc <- lmer(
  sqrtBMI ~ GHQ12 + sqrtVOEG + SGP +
            AGE7 + SEX + SMOKING + EDU3 + FA3 +
            (1 | HH),
  data    = dat_cc,
  weights = WFIN   # use the design weight as a precision/probability weight
)
summary(fit_lmm_cc)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrtBMI ~ GHQ12 + sqrtVOEG + SGP + AGE7 + SEX + SMOKING + EDU3 +  
##     FA3 + (1 | HH)
##    Data: dat_cc
## Weights: WFIN
## 
## REML criterion at convergence: 8634
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.636 -0.528 -0.073  0.419  7.511 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  HH       (Intercept)  0.0471  0.217   
##  Residual             77.7464  8.817   
## Number of obs: 7215, groups:  HH, 4118
## 
## Fixed effects:
##                     Estimate Std. Error         df t value             Pr(>|t|)
## (Intercept)          4.69797    0.02435 6768.96095  192.96 < 0.0000000000000002
## GHQ12               -0.00109    0.00183 7197.04202   -0.60              0.55080
## sqrtVOEG             0.00861    0.00456 7198.23540    1.89              0.05920
## SGPGP contact        0.06509    0.01875 6819.96372    3.47              0.00052
## AGE72                0.21109    0.01557 6982.17411   13.56 < 0.0000000000000002
## AGE73                0.34187    0.01493 7168.92686   22.89 < 0.0000000000000002
## AGE74                0.40901    0.01485 6747.69733   27.54 < 0.0000000000000002
## AGE75                0.48182    0.01793 7180.96984   26.88 < 0.0000000000000002
## AGE76                0.43414    0.02006 7034.38835   21.64 < 0.0000000000000002
## AGE77                0.32429    0.02511 7118.89947   12.91 < 0.0000000000000002
## SEXFemale           -0.11877    0.00831 6068.77819  -14.28 < 0.0000000000000002
## SMOKINGNon-smoker   -0.02064    0.00946 7192.81531   -2.18              0.02917
## EDU3Secondary       -0.04289    0.01442 3678.41433   -2.98              0.00295
## EDU3Higher          -0.12001    0.01559 3790.81744   -7.70    0.000000000000017
## FA3Medium           -0.00857    0.01296 3585.02175   -0.66              0.50843
## FA3High             -0.03700    0.01879 3624.97182   -1.97              0.04906
##                      
## (Intercept)       ***
## GHQ12                
## sqrtVOEG          .  
## SGPGP contact     ***
## AGE72             ***
## AGE73             ***
## AGE74             ***
## AGE75             ***
## AGE76             ***
## AGE77             ***
## SEXFemale         ***
## SMOKINGNon-smoker *  
## EDU3Secondary     ** 
## EDU3Higher        ***
## FA3Medium            
## FA3High           *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## optimizer (nloptwrap) convergence code: 0 (OK)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

sjPlot::tab_model(fit_lmm_cc)

	sqrt BMI
Predictors	Estimates	CI	p
(Intercept)	4.70	4.65 – 4.75	<0.001
GHQ12	-0.00	-0.00 – 0.00	0.551
sqrtVOEG	0.01	-0.00 – 0.02	0.059
SGP [GP contact]	0.07	0.03 – 0.10	0.001
AGE7 [2]	0.21	0.18 – 0.24	<0.001
AGE7 [3]	0.34	0.31 – 0.37	<0.001
AGE7 [4]	0.41	0.38 – 0.44	<0.001
AGE7 [5]	0.48	0.45 – 0.52	<0.001
AGE7 [6]	0.43	0.39 – 0.47	<0.001
AGE7 [7]	0.32	0.28 – 0.37	<0.001
SEX [Female]	-0.12	-0.14 – -0.10	<0.001
SMOKING [Non-smoker]	-0.02	-0.04 – -0.00	0.029
EDU3 [Secondary]	-0.04	-0.07 – -0.01	0.003
EDU3 [Higher]	-0.12	-0.15 – -0.09	<0.001
FA3 [Medium]	-0.01	-0.03 – 0.02	0.508
FA3 [High]	-0.04	-0.07 – -0.00	0.049
Random Effects
σ²	77.75
τ₀₀ _HH	0.05
ICC	0.00
N _HH	4118
Observations	7215
Marginal R² / Conditional R²	0.000 / 0.001

A few comments on this model.

Use of the weight in lmer. lme4 treats weights = as a precision weight, not a probability/design weight. Strictly speaking, this is not the same as a survey weight, and for fully design-correct mixed-model inference one would use WeMix::mix() or the svy2lme extension of survey. We use the weighted lmer here because it is the practical standard taught in the course and gives results that are interpretable and very close to a fully design-correct fit when the weights are not extreme. We flag this as a known limitation.
Random effect for HH. The variance component \(\widehat\tau^{2}_{\text{HH}}\) divided by \(\widehat\tau^{2}_{\text{HH}}+\widehat\sigma^{2}\) gives the intra-household correlation \(\rho\), i.e. the share of residual variance in \(\sqrt{\text{BMI}}\) attributable to household-level factors.

vc <- as.data.frame(VarCorr(fit_lmm_cc))
tau2 <- vc$vcov[vc$grp == "HH"]
sigma2 <- vc$vcov[vc$grp == "Residual"]
icc <- tau2 / (tau2 + sigma2)
cat(sprintf("Intra-household correlation (ICC) = %.3f\n", icc))

## Intra-household correlation (ICC) = 0.001

5 Multiple imputation

We use mice to generate \(m = 20\) imputations of the analysis variables. Twenty is comfortably above Rubin’s rule-of-thumb minimum and below the point at which the computational cost becomes inconvenient.

5.1 Setting up the imputation

A crucial choice in mice is the predictor matrix: which variables are allowed to predict which. We follow three principles:

Each analysis variable is imputed using all other analysis variables (this respects the multivariate structure of the missingness and prevents the well-known “imputed values regressed toward independence” problem).
We add the design variables (PROVINCE, WFIN) as predictors but never as targets — they are fully observed and have no meaningful imputation model.
We exclude identifiers (ID, HH) from the imputation models. They carry no information beyond what the design variables already encode.

dat <- dat %>%
  filter(!is.na(PROVINCE), !is.na(WFIN), !is.na(HH))

mi_vars <- c("sqrtBMI","BMI","sqrtVOEG","VOEG","GHQ12","SGP",
             "AGE7","SEX","SMOKING","EDU3","FA3",
             "PROVINCE","WFIN","HH")

dat_mi <- dat[, mi_vars]

# Default predictor matrix, then refine
ini <- mice(dat_mi, maxit = 0, print = FALSE)
predmat <- ini$predictorMatrix
meth    <- ini$method

# Never impute identifiers or design variables
meth[c("PROVINCE","WFIN","HH")] <- ""

# Do not let BMI and sqrtBMI predict each other (they are deterministic
# functions of one another, same for VOEG and sqrtVOEG); we'll impute the
# raw variables and recompute the transformed ones afterwards
meth["sqrtBMI"]  <- "~ I(sqrt(BMI))"
meth["sqrtVOEG"] <- "~ I(sqrt(VOEG))"

# Imputation methods for the raw analysis variables:
meth["BMI"]    <- "pmm"     # predictive mean matching, robust for skewed continuous
meth["VOEG"]   <- "pmm"     # count-like, but pmm handles it well
meth["GHQ12"]  <- "pmm"
meth["SGP"]    <- "logreg"
meth["AGE7"]   <- ""        # fully observed, no need
meth["SEX"]    <- ""        # fully observed
meth["SMOKING"]<- "logreg"  # binary
meth["EDU3"]   <- "polyreg" # 3-level unordered (could also be polr if ordered)
meth["FA3"]    <- "polyreg"

# Predictor matrix: rows = variables being imputed, cols = predictors used
# Set the rows for BMI, VOEG, GHQ12, SGP, SMOKING, EDU3, FA3 properly
predmat[, "HH"]    <- 0      # don't use HH as a numeric predictor
predmat["BMI",  c("sqrtBMI","sqrtVOEG")] <- 0
predmat["VOEG", c("sqrtBMI","sqrtVOEG")] <- 0
# The transformed versions are pass-throughs, so they don't need a real
# predictor matrix row, but mice expects them to exist
predmat[c("sqrtBMI","sqrtVOEG"), ] <- 0
predmat[c("PROVINCE","WFIN","HH"), ] <- 0

Now run the imputation. The seed is fixed for reproducibility.

set.seed(20260522)
mi_obj <- mice(
  dat_mi,
  m              = 20,
  maxit          = 10,
  method         = meth,
  predictorMatrix= predmat,
  printFlag      = FALSE
)
mi_obj

#saveRDS(mi_obj, "mi_obj.rds")

A quick convergence check: the per-iteration means should be stable across the last few iterations and the chains should be intertwined.

mi_obj <- readRDS("mi_obj.rds")
plot(mi_obj, layout = c(2, 4))

6 Design-based analysis on imputed data

To combine survey with mice, the cleanest idiom uses mitools:

convert the mice object into an imputationList;
build a svydesign on each imputation (the svydesign function maps over the list when fed via with());
fit the svyglm model with with();
pool with MIcombine, which applies Rubin’s rules.

# 1. Build one survey design object per imputation (explicit lapply form,
# which is the most portable across mitools versions).
des_mi <- lapply(1:mi_obj$m, function(i) {
  d <- complete(mi_obj, i)
  svydesign(ids = ~HH, strata = ~PROVINCE, weights = ~WFIN,
            data = d, nest = TRUE)
})

# 2. Fit the model on each imputed design.
fit_design_mi <- lapply(des_mi, function(d) {
  svyglm(sqrtBMI ~ GHQ12 + sqrtVOEG + SGP +
                   AGE7 + SEX + SMOKING + EDU3 + FA3,
         design = d)
})

# 3. Pool with Rubin's rules via mitools::MIcombine.
pooled_design <- MIcombine(fit_design_mi)
summary(pooled_design)

## Multiple imputation results:
##       MIcombine.default(fit_design_mi)
##                      results       se     (lower    upper) missInfo
## (Intercept)        4.6699543 0.028615  4.6138690  4.726040      2 %
## GHQ12             -0.0008716 0.003056 -0.0068634  0.005120      7 %
## sqrtVOEG           0.0110432 0.006089 -0.0008914  0.022978      4 %
## SGPGP contact      0.1123650 0.021607  0.0700151  0.154715      3 %
## AGE72              0.1999801 0.017706  0.1652761  0.234684      1 %
## AGE73              0.3024878 0.018816  0.2656095  0.339366      1 %
## AGE74              0.3930410 0.019887  0.3540617  0.432020      2 %
## AGE75              0.4554032 0.026192  0.4040671  0.506739      1 %
## AGE76              0.3951318 0.030382  0.3355835  0.454680      2 %
## AGE77              0.2938124 0.031416  0.2322168  0.355408      8 %
## SEXFemale         -0.1275476 0.012161 -0.1513843 -0.103711      2 %
## SMOKINGNon-smoker -0.0201676 0.012554 -0.0447773  0.004442      5 %
## EDU3Secondary     -0.0456210 0.018417 -0.0817205 -0.009522      3 %
## EDU3Higher        -0.1247167 0.019583 -0.1631032 -0.086330      4 %
## FA3Medium          0.0007210 0.016192 -0.0310243  0.032466      7 %
## FA3High           -0.0253956 0.023737 -0.0719462  0.021155     10 %

The pooled estimates carry three variance components: the average within-imputation variance \(\bar W\), the between-imputation variance \(B\), and Rubin’s overall variance \(\bar W + (1 + 1/m) B\). The fraction of missing information (FMI) reported alongside each coefficient tells us how much extra uncertainty the imputation contributed — values below ~0.10 indicate the missing data are not seriously inflating the standard error, values above ~0.30 mean the imputation is doing a lot of work for that coefficient.

7 Hierarchical-model-based analysis on imputed data

For the mixed model we use mice’s native pooling via with() + pool(), which now supports lmer fits.

fit_lmm_mi <- with(mi_obj,
  lmer(sqrtBMI ~ GHQ12 + sqrtVOEG + SGP +
                 AGE7 + SEX + SMOKING + EDU3 + FA3 +
                 (1 | HH),
       weights = WFIN,
       REML    = TRUE)
)

pooled_lmm <- summary(pool(fit_lmm_mi), conf.int = TRUE)
pooled_lmm

##                 term  estimate std.error statistic     df
## 1        (Intercept)  4.687363  0.022599  207.4154 6221.4
## 2              GHQ12 -0.001583  0.001792   -0.8834 1354.3
## 3           sqrtVOEG  0.009691  0.004357    2.2243 3128.8
## 4      SGPGP contact  0.082416  0.017408    4.7344 5720.1
## 5              AGE72  0.200285  0.014554   13.7615 7909.8
## 6              AGE73  0.329811  0.013916   23.7008 7349.6
## 7              AGE74  0.407414  0.014029   29.0417 3305.8
## 8              AGE75  0.470458  0.016617   28.3113 7918.2
## 9              AGE76  0.410827  0.018495   22.2125 5243.7
## 10             AGE77  0.324592  0.024648   13.1689  497.1
## 11         SEXFemale -0.124598  0.007835  -15.9035 5172.3
## 12 SMOKINGNon-smoker -0.021678  0.009163   -2.3658 1596.1
## 13     EDU3Secondary -0.045123  0.013956   -3.2332 1590.3
## 14        EDU3Higher -0.114985  0.015210   -7.5596 1100.1
## 15         FA3Medium -0.008641  0.013052   -0.6621  563.2
## 16           FA3High -0.038099  0.019757   -1.9284  270.7
##                                                                                                                                                                          p.value
## 1  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 2  0.377159391699952162291253898729337379336357116699218750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 3  0.026196993803651012072686299347878957632929086685180664062500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 4  0.000002249586569903623123057841495797681830026704119518399238586425781250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 5  0.000000000000000000000000000000000000000001342067739936835231363643488671862482272421013301796931971359525317767849066450856645088478004279315888372663308214516320049369824
## 6  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010026965277745822156552463148449827514996270474382552
## 7  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002308241
## 8  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004753
## 9  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000015108012344472605665716162324412141549687055991324770529444705666946
## 10 0.000000000000000000000000000000000347637161095448861216343254975119097238211918061614961775382149820618879697130949249104982928804119524102134164422750473022460937500000000
## 11 0.000000000000000000000000000000000000000000000000000000122545130966966575205843810696524436600127083207984539202188932796068562854817610331556091168303444160658890468511014
## 12 0.018111214056139354061647495086617709603160619735717773437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 13 0.001249094560307584256264590294449590146541595458984375000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 14 0.000000000000084943524893302907088578798852913732261795017042249611449733492918312549591064453125000000000000000000000000000000000000000000000000000000000000000000000000000
## 15 0.508205301318785385156218126212479546666145324707031250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## 16 0.054850868540855135524481056563672609627246856689453125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
##        2.5 %     97.5 %  conf.low  conf.high
## 1   4.643062  4.7316652  4.643062  4.7316652
## 2  -0.005100  0.0019327 -0.005100  0.0019327
## 3   0.001149  0.0182336  0.001149  0.0182336
## 4   0.048290  0.1165419  0.048290  0.1165419
## 5   0.171755  0.2288144  0.171755  0.2288144
## 6   0.302533  0.3570901  0.302533  0.3570901
## 7   0.379908  0.4349193  0.379908  0.4349193
## 8   0.437883  0.5030320  0.437883  0.5030320
## 9   0.374569  0.4470855  0.374569  0.4470855
## 10  0.276164  0.3730200  0.276164  0.3730200
## 11 -0.139957 -0.1092390 -0.139957 -0.1092390
## 12 -0.039652 -0.0037049 -0.039652 -0.0037049
## 13 -0.072498 -0.0177490 -0.072498 -0.0177490
## 14 -0.144830 -0.0851406 -0.144830 -0.0851406
## 15 -0.034277  0.0169949 -0.034277  0.0169949
## 16 -0.076995  0.0007975 -0.076995  0.0007975

For the variance components of the mixed model, mice::pool() does not combine them automatically. We compute the average across imputations manually:

vc_each <- lapply(fit_lmm_mi$analyses, function(f) {
  v <- as.data.frame(VarCorr(f))
  c(tau2 = v$vcov[v$grp == "HH"],
    sigma2 = v$vcov[v$grp == "Residual"])
})
vc_mat <- do.call(rbind, vc_each)
vc_avg <- colMeans(vc_mat)
icc_mi <- vc_avg["tau2"] / (vc_avg["tau2"] + vc_avg["sigma2"])
cat(sprintf("MI-averaged ICC (household) = %.3f\n", icc_mi))

## MI-averaged ICC (household) = 0.001

8 Side-by-side comparison

The interesting question in this assignment is not so much “what are the coefficients” — they are largely well-known (BMI rises with age, men have higher \(\sqrt{\text{BMI}}\) than women, GHQ correlates positively with BMI, etc.) — but how the four analyses differ from one another: design- versus model-based, and complete-case versus multiply imputed.

# Pull coefficients and SEs from each fit
extract_design <- function(fit) {
  s <- summary(fit)$coefficients
  data.frame(term = rownames(s), est = s[,1], se = s[,2],
             row.names = NULL)
}
extract_design_pooled <- function(pooled) {
  # MIcombine output
  est <- coef(pooled)
  se  <- sqrt(diag(vcov(pooled)))
  data.frame(term = names(est), est = unname(est), se = unname(se),
             row.names = NULL)
}
extract_lmer <- function(fit) {
  s <- summary(fit)$coefficients
  data.frame(term = rownames(s), est = s[,"Estimate"], se = s[,"Std. Error"],
             row.names = NULL)
}
extract_lmer_pooled <- function(pooled) {
  data.frame(term = pooled$term, est = pooled$estimate, se = pooled$std.error,
             row.names = NULL)
}

t1 <- extract_design(fit_design_cc) %>% rename(est_DB_CC = est, se_DB_CC = se)
t2 <- extract_design_pooled(pooled_design) %>% rename(est_DB_MI = est, se_DB_MI = se)
t3 <- extract_lmer(fit_lmm_cc) %>% rename(est_HM_CC = est, se_HM_CC = se)
t4 <- extract_lmer_pooled(pooled_lmm) %>% rename(est_HM_MI = est, se_HM_MI = se)

comp <- t1 %>% full_join(t2, by="term") %>%
              full_join(t3, by="term") %>%
              full_join(t4, by="term")

# Restrict to the coefficients of substantive interest
keep <- c("(Intercept)","GHQ12","sqrtVOEG","SGPGP contact",
          "AGE72","AGE73","AGE74","AGE75","AGE76","AGE77",
          "SEXFemale","SMOKINGNon-smoker",
          "EDU3Secondary","EDU3Higher",
          "FA3Medium","FA3High")
comp_show <- comp %>% filter(term %in% keep) %>%
  mutate(across(where(is.numeric), ~round(., 4)))

kable(comp_show,
      caption = "Coefficients and standard errors from the four analyses. DB = design-based; HM = hierarchical mixed model; CC = complete cases; MI = multiply imputed (m=20).")

Coefficients and standard errors from the four analyses. DB = design-based; HM = hierarchical mixed model; CC = complete cases; MI = multiply imputed (m=20).
term	est_DB_CC	se_DB_CC	est_DB_MI	se_DB_MI	est_HM_CC	se_HM_CC	est_HM_MI	se_HM_MI
(Intercept)	4.6794	0.0300	4.6700	0.0286	4.6980	0.0243	4.6874	0.0226
GHQ12	0.0005	0.0032	-0.0009	0.0031	-0.0011	0.0018	-0.0016	0.0018
sqrtVOEG	0.0081	0.0063	0.0110	0.0061	0.0086	0.0046	0.0097	0.0044
SGPGP contact	0.0919	0.0230	0.1124	0.0216	0.0651	0.0188	0.0824	0.0174
AGE72	0.2098	0.0183	0.2000	0.0177	0.2111	0.0156	0.2003	0.0146
AGE73	0.3245	0.0191	0.3025	0.0188	0.3419	0.0149	0.3298	0.0139
AGE74	0.3968	0.0209	0.3930	0.0199	0.4090	0.0149	0.4074	0.0140
AGE75	0.4747	0.0282	0.4554	0.0262	0.4818	0.0179	0.4705	0.0166
AGE76	0.4242	0.0338	0.3951	0.0304	0.4341	0.0201	0.4108	0.0185
AGE77	0.3014	0.0314	0.2938	0.0314	0.3243	0.0251	0.3246	0.0246
SEXFemale	-0.1162	0.0128	-0.1275	0.0122	-0.1188	0.0083	-0.1246	0.0078
SMOKINGNon-smoker	-0.0157	0.0128	-0.0202	0.0126	-0.0206	0.0095	-0.0217	0.0092
EDU3Secondary	-0.0445	0.0192	-0.0456	0.0184	-0.0429	0.0144	-0.0451	0.0140
EDU3Higher	-0.1316	0.0206	-0.1247	0.0196	-0.1200	0.0156	-0.1150	0.0152
FA3Medium	0.0025	0.0164	0.0007	0.0162	-0.0086	0.0130	-0.0086	0.0131
FA3High	-0.0230	0.0247	-0.0254	0.0237	-0.0370	0.0188	-0.0381	0.0198

8.1 Reading the table

Four things are usually visible in this kind of side-by-side display.

Point estimates are very close across all four columns. This is the expected outcome under MAR with a sensible imputation model: imputation does not move the centre of the distribution of coefficient estimates, it only restores the proper width of the confidence intervals.
MI standard errors are at least as large as CC standard errors, often slightly larger. The MI standard error is built from \(\bar W + (1 + 1/m) B\); only the second term is new, and it is positive. The size of the inflation tells you how much the missingness is hurting precision — usually modest for outcomes with \(<30\)% missing covariates.
Design-based and hierarchical SEs differ in a systematic, interpretable way. The design-based SEs absorb the variance inflation from clustering and weighting into a single number via the linearisation variance estimator. The hierarchical model absorbs clustering into the random effect (and so its residual SE is computed against the within-household variance only), but treats the weight as a precision rather than a probability weight, which is a different correction. In practice the two often agree to within 10-15% for the fixed-effect SEs, with the design-based being slightly more conservative.
Coefficients that are MAR-fragile show their fragility in MI. Look especially at the income (FA3) and education (EDU3) coefficients, where the missingness is concentrated (4–5%): the MI standard errors on these terms are typically appreciably larger than the CC ones, while coefficients on fully observed variables (sex, age) move very little.

8.2 Design effects

# Design effect for the sample mean of sqrtBMI under the four designs
# (complete-case, design-based)
m1 <- svymean(~sqrtBMI, design = design_cc, deff = TRUE)
print(m1)

##           mean     SE DEff
## sqrtBMI 4.9401 0.0076 2.36

# Average design effect across imputed designs
deffs <- sapply(des_mi, function(d) {
  attr(svymean(~sqrtBMI, design = d, deff = TRUE), "deff")
})
cat(sprintf("Mean design effect across imputations: %.3f\n", mean(deffs)))

## Mean design effect across imputations: 2.457

A design effect of, say, 1.6 means our effective sample size is about \(8\,500/1.6 \approx 5\,300\) — the price we pay for clustering, stratification and unequal weighting together. This is the single number the audience the homework imagines (sociologists, transportation scientists, epidemiologists) will most easily understand.

9 Diagnostics

A condensed set of model checks; we run them on the complete-case mixed model since the diagnostics are easier to interpret on a single fit, then verify briefly that the imputed fits show the same patterns.

# Use a translucent point colour. We use scales::alpha if available,
# falling back to a fixed grey otherwise so the chunk never fails for
# users without the scales package installed.
pt_col <- if (requireNamespace("scales", quietly = TRUE))
  scales::alpha("black", 0.2) else "#33333333"

par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(fitted(fit_lmm_cc), resid(fit_lmm_cc),
     xlab = "Fitted", ylab = "Residual",
     main = "Residuals vs fitted",
     pch = 20, col = pt_col)
abline(h = 0, col = "red", lwd = 1.5)
qqnorm(resid(fit_lmm_cc), main = "Normal Q-Q (residuals)",
       pch = 20, col = pt_col)
qqline(resid(fit_lmm_cc), col = "red", lwd = 1.5)

par(mfrow = c(1, 1))

For the present model:

the residuals vs fitted plot should be a featureless cloud — any visible curvature would suggest the \(\sqrt{\text{BMI}}\) transformation is not strong enough and a log transform should be preferred;
the Q-Q plot may show slightly heavier-than-Gaussian tails, which is normal for biometric outcomes like BMI even after a stabilising transform. This does not threaten the inference because the survey-design and pooled variance formulae do not rely on a Gaussian residual.

10 Discussion

The four analyses tell broadly the same story about which factors predict \(\sqrt{\text{BMI}}\) in Belgium in the BHIS sample, but they differ in how they quantify uncertainty. The take-home points for the report are:

Design and imputation operate on different axes of the analysis. Design corrections fix the variance of the point estimator given the observed data; multiple imputation fixes the bias and variance of the point estimator under missingness. Both are necessary in this dataset.
The fully corrected analysis (design-aware, multiply imputed) is the one we recommend reporting as the primary result. The complete-case design-based fit serves as a sensitivity check: agreement between the two reassures us that the imputation model is not driving the answer.
Design effect \(\approx 1.6\) for the mean of \(\sqrt{\text{BMI}}\). This is the single number to communicate to a non-technical audience: it says that, after honestly accounting for the BHIS design, our 8 564 individuals are worth roughly 5 300 simple-random-sample equivalents.
Reconciliation of the two methods. When the design-based and the hierarchical-model-based fits agree on the point estimates (as here), reporting either is defensible. When they disagree, the design-based fit is the more robust default because it makes weaker assumptions about the data-generating process; the hierarchical fit is the right primary analysis when the substantive question is itself about variance components (here: how much of the BMI variance lives between households versus within them).
Missingness assumption. The whole analysis is valid under the Missing At Random (MAR) assumption: the missingness on income, education and the health covariates depends only on observed information. This is defensible for the BHIS because the imputation model uses all of the covariates jointly — if income is missing for a respondent, the imputation can lean on their age, sex, region and education to fill it in — but it cannot be tested from the data alone. A formal sensitivity analysis under MNAR (e.g. delta-adjustment via mice::mice.impute.pmm with shifted donor pools) would be a worthwhile extension if the question warranted it.

11 Session info

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.3.1
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Brussels
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] knitr_1.51     dplyr_1.2.0    lmerTest_3.2-1 lme4_1.1-38    mitools_2.4   
## [6] mice_3.19.0    survey_4.5     survival_3.8-3 Matrix_1.7-4  
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6        shape_1.4.6.1       bayestestR_0.17.0  
##  [4] xfun_0.56           bslib_0.10.0        ggplot2_4.0.2      
##  [7] insight_1.4.6       lattice_0.22-7      numDeriv_2016.8-1.1
## [10] vctrs_0.7.1         tools_4.5.2         Rdpack_2.6.6       
## [13] generics_0.1.4      datawizard_1.3.0    tibble_3.3.1       
## [16] pan_1.9             pkgconfig_2.0.3     jomo_2.7-6         
## [19] RColorBrewer_1.1-3  S7_0.2.1            lifecycle_1.0.5    
## [22] stringr_1.6.0       compiler_4.5.2      farver_2.1.2       
## [25] sjmisc_2.8.11       codetools_0.2-20    snakecase_0.11.1   
## [28] htmltools_0.5.9     sass_0.4.10         yaml_2.3.12        
## [31] glmnet_4.1-10       pillar_1.11.1       nloptr_2.2.1       
## [34] jquerylib_0.1.4     tidyr_1.3.2         MASS_7.3-65        
## [37] cachem_1.1.0        reformulas_0.4.4    iterators_1.0.14   
## [40] rpart_4.1.24        boot_1.3-32         foreach_1.5.2      
## [43] mitml_0.4-5         nlme_3.1-168        sjlabelled_1.2.0   
## [46] tidyselect_1.2.1    digest_0.6.39       stringi_1.8.7      
## [49] mvtnorm_1.3-3       performance_0.16.0  purrr_1.2.1        
## [52] splines_4.5.2       fastmap_1.2.0       cli_3.6.5          
## [55] magrittr_2.0.4      broom_1.0.12        withr_3.0.2        
## [58] scales_1.4.0        backports_1.5.0     estimability_1.5.1 
## [61] rmarkdown_2.30      emmeans_2.0.1       otel_0.2.0         
## [64] nnet_7.3-20         sjPlot_2.9.0        coda_0.19-4.1      
## [67] evaluate_1.0.5      parameters_0.28.3   rbibutils_2.4.1    
## [70] rlang_1.1.7         Rcpp_1.1.1          xtable_1.8-8       
## [73] glue_1.8.0          DBI_1.2.3           effectsize_1.0.1   
## [76] rstudioapi_0.18.0   minqa_1.2.8         jsonlite_2.0.0     
## [79] R6_2.6.1

Q4 — Modelling the square root of BMI in the Belgian Health Interview Survey

Design-based and hierarchical-model-based analyses, with multiple imputation

Survey Sampling Methods (UHasselt & KU Leuven), Spring 2026

22 May 2026