Missing data frequently occurs in real-world datasets, 
leading to biased estimates and reduced statistical power.

## Primary types of missingness

- Missing completely at random (MCAR): 
The likelihood of missingness is unrelated to any observed or unobserved data. Examples:
  - A participant accidentally skips a survey question due to a technical glitch.
  - A random subset of data is lost due to equipment malfunction during data collection.
- Missing at random (MAR): The probability of missingness is related to observed data but not to the missing values themselves. Examples:
  - Participants with lower educational backgrounds might not report their income due to misunderstanding the survey question.
  - Older respondents might skip answering online survey questions due to difficulty using the technology.
- Missing not at random (MNAR): The missingness is related directly to the missing data itself, creating more challenging conditions for imputation. Examples:
  - Individuals with high incomes may choose not to disclose their exact income in a survey.
  - Patients experiencing severe health issues may systematically refuse to answer questions about their health condition.

## Commonly used methods

### Complete-case analysis

- Remove records with missing values. 

- Easy but may result in substantial data loss and biased results if data aren't MCAR.

### Simple imputation (with mean/median/mode)

- Replace missing numeric values with the mean or median of the observed values.

- Replace missing categorical values with the mode (the most frequent value).

- Easy but can lead to biased results and underestimated variance.

- May produce unbiased point estimates under MCAR

### Reweighting

- More commonly adopted in analyzing survey data

- Procedure:
  1. Confirm auxiliary variables that 
  are observed for both respondents and nonrespondents,
  related to the likelihood of response,
  and related to the survey variables of interest
  1. Estimate the probability of response for each unit in the sample 
  using a logistic regression model with auxiliary variables as predictors
  1. Focus on complete cases
  and take the reciprocal of the probability of response as the weight
  for each record
  
- Ineffective if missingness is MNAR 
(reweighting can't adjust for unobserved factors)

### Iterative regression imputation (single imputation)

- Predict missing values using regression, iteratively updating until stable.

- Better than the mean/median/mode imputation, but doesn't fully account for uncertainty.

- Can work for MCAR and MAR if the regression model is correctly specified

```{r, eval = F}
# Illustrative example: impute the airquality data
## Method I: Frequentist
### Initial mean imputation
df <- airquality
df$Ozone[is.na(df$Ozone)] <- mean(df$Ozone, na.rm = TRUE)
df$Solar.R[is.na(df$Solar.R)] <- mean(df$Solar.R, na.rm = TRUE)

### Set convergence criteria
tolerance <- 1e-4
max_iterations <- 20

### Iterative process with convergence check
prev_df <- df
for (i in 1:max_iterations) {
  #### Impute Ozone
  model_ozone <- lm(Ozone ~ Solar.R + Wind + Temp + as.factor(Month) + as.factor(Day), data = df)
  missing_ozone <- which(is.na(airquality$Ozone))
  if (length(missing_ozone) > 0) {
    df$Ozone[missing_ozone] <- predict(model_ozone, newdata = df[missing_ozone, ])
  }
  
  #### Impute Solar.R
  model_solar <- lm(Solar.R ~ Ozone + Wind + Temp + as.factor(Month) + as.factor(Day), data = df)
  missing_solar <- which(is.na(airquality$Solar.R))
  if (length(missing_solar) > 0) {
    df$Solar.R[missing_solar] <- predict(model_solar, newdata = df[missing_solar, ])
  }
  
  #### Check the mean squared difference
  diff <- mean((df$Ozone - prev_df$Ozone)^2 + (df$Solar.R - prev_df$Solar.R)^2)
  cat("Iteration:", i, "Difference:", diff, "\n")
  
  if (diff < tolerance) {
    cat("Convergence achieved after", i, "iterations.\n")
    break
  }
  
  #### Update previous dataframe for next iteration
  prev_df <- df
}
```

```{r, eval = F}
# Illustrative example: impute the airquality data
## Method II: Bayes
### Default regression models include:
### predictive mean matching (pmm, for numeric data) 
### logistic regression (logreg, for factors with 2 levels) 
### polytomous regression (polyreg, for unordered factor > 2 levels) 
### proportional odds model for (polr, ordered factors > 2 levels)
library(mice)
colSums(is.na(airquality))
imp_obj <- mice(
  airquality,
  m = 1, # single imputation
  maxit = 20, 
  method = c("pmm","pmm","", "", "", ""), 
  seed = 123
)
### Extract the imputed data
data_imp <- complete(imp_obj, 1)
plot(imp_obj) # Each trajectory should stabilize; no upward/downward trends = good convergence

## Could be problematic if using the convergence rule for frequentist
maxit <- 20
### Initialize the imputation model
ini <- mice(airquality, m = 1, maxit = 0, seed = 123)
prev_imp <- ini
prev_df = complete(prev_imp, 1)
### One iteration at a time
for (i in 1:max_iterations) {
  imp <- mice.mids(prev_imp, maxit = 1, printFlag = FALSE)  # perform one iteration
  df <- complete(imp, 1) # extract current imputed dataset
  
  diff <- mean((df$Ozone - prev_df$Ozone)^2 + (df$Solar.R - prev_df$Solar.R)^2)
  cat("Iteration:", i, "Difference:", diff, "\n")
  if (diff < 1e-4) {
    cat("Convergence achieved after", i, "iterations.\n")
    break
  }
  prev_df <- df
  prev_imp = imp
}
```

### Multiple imputation (MI)

- Designed to produce unbiased and efficient estimates under MCAR and MAR, 
provided the imputation model is correctly specified.

- Procedure:
  1. Create multiple imputed datasets.
  1. Analyze each separately.
  1. Pool the results using the Rubin's rule [@Rubin1987].

- Rubin's rule

$$
  \text{Pooled point estimate: }
  \bar Q = m^{-1}\sum_{j=1}^m Q_j,\quad\text{ where } Q_j\text{ is a point estimate}
$$
$$
  \text{Average within-imputation variance: } V_w = m^{-1}\sum_{j=1}^m{\rm var}(Q_j)
$$
$$
  \text{Between-imputation variance: }
  V_b = \frac{1}{m-1}\sum_{j=1}^m(Q_j-\bar Q)^2
$$
$$
  \text{Total variance: }V_t = V_w+(1+m^{-1})V_b
$$
$$
  \text{Degrees of freedom: }\nu =(m-1)[1+ V_w/\{(1+m^{-1})V_b\}]^2
$$
$$
  \text{95\% confidence interval: }
  \bar Q \pm t_{\nu, 0.975}\sqrt{V_t}
$$

```{r, eval = F}
# Illustrative example: impute the airquality data
library(mice)
colSums(is.na(airquality))
imp_obj <- mice(
  airquality,
  m = 2, maxit = 10, 
  method = c("pmm","pmm","", "", "", ""), 
  seed = 123
)
## Check convergence visually
## Each trajectory should stabilize and mix well over iterations
## No upward/downward trends = good convergence.
## Wiggly lines or divergence = possible convergence issues.
plot(imp_obj)
## Extract imputed data
for (i in 1:imp_obj$m) {
  data_name <- paste0("data_imp_", i)
  assign(data_name, complete(imp_obj, i))
}
fit <- with(data = imp_obj, exp = lm(Ozone ~ Solar.R + Wind + Temp))
pooled <- pool(fit, rule = "rubin1987") # Pool the results using Rubin's rules
summary(pooled, conf.int = T) # Summarize the pooled result
```

```{r, eval = F}
# Illustrative example: impute the mice::nhanes2 data
library(mice)
colSums(is.na(nhanes2))
imp_obj <- mice(
  nhanes2, m = 2, maxit = 10, 
  method = c("", "pmm", "logreg", "pmm"), 
  seed = 123
)
plot(imp_obj)
for (i in 1:imp_obj$m) {
  data_name <- paste0("data_imp_", i)
  assign(data_name, complete(imp_obj, i))
}
fit <- with(data = imp_obj, exp = lm(chl ~ age + bmi + hyp))
pooled <- pool(fit, rule = "rubin1987")
summary(pooled, conf.int = T)
```

### Maximum likelihood estimation

- Incorporates the density of missingness into the likelihood when estimating unknown parameters.
Sophistication is needed in constructing and optimizing complex likelihood functions.

- Required for MNAR.

## Bibliography
