stonegold546

### Block 02 -- Load packages

library(ordinalCont)  # Just to check frequentist
library(splines2)
library(brms)

### Block 03 -- Load helper functions

source("helper_file_ord_generic.R")

stonegold546

library(data.table)
library(ggplot2)
library(scales)

set.seed(12345)
post_df <- data.table(
  est = rnorm(13 * 12 / 2, 0, .05),
  se = rlnorm(13 * 12 / 2, log(.09), log(.20 / .09) / qnorm(.9999))
)
post_df[, lo := qnorm(.025, est, se)]

stonegold546

library(rstan)
library(lavaan)

cmdstanr::set_cmdstan_path("~/cmdstan/")

# https://rfarouni.github.io/assets/projects/BayesianFactorAnalysis/BayesianFactorAnalysis.html
# L is N_items BY N_factor loading matrix, constrained lower-triangular
# Cov = LL' + diag(res_vars), factors constrained orthogonal

# Chose items 3, 6 and 8 as markers of factors 1, 2 and 3

stonegold546

functions {
  int sum0(vector x) {
    int len = num_elements(x);
    int res = 0;
    for (i in 1:len) if (x[i] == 0) res += 1;
    return res;
  }
}
data {
  int N;

stonegold546

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# Simulate an example
set.seed(1234)
hist(x_true <- rt(2e2, 3) * 1.5)
# Create outcome
hist(y <- scale(rnorm(length(x_true), -.075 * x_true))[, 1])

stonegold546

library(lavaan)
library(rstan)
library(cmdstanr)

set_cmdstan_path("~/cmdstan/")

source("create_poly_stan_missing.R")

HS9 <- HolzingerSwineford1939[, paste0("x", 1:9)]
# Pearson correlations

stonegold546

library(lavaan)

X <- HolzingerSwineford1939[, paste0("x", 1:9)]
colMeans(X)
apply(X, 2, sd)

# min-max scaling version
X.p <- as.data.frame(apply(X, 2, function (x) (x - min(x)) / (max(x) - min(x))))
colMeans(X.p)
apply(X.p, 2, sd)

stonegold546

progesterone <- c(1513, 2025)
placebo <- c(1459, 2013)
(tab <- as.data.frame(rbind(progesterone, placebo)))
names(tab) <- c("pass", "total")
tab$fail <- tab$total - tab$pass
tab$treatment <- c(1, 0)
chisq.test(tab[c("pass", "fail")], correct = FALSE)
# chisq.test(long.tab$treatment, long.tab$pass, correct = FALSE)
summary(glm(cbind(pass, fail) ~ treatment, binomial, tab))

stonegold546

data {
  real<lower = 0> sd_m;
  real<lower = 0> sd_m_diff;
  real<lower = 0> sd_st;
  real<lower = 0> sd_st_r;
  int<lower = 0, upper = 1> nu_choice;
  int<lower = 0> N;
  vector<lower = 0, upper = 1>[N] x;
  vector[N] y;
}

stonegold546

library(lavaan)

summary(m0 <- cfa(
  "F1 =~ x1 + x2 + x3 + x9\n F2 =~ x4 + x5 + x6\n F3 =~ x7 + x8 + x9\n x3 ~~ x5\n x2 ~~ x7\n x4 ~~ x7",
  HolzingerSwineford1939, std.lv = TRUE, meanstructure = TRUE))

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)