inst/doc/mc-pls.R

In plssem: Complex Partial Least Squares Structural Equation Modeling

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  eval = identical(Sys.getenv("IN_PKGDOWN"), "true")
)

library(plssem)

## -----------------------------------------------------------------------------
# set.seed(7208094)
# n <- 10000
# r <- 0.6
# x <- rnorm(n)
# y <- r * x + rnorm(n, sd = sqrt(1 - r^2))
# cor(x, y)

## -----------------------------------------------------------------------------
# z <- as.integer(x - mean(x) > 0)
# w <- as.integer(y - mean(y) > 0)

## -----------------------------------------------------------------------------
# cor(z, w)

## -----------------------------------------------------------------------------
# g <- function(r_hat, R = 5000) {
#   # Generate continuous data
#   x <- rnorm(R)
#   y <- r_hat * x + rnorm(R, sd = sqrt(1 - r_hat^2))
# 
#   # Generate binary data
#   z <- as.integer(x - mean(x) > 0)
#   w <- as.integer(y - mean(y) > 0)
# 
#   # Return a data.frame
#   data.frame(z, w)
# }
# 
# f <- function(df) {
#   # Biased correlation between binary data
#   cor(df$z, df$w)
# }

## -----------------------------------------------------------------------------
# print(f(g(0.0)))
# print(f(g(0.3)))
# print(f(g(0.5)))
# print(f(g(0.6)))

## -----------------------------------------------------------------------------
# r_biased_observed <- cor(z, w)
# h <- function(r_hat) {
#   f(g(r_hat)) - r_biased_observed
# }

## -----------------------------------------------------------------------------
# r_hat <- r_biased_observed # reasonable starting value
# iter <- 20
# history <- r_hat
# 
# for (i in seq_len(iter)) {
#   temperature <- 1 / sqrt(i)
#   epsilon <- h(r_hat)
# 
#   cat(sprintf("Iteration: %2d, r_hat: %.3f, epsilon: %6.3f\n", i, r_hat, epsilon))
# 
#   r_hat <- r_hat - temperature * epsilon
#   history <- c(history, r_hat)
# }

## ----echo=FALSE---------------------------------------------------------------
# t <- seq_along(history)
# 
# fit <- stats::nls(
#   history ~ c + a * exp(-k * t),
#   start = list(
#     c = mean(utils::tail(history, 3)),
#     a = history[1] - mean(utils::tail(history, 3)),
#     k = 0.1
#   ),
#   algorithm = "port",
#   lower = c(c = -Inf, a = -Inf, k = 0)
# )
# 
# fpredict <- function(t) {
#   stats::predict(fit, newdata = data.frame(t = t))
# }
# 
# df <- data.frame(t = t, r_hat = history)
# 
# ggplot2::ggplot(df, ggplot2::aes(x = t, y = r_hat)) +
#   ggplot2::geom_hline(
#     yintercept = r,
#     linetype = 2,
#     linewidth = 0.8,
#     color = "grey40"
#   ) +
#   ggplot2::geom_function(
#     fun = fpredict,
#     linewidth = 1.2,
#     color = "#0072B2"
#   ) +
#   ggplot2::geom_point(size = 2, color = "#D55E00") +
#   ggplot2::labs(
#     x = "Iteration",
#     y = "r_hat",
#     title = "Convergence (with exponential fit)"
#   ) +
#   ggplot2::theme_minimal()

## -----------------------------------------------------------------------------
# set.seed(2346259)
# pls('y ~ x', data = data.frame(x = z, y = w),
#     ordered = c("y", "x"), mcpls = TRUE)