refitME_tutorial.Rmd

---
title: '`refitME`: Measurement Error Modelling using Monte Carlo Expectation Maximization in `R`'
author: "Jakub Stoklosa, Hwang W-H., & David Warton"
date: "`r format(Sys.time(), '%d %B, %Y')`"
header-includes: # Allows you to add in your own Latex packages.
- \usepackage{float} # Use the 'float' package.
- \floatplacement{figure}{H} # Make every figure with caption = h
- \usepackage{ragged2e}
- \usepackage{geometry}
- \geometry{verbose,tmargin=3cm,bmargin=3cm,lmargin=3cm,rmargin=3cm}
output:
  pdf_document: default
  fig_caption: yes
  html_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.pos = 'H')
```

This tutorial documents fitting an MCEM algorithm via the `refitME` `R`-package. For more specific details see: *`refitME`: Measurement Error Modelling using Monte Carlo Expectation Maximization in `R`.* Also, see ```?refitME``` for further details on the fitting function, input arguments and output.

## Example 1: A simple GLM example taken from Carroll _et al._ (2006). 

The Framingham heart study data set. We also fit `simex` models and compare them with MCEM. Computational times for both models are also reported.

### Load data and `R`-packages.

```{r, echo=TRUE}
suppressWarnings(suppressMessages(library(refitME)))
suppressWarnings(suppressMessages(library(simex)))

set.seed(2020)

B <- 100  # The number of Monte Carlo replication values/SIMEX simulations.
data(Framinghamdata)
```

### Setup all variables (the construction below follows the Carroll et al. (2006) monograph).

```{r, echo=TRUE}
W <- as.matrix(Framinghamdata$w1) # Matrix of error-contaminated covariate.
sigma.sq.u <- 0.01259/2 # ME variance, obtained from Carroll et al. (2006) monograph.
```

### Fit the naive model.

The first stored variable `w1` is the error contaminated variable used in the analysis.

```{r, echo=TRUE}
mod_naiv1 <- glm(Y ~ w1 + z1 + z2 + z3, x = TRUE, family = binomial, 
                 data = Framinghamdata)
```

\clearpage

### Fit the SIMEX model.

```{r, echo=TRUE}
start <- Sys.time()
mod_simex1 <- simex(mod_naiv1, SIMEXvariable = c("w1"),
                  measurement.error = cbind(sqrt(sigma.sq.u)), B = B) # SIMEX.
end <- Sys.time()
t1 <- difftime(end,start, units = "secs")
comp.time <- c(t1)
```

### Fit the MCEM model.

```{r, echo=TRUE}
start <- Sys.time()
est <- refitME(mod_naiv1, sigma.sq.u, W, B)
end <- Sys.time()
t2 <- difftime(end, start, units = "secs")
comp.time <- c(comp.time, t2)
```

### Report model estimates and compare computational times.

```{r, echo=TRUE}
est.beta <- rbind(coef(mod_naiv1), coef(mod_simex1), est$beta)
est.beta.se <- rbind(sqrt(diag(vcov(mod_naiv1))),
                   sqrt(diag(mod_simex1$variance.jackknife)), est$beta.se2)
row.names(est.beta) = row.names(est.beta.se) <- c("Naive GLM", "SIMEX", "MCEM")
colnames(est.beta) = colnames(est.beta.se) <- c("(Intercept)", "SBP", "chol. level", 
                                                "age", "smoke")
round(est.beta, digits = 3)
round(est.beta.se, digits = 3)  # Standard error estimates.

names(comp.time) <- c("SIMEX", "MCEM")
comp.time  # SIMEX and MCEM.
```

## Example 2: A GAM example taken from Ganguli _et al._ (2005).

The Milan mortality air pollution data set. Here, we fit GAM models via the `mgcv` package where one covariate is error-contaminated.

### Load data and `R`-packages.

```{r, echo=TRUE}
suppressWarnings(suppressMessages(library(refitME)))
suppressWarnings(suppressMessages(library(SemiPar)))

set.seed(2020)

epsilon <- 0.00001  # A set convergence threshold.
B <- 50  # The number of Monte Carlo replication values.

family <- "poisson"

data(milan.mort)
dat.air <- milan.mort
``` 

### Setup all variables.

```{r, echo=TRUE}
Y <- dat.air[, 6]  # Mortality counts.
n <- length(Y)

z1 <- (dat.air[, 1])
z2 <- (dat.air[, 4])
z3 <- (dat.air[, 5])
w1 <- log(dat.air[, 9])
W <- as.matrix(w1)
dat <- data.frame(cbind(Y, z1, z2, z3, w1))

sigma.sq.u <- 0.0915 # This gives a reliability ratio of 0.7.
rel.rat <- round((1 - sigma.sq.u/var(dat$w1))*100, digits = 0)
```

\clearpage

### Fit the naive model.

```{r, echo=TRUE}
mod_naiv1 <- gam(Y ~ s(w1) + s(z1, k = 25) + s(z2) + s(z3), family = "poisson", data = dat)
```

```{r, echo=FALSE, fig.show='hide'}
plot_mod_naiv_P <- plot(mod_naiv1, select = 1)
```

### Fit the MCEM model.

```{r, echo=TRUE}
est <- refitME(mod_naiv1, sigma.sq.u, W, B, se.comp = FALSE)
```

```{r fig1, fig.height=6, fig.pos='H', echo=FALSE, fig.cap="\\label{fig:fig1} *Plots of smooths against covariate. TSP (top left is the error contaminated variable).*"}
xlab.names <- c("log(TSP)", "Day", "Temp", "Humidity")

op <- par(mfrow = c(2, 2), las = 1)

for(i in 1:4) {
  if (i == 1) {
    plot(est$mod, select = i, ylim = c(-0.35, 0.1), xlim = range(plot_mod_naiv_P[[1]]$x), 
         rug = FALSE, col = "blue", all.terms = TRUE,
         xlab = xlab.names[i], ylab = "s(Mortaity counts)", lwd = 2, cex.lab = 1.3, cex.axis = 1.3,
         cex.main = 2, font.lab = 1.1, cex = 1.4, shade = T)
    lines(plot_mod_naiv_P[[1]]$x, plot_mod_naiv_P[[1]]$fit, type = "l", col = "red", lwd = 2, lty = 2)
    title(main = bquote("Reliability ratio of predictor is"~.(rel.rat) ~ "%"),
          outer = F, line = 1, cex = 1.4)
    legend("bottomright", c("Naive GAM", "MCEM GAM"), col = c("red", "blue"), 
           lty = c(2, 1), lwd = 2, bty = "n")
    for(j in 1:2) {
      axis(j, labels = FALSE)
    }
  }
  if (i == 2) {
    plot(est$mod, select = i, ylim = c(-0.25, 0.3), rug = FALSE, col = "blue", all.terms = TRUE,
         xlab = xlab.names[i], ylab = "s(Mortaity counts)", lwd = 2, cex.lab = 1.3, cex.axis = 1.3,
         cex.main = 2, font.lab = 1.1, cex = 1.4, shade = T)
    lines(plot_mod_naiv_P[[2]]$x, plot_mod_naiv_P[[2]]$fit, type = "l", col = "red", lwd = 2, lty = 2)
    for(j in 1:2) {
      axis(j, labels = FALSE)
    }
  }
  if (i == 3) {
    plot(est$mod, select = i, ylim = c(-0.2, 0.4), rug = FALSE, col = "blue", all.terms = TRUE,
         xlab = xlab.names[i], ylab = "s(Mortaity counts)", lwd = 2, cex.lab = 1.3, cex.axis = 1.3,
         cex.main = 2, font.lab = 1.1, cex = 1.4, shade = T)
    lines(plot_mod_naiv_P[[3]]$x, plot_mod_naiv_P[[3]]$fit, type = "l", col = "red", lwd = 2, lty = 2)
    for(j in 1:2) {
      axis(j, labels = FALSE)
    }
  }
  if (i == 4) {
    plot(est$mod, select = i, ylim = c(-0.06, 0.08), rug = FALSE, col = "blue", all.terms = TRUE,
         xlab = xlab.names[i], ylab = "s(Mortaity counts)", lwd = 2, cex.lab = 1.1, cex.axis = 1.1,
         cex.main = 2, font.lab = 1.1, cex = 1.4, shade = T)
    lines(plot_mod_naiv_P[[4]]$x, plot_mod_naiv_P[[4]]$fit, type = "l", col = "red", lwd = 2, lty = 2)
    for(j in 1:2) {
      axis(j, labels = FALSE)
    }
  }
}

title(main = "MCEM (Poisson GAM) fitted to the air pollution data.", outer = T, line = -2)
par(op)
```