## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
library(pinsearch)
library(lavaan)
set.seed(1131)

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
# 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)

## ----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)

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

## -----------------------------------------------------------------------------
hist(resamp_dmacs)

## -----------------------------------------------------------------------------
# Estimate bias
bias <- mean(resamp_dmacs) - my_dmacs(input_coef)
# Bias corrected dMACS
my_dmacs(input_coef) - bias

## -----------------------------------------------------------------------------
# 95% Percentile CI
quantile(resamp_dmacs, c(.025, .975))

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

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

## -----------------------------------------------------------------------------
# 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")
)