library(EPSY905R) # to run this script, you will first have to install JAGS from http://mcmc-jags.sourceforge.net dataIQ = read.csv(file = "iqdata.csv") #AUTOMATING PACKAGES NEEDED FOR ANALYSES-------------------------------------------------------------------- haspackage = require("lavaan") if (haspackage==FALSE){ install.packages("lavaan") } library("lavaan") haspackage = require("R2jags") if (haspackage==FALSE){ install.packages("R2jags") } library("R2jags") haspackage = require("nlme") if (haspackage==FALSE){ install.packages("nlme") } library("nlme") # IQ Data Examples ----------------------------------------------------------------------------------------- # estimation with Least Squares: model1LS = lm(y~1, data = dataIQ) summary(model1LS) # estimation with REML: model1REML = gls(y~1, data = dataIQ, method= "REML") summary(model1REML) summary(model1REML)$sigma^2 # estimation with ML: model1ML = gls(y~1, data = dataIQ, method= "ML") summary(model1ML) summary(model1ML)$sigma^2 # estimation with Bayesian model01Bayes = function(){ # likelihood for (i in 1:n){ y[i] ~ dnorm(mu, tau) } #priors mu ~ dnorm(100, 1/2.13^2) tau ~ dunif(1/400, 1/200) sigma2 = 1/tau } data = list(y = dataIQ$y, n = nrow(dataIQ)) jags.param = c("mu", "tau", "sigma2") fit <- jags.parallel(data=data, parameters.to.save=jags.param, n.iter=50000, n.chains=2,n.thin=2,n.burnin=40000, model.file=model01Bayes) chain = as.mcmc(fit) plot(chain) fit model02Bayes = function(){ # likelihood for (i in 1:n){ y[i] ~ dnorm(mu, tau) } #priors mu ~ dunif(-10000, 10000) tau ~ dunif(1/5000, 100000) sigma2 = 1/tau } data = list(y = dataIQ$y, n = nrow(dataIQ)) jags.param = c("mu", "tau", "sigma2") fit2 <- jags.parallel(data=data, parameters.to.save=jags.param, n.iter=50000, n.chains=2,n.thin=2,n.burnin=40000, model.file=model02Bayes) chain2 = as.mcmc(fit2) plot(chain2) fit2 model03Bayes = function(){ # likelihood for (i in 1:n){ y[i] ~ dnorm(mu, tau) } #priors mu ~ dunif(-10000, 10000) tau ~ dgamma(.01, .01) sigma2 = 1/tau } data = list(y = dataIQ$y, n = nrow(dataIQ)) jags.param = c("mu", "tau", "sigma2") fit3 <- jags.parallel(data=data, parameters.to.save=jags.param, n.iter=50000, n.chains=2,n.thin=2,n.burnin=40000, model.file=model03Bayes) chain3 = as.mcmc(fit3) plot(chain3) fit3 model04Bayes = function(){ # likelihood for (i in 1:n){ y[i] ~ dnorm(mu, tau) } #priors mu ~ dunif(102, 103) tau ~ dunif(1/242, 1/238) sigma2 = 1/tau } data = list(y = dataIQ$y, n = nrow(dataIQ)) jags.param = c("mu", "tau", "sigma2") fit4 <- jags.parallel(data=data, parameters.to.save=jags.param, n.iter=50000, n.chains=2,n.thin=2,n.burnin=40000, model.file=model04Bayes) chain4 = as.mcmc(fit4) plot(chain4) fit4 # BLAVAAN EXAMPLE (RUN AT YOUR OWN RISK) ----------------------------------------------------------------------- haspackage = require("blavaan") if (haspackage==FALSE){ install.packages("blavaan") } library("blavaan") data01 = EPSY905R::dataMath model01.syntax = " #endogenous variable equations perf ~ hsl + msc + mse use ~ mse mse ~ hsl + cc + female msc ~ mse + cc + hsl cc ~ hsl hsl ~ female #endogenous variable intercepts perf ~ 1 use ~ 1 mse ~ 1 msc ~ 1 cc ~ 1 hsl ~ 1 #endogenous variable residual variances perf ~~ perf use ~~ use mse ~~ mse msc ~~ msc cc ~~ cc hsl ~~ hsl #endogenous variable residual covariances #none specfied in the original model so these have zeros: perf ~~ 0*use + 0*mse + 0*msc + 0*cc + 0*hsl use ~~ 0*mse + 0*msc + 0*cc + 0*hsl mse ~~ 0*msc + 0*cc + 0*hsl msc ~~ 0*cc + 0*hsl cc ~~ 0*hsl " #estimate model model01.fit = sem(model01.syntax, data=data01, mimic = "MPLUS", estimator = "MLR") model01.fit_mcmc = bsem(model01.syntax, data=data01) plot(model01.fit_mcmc) summary(model01.fit_mcmc) summary(model01.fit) #see if model converged inspect(model01.fit, what="converged") #show summary of model fit statistics and parameters summary(model01.fit, standardized=TRUE, fit.measures=TRUE) #show normalized residual covariance matrix residuals(model01.fit ,type="normalized") #calculate modification indices model01.mi=inspect(model01.fit, what="mi") #order from highest to lowest model01.mi=model01.mi[with(model01.mi,order(-mi)),] #display values of indices model01.mi