---
title: "Confidence Intervals for dMACS"
author: "Mark Lai"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Confidence Intervals for dMACS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r}
library(pinsearch)
library(lavaan)
set.seed(1131)
```

## Effect Size for Noninvariance

```{r}
mod <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
# Output the final partial invariance model, and the noninvariant items
ps1 <-
    pinSearch(mod,
        data = HolzingerSwineford1939,
        group = "school", type = "residuals",
        effect_size = TRUE
    )
ps1
```

## Confidence Intervals with Resampling

### Monte Carlo Methods

```{r}
# Extract the parameter estimates for dmacs input
lav_free <- lavInspect(ps1[[1]])
# Index of the input parameters (item 3)
input_ind <- c(
    lav_free[[1]]$lambda[3, 1],
    lav_free[[1]]$nu[3, 1],
    lav_free[[2]]$nu[3, 1],
    lav_free[[1]]$theta[3, 3]
)
input_coef <- coef(ps1[[1]])[input_ind]
# Wrapper for computing dmacs with the input
my_dmacs <- function(x) {
    dmacs(
        intercepts = matrix(x[2:3]),
        loadings = matrix(x[1], nrow = 2),
        uniqueness = matrix(x[4], nrow = 2),
        ns = c(156, 145)
    )
}
# Note the effect size is slightly different
# using implied SD
my_dmacs(input_coef)
```

```{r, eval=!file.exists("resamp_dmacs.rds"), message=FALSE}
# Resample coefficients
num_resamp <- 10000
resamp_coef <- MASS::mvrnorm(
    n = num_resamp,
    mu = input_coef,
    Sigma = vcov(ps1[[1]])[input_ind, input_ind]
)
# Resampled distribution of dmacs
resamp_dmacs <- apply(resamp_coef, 1, my_dmacs)
```

```{r, include=FALSE}
if (file.exists("resamp_dmacs.rds")) {
    resamp_dmacs <- readRDS("resamp_dmacs.rds")
} else {
    saveRDS(resamp_dmacs, file = "resamp_dmacs.rds")
}
```

```{r}
hist(resamp_dmacs)
```

Bias correction

```{r}
# Estimate bias
bias <- mean(resamp_dmacs) - my_dmacs(input_coef)
# Bias corrected dMACS
my_dmacs(input_coef) - bias
```

Confidence Intervals

```{r}
# 95% Percentile CI
quantile(resamp_dmacs, c(.025, .975))
```

### Nonparametric Bootstrap

```{r, eval=!file.exists("boot_dmacs.rds")}
boot_dmacs <- bootstrapLavaan(ps1[[1]], R = 1000, FUN = pin_effsize)
```

```{r, include=FALSE}
if (file.exists("boot_dmacs.rds")) {
    boot_dmacs <- readRDS("boot_dmacs.rds")
} else {
    saveRDS(boot_dmacs, file = "boot_dmacs.rds")
}
```

```{r}
# Bias
boot_bias <- colMeans(boot_dmacs, na.rm = TRUE) - ps1$effect_size
# Bias corrected dMACS
ps1$effect_size - boot_bias
# 95% Percentile CI
boot::boot.ci(
    structure(list(t0 = ps1$effect_size,
                   t = boot_dmacs,
                   R = 1000, class = "boot")),
    type = c("norm", "basic", "perc")
)
```