cdriveraus

# Multiple subjects -- correlated intercept and slope ---------------------

N <- 50

times <- seq(0,10,1)

interceptmu <- 4 #intercept / starting point

slopemu <- .3 #slope / growth rate

cdriveraus

#generate data from python script

library('reticulate')
# reticulate::install_python()

a=py_run_string("
import matplotlib.pyplot as plt
import numpy as np
from random import random
from random import seed

cdriveraus

# The current manual is here:
# https://cran.r-project.org/web/packages/ctsem/vignettes/hierarchicalmanual.pdf
# A few blog posts with examples are here: https://cdriver.netlify.app/

# Here's a bivariate panel example script:

library(ctsem)
data(AnomAuth)
longdat <- ctDeintervalise(ctWideToLong(AnomAuth,Tpoints = 5,n.manifest = 2))

cdriveraus

invlogit <- function(x) {
  exp(x)/(1+exp(x));
}

corbase <- matrix(c(1,2,0,1),2,2)
mcor <- cov2cor(corbase %*% t(corbase))
print(mcor)
corchol <- t(chol(mcor))

r <- matrix(rnorm(2000000),ncol=2)

cdriveraus

stn_md> stancode <- 'data {real y_mean;} parameters {real y;} model {y ~ normal(y_mean,1);}'

stn_md> mod <- stan_model(model_code = stancode, verbose = TRUE)

TRANSLATING MODEL '73fc79f8b1915e8208c736914c86d1a1' FROM Stan CODE TO C++ CODE NOW.
successful in parsing the Stan model '73fc79f8b1915e8208c736914c86d1a1'.
The NEXT version of Stan will not be able to pre-process your Stan program.
Please open an issue at
 https://github.com/stan-dev/stanc3/issues 
if you can share or at least describe your Stan program. This will help ensure that Stan

cdriveraus

n <- 100000
pe <- .2 #patriarchy benefit to male ability
ba <- rnorm(n) #baseline ability
m <- rbinom(n,1,.5) #male
h <- ba + pe * m + rnorm(n)#h index
e1 <- ba + pe *m + rnorm(n)#earnings assuming caused by ability
e2 <- h + rnorm(n)#earnings assuming caused by h-index

lm(e1~h+m)
lm(e2~h+m)

cdriveraus

#' ctStanFit
#'
#' Fits a ctsem model specified via \code{\link{ctModel}} with type either 'stanct' or 'standt', using Bayseian inference software
#' Stan. 
#' 
#' @param datalong long format data containing columns for subject id (numeric values, 1 to max subjects), manifest variables, 
#' any time dependent (i.e. varying within subject) predictors, 
#' and any time independent (not varying within subject) predictors.
#' @param ctstanmodel model object as generated by \code{\link{ctModel}} with type='stanct' or 'standt', for continuous or discrete time
#' models respectively.

cdriveraus

install.packages('ctsem')
library(ctsem)

#### example 1 
scode <- "
parameters {
  real y[2]; 
} 
model {
  y[1] ~ normal(0, .5);

cdriveraus

##install ctsem (If after May 2019, CRAN is recommended)
install.packages("devtools")
library(devtools)
install_github("cdriveraus/ctsem")
library(ctsem)


#specify generative model (linear sde only at present for data generation -- on the to do list)
gm <- ctModel(
  type='omx', #omx is older model type still needed for data generation

cdriveraus

n=100
a<-rnorm(n)
b<-a*.5 + rnorm(n)
bs <- scale(b)
as<-scale(a)

summary(lm(b~a))
summary(lm(a~b))
summary(lm(bs~as))
summary(lm(as~bs))