Skip to content

Instantly share code, notes, and snippets.

@cdriveraus
cdriveraus / regressionci.R
Created October 30, 2016 19:39
Does a 95% confidence interval contain the true parameter 95% of the time?
b=.8
n=100
nruns=10000
out<-rep(NA,nruns)
for(i in 1:nruns){
x=rnorm(n,3,2)
y=b*x + rnorm(n, 0, 3)
fit=lm(y~x)
library("qgraph")
library("igraph")
library("IsingSampler")
library("IsingFit")
library('rstan')
set.seed(1337)
Kappa <- as.matrix(get.adjacency(watts.strogatz.game(1,10,1,0)))
Kappa[1,5] <- Kappa[5,1] <- Kappa[3,9] <- Kappa[9,3] <-Kappa[7,1] <-Kappa[1,7] <- Kappa[5,9] <-Kappa[9,5] <-Kappa[7,3] <-Kappa[3,7] <-.5
n=10000
beta=.4
#no effect
p=c()
for(i in 1:n){
s=rnorm(30,0,1)
sm=s+rnorm(30,0,1) #with measurement error
x=rnorm(30,0,1)
p=c(p,summary(lm(sm~x))$coefficients[2,4])
@cdriveraus
cdriveraus / effect of measurement error3.R
Last active February 17, 2017 16:13
Effect of measurement error mk2
N=10000 #iterations
obs=20 #observations per iteration
beta=exp(rnorm(N,-2,1))
Bs=c() #bias in estimated Beta (standardised) for no error condition
Bsm=c() #bias in estimated Beta (standardised) for with error condition
betas2=c() #unstandardised true beta for no error iterations that came out significant
betasm2=c() #unstandardised true beta for with error iterations that came out significant
for(i in 1:N){
@cdriveraus
cdriveraus / cor_stdreg.R
Last active March 14, 2019 16:02
Correlation and standardised regression
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))
@cdriveraus
cdriveraus / gist:27fa9048f6f683f7548929d6dfe38975
Created March 14, 2019 18:56
Common and unique process model
##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
install.packages('ctsem')
library(ctsem)
#### example 1
scode <- "
parameters {
real y[2];
}
model {
y[1] ~ normal(0, .5);
@cdriveraus
cdriveraus / Function code
Created November 26, 2019 14:26
devtools::run_examples problem
#' 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
cdriveraus / inequality.R
Created February 23, 2020 16:44
Earnings inequality
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)
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