To conduct model checking and comparison analyses on the CFA model with one latent variable (factor):
Run the model in “One Factor CFA.bug” for with 3 chains for 5500 iterations, monitoring the parameters
Save the output via the “coda” function in BUGS. Place the files “coda1.txt”, “coda2.txt”, “coda3.txt”, and “codaIndex.txt” in the “coda files” folder. These are used in the model checking component.
Run the model in “One Factor CFA with CPO.bug” for with 3 chains for 5500 iterations, monitoring “inv.p.x”
Save the output via the “coda” function in BUGS. Rename the files “codaforCPO1.txt”, “codaforCPO2.txt”, “codaforCPO3.txt”, and “codaforCPOIndex.txt” and place them in the “coda files” folder. These are used in the model comparison component.
Run “Governing Code for CFA One Latent Variable Model Checking and Comparison.R”, which will call other R code.
# Read in the raw text file as a single string - this is in the WinBUGS format
data_text <- readLines("CFA-One-Latent-Variable/WinBUGS/One Factor CFA data.txt")
# Evaluate the text as R code to reconstruct the list
# This assigns values to n and J and all the x values in a way that R can deal with.
data_list <- eval(parse(text = data_text))
# Define initial values for each chain
# These have been copied from ints1.txt, inits2.txt, and inits3.txt and pasted here.)
inits_list <- list(
list(
tau = c(0.1, 0.1, 0.1, 0.1, 0.1),
lambda = c(NA, 0, 0, 0, 0),
inv.phi = 1,
inv.psi = c(1, 1, 1, 1, 1)
),
list(
tau = c(3, 3, 3, 3, 3),
lambda = c(NA, 3, 3, 3, 3),
inv.phi = 2,
inv.psi = c(2, 2, 2, 2, 2)
),
list(
tau = c(5, 5, 5, 5, 5),
lambda = c(NA, 6, 6, 6, 6),
inv.phi = 0.5,
inv.psi = c(0.5, 0.5, 0.5, 0.5, 0.5)
)
)
############################################################################################
# Define the
# main analysis folder
# R folder
# Data folder
# Coda folder
############################################################################################
# Using file.path() throughout the sourced script now, file.path() automatically inserts the correct separator between
main.folder <- "CFA-One-Latent-Variable"
R.folder <- file.path(main.folder, "R")
rjags.folder <- file.path(main.folder, "Rjags")
data.folder <- file.path(main.folder, "Data")
coda.folder <- file.path(main.folder, "coda_files")Initialize the model
model1 <- jags.model(file = file.path(rjags.folder, "One Factor CFA.jags"), data = data_list, inits = inits_list, n.chains = 3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2500
## Unobserved stochastic nodes: 515
## Total graph size: 8029
##
## Initializing model# Draw samples from iteration 1 to 5500
samples1 <- coda.samples(model1, variable.names = c("lambda", "tau", "phi", "ksi", "psi"), n.iter = 5500)Where instructed to use the coda function to create coda files in coda folder, instead build them in a way RJAGS can work with.
# Write individual chain outputs
# This has been modified to work in R from the original WinBUGS coda function
# Rebuild chain files without formatting issues
for (i in seq_along(samples1)) {
write.table(samples1[[i]],
file = paste0("CFA-One-Latent-Variable/coda_files/coda", i, ".txt"),
sep = " ", # force space
eol = "\n", # force Unix line endings
row.names = FALSE,
col.names = FALSE,
quote = FALSE,
fileEncoding = "UTF-8") # clean encoding
}
# Variable names include the 500 individual ksi[n], the five factors' lambda[J] (loading), psi[J] (variance of factor 1), and tau[J] (intercept), and phi.
var_names <- varnames(samples1)
index_lines <- sapply(seq_along(var_names), function(i) {
paste(var_names[i], i, i)
})
writeLines(index_lines, file.path(coda.folder, "codaIndex.txt"), useBytes = TRUE)Initialize the model
model2 <- jags.model(file = "CFA-One-Latent-Variable/Rjags/One Factor CFA with CPO.jags", data = data_list, inits = inits_list, n.chains = 3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2500
## Unobserved stochastic nodes: 515
## Total graph size: 10529
##
## Initializing modelsamples2 <- coda.samples(model2, variable.names = c("inv.psi"), n.iter = 5500)# Write individual chain outputs
for (i in 1:length(samples2)) {
write.table(samples2[[i]],
file = paste0("CFA-One-Latent-Variable/coda_files/codaforCPO", i, ".txt"),
sep = " ", # force space
eol = "\n", # force Unix line endings
row.names = FALSE,
col.names = FALSE,
quote = FALSE,
fileEncoding = "UTF-8") # clean encoding
}
# Create codaforCPOIndex.txt
varnames <- varnames(samples2)
index_lines <- character(length(varnames))
for (i in seq_along(varnames)) {
index_lines[i] <- paste(varnames[i], i, i)
}
writeLines(index_lines, "CFA-One-Latent-Variable/coda_files/codaforCPOIndex.txt")############################################################################################
# Read in the observed data
############################################################################################
file.name <- "IIS.dat"
x <- as.matrix(read.table(file.path(data.folder, file.name), header=TRUE))
##########################################################################################
# Extract features
# n is the number of examinees
# J is the number of observables
##########################################################################################
n <- nrow(x)
J <- ncol(x)############################################################################################
# Read in the draws from the coda folder
###########################################################################
# Read chains manually
c1 <- read.table(file.path(coda.folder, "coda1.txt"))
c2 <- read.table(file.path(coda.folder, "coda2.txt"))
c3 <- read.table(file.path(coda.folder, "coda3.txt"))
# Wrap them into mcmc objects
mcmc1 <- mcmc(c1, start = 1, thin = 1)
mcmc2 <- mcmc(c2, start = 1, thin = 1)
mcmc3 <- mcmc(c3, start = 1, thin = 1)
# Combine into mcmc.list
draws.from.bugs <- mcmc.list(mcmc1, mcmc2, mcmc3)
###########################################################################
# Declare burnin and thin if needed
###########################################################################
n.burnin=500
# Fix variable names issue by assigning colnames to match
lambda_names <- paste0("lambda[", 1:5, "]")
phi_name <- "phi"
tau_names <- paste0("tau[", 1:5, "]")
ksi_names <- paste0("ksi[", 1:500, "]")
psi_names <- paste0("psi[", 1:5, "]")
var_names <- c(lambda_names, phi_name, tau_names, ksi_names, psi_names)
for (i in 1:length(draws.from.bugs)) {
colnames(draws.from.bugs[[i]]) <- var_names
}
############################################################################################
# Select the iterations to summarize
############################################################################################
draws.to.analyze <- window(draws.from.bugs, start = n.burnin + 1)
############################################################################################
# Convert the draws to a matrix
############################################################################################
draws.to.analyze.as.matrix <- as.matrix(draws.to.analyze)Structure the draws performs reorganization and extraction of posterior samples from the MCMC output into structured arrays that mirror the parameters of a confirmatory factor analysis (CFA) model with one latent variable.
Each parameter is reshaped into arrays with meaningful dimensions:
| Parameter | Array Shape | Purpose |
|---|---|---|
lambda.iterations |
[n.iterations, J, M] |
Factor loadings per item |
tau.iterations |
[n.iterations, J] |
Intercepts per item |
psi.iterations |
[n.iterations, J, J] |
Diagonal matrix for residual variances |
kappa.iterations |
[n.iterations, M] |
Latent variable means (typically fixed) |
phi.iterations |
[n.iterations, M, M] |
Latent variable covariances |
ksi.iterations |
[n.iterations, n, M] |
Individual latent scores |
############################################################################################
# Structure the draws.
# This code has been brought in from the 'Structure the Draws.R' source code
############################################################################################
##########################################################################################
# Define the number of latent variables
##########################################################################################
M <- 1
##########################################################################################
# Define the total number of iterations
##########################################################################################
n.iters.total <- nrow(draws.to.analyze.as.matrix)
##########################################################################################
# Define empty arrays to be filled in
##########################################################################################
lambda.iterations <- array(NA, c(n.iters.total, J, M))
tau.iterations <- array(NA, c(n.iters.total, J))
psi.iterations <- array(NA, c(n.iters.total, J, J))
kappa.iterations <- array(NA, c(n.iters.total, M))
phi.iterations <- array(NA, c(n.iters.total, M, M))
ksi.iterations <- array(NA, c(n.iters.total, n, M))
##########################################################################################
# Code for a one-factor model
##########################################################################################
##########################################################################################
# Define
# 'lambda.fixed.and.free'
# 'tau.fixed.and.free'
# 'psi.fixed.and.free'
# 'kappa.fixed.and.free'
# 'phi.fixed.and.free'
# 'ksi.fixed.and.free'
# with
# numeric values for fixed values
# NA's for unknown parameters
##########################################################################################
lambda.fixed.and.free <- matrix(NA, nrow=J, ncol=M) # Empty matrix to hold the lambda values
lambda.fixed.and.free[1,1] = 1 # fixes the first loading factor to 1
tau.fixed.and.free <- rep(NA, J) # Empty vector of length J to hold all freely estimated loading intercepts
psi.fixed.and.free <- diag(as.numeric(NA), J, J) # Empty diagonal matrix with NA on the diag all others fixed at 0
kappa.fixed.and.free <- rep(0, M) # Fix this vector of latent var means to 0 for as many latent variables as are in the model (1)
phi.fixed.and.free <- matrix(NA, M, M) # This matrix of latent factor variance is freely estimated
ksi.fixed.and.free <- matrix(NA, n, M) # This matrix of individual latent scores is freely estimated
##########################################################################################
# Loop over iterations
##########################################################################################
for(which.iter in 1:n.iters.total){
##########################################################################################
# Define the parameters for this matrix based on the fixed and free versions
##########################################################################################
lambda.iterations[which.iter, , ] <- lambda.fixed.and.free
tau.iterations[which.iter, ] <- tau.fixed.and.free
psi.iterations[which.iter, , ] <- psi.fixed.and.free
kappa.iterations[which.iter, ] <- kappa.fixed.and.free
phi.iterations[which.iter, , ] <- phi.fixed.and.free
ksi.iterations[which.iter, , ] <- ksi.fixed.and.free
##########################################################################################
# Fill in the loadings
##########################################################################################
for(j in 1:J){
for(m in 1:M){
if(is.na(lambda.fixed.and.free[j,m])){
lambda.iterations[which.iter,j,m] = draws.to.analyze.as.matrix[which.iter, paste("lambda[", j, "]", sep="")]
}
} # closes loop over m
} # closes loop over j
##########################################################################################
# Fill in the intercepts
##########################################################################################
for(j in 1:J){
if(is.na(tau.fixed.and.free[j])){
tau.iterations[which.iter,j] = draws.to.analyze.as.matrix[which.iter, paste("tau[", j, "]", sep="")]
}
} # closes loop over j
##########################################################################################
# Fill in the error variances
##########################################################################################
for(j in 1:J){
#for(jj in 1:J){
if(is.na(psi.fixed.and.free[j,j])){
psi.iterations[which.iter,j,j] = draws.to.analyze.as.matrix[which.iter, paste("psi[", j, "]", sep="")]
}
#} # closes loop over jj
} # closes loop over j
##########################################################################################
# Fill in the latent variable means
##########################################################################################
for(j in 1:M){
if(is.na(kappa.fixed.and.free[m])){
kappa.iterations[which.iter,m] = draws.to.analyze.as.matrix[which.iter, paste("kappa", sep="")]
}
} # closes loop over j
##########################################################################################
# Fill in the covariance matrix of latent variables
##########################################################################################
for(m in 1:M){
for(mm in 1:M){
if(is.na(phi.fixed.and.free[m,mm])){
phi.iterations[which.iter,m,mm] = draws.to.analyze.as.matrix[which.iter, paste("phi", sep="")]
}
} # closes loop over mm
} # closes loop over m
##########################################################################################
# Fill in the latent variable values
##########################################################################################
for(i in 1:n){
for(m in 1:M){
if(is.na(ksi.fixed.and.free[i,m])){
ksi.iterations[which.iter,i,m] = draws.to.analyze.as.matrix[which.iter, paste("ksi[", i, "]", sep="")]
}
} # closes loop over m
} # closes loop over i
} # closes loop over iterations############################################################################################
# Conduct the PPMC analyses
############################################################################################
# Create a subfolder if one isn't there
PPMC.folder <- file.path(main.folder, "PPMC_Analyses")
PP.datasets.folder <- file.path(PPMC.folder, "PP_Datasets")
dir.create(PPMC.folder, showWarnings = FALSE, recursive = TRUE)
dir.create(PP.datasets.folder, recursive = TRUE, showWarnings = FALSE)
# Set global output folder for plots and tables
output.folder <- PPMC.folder
write.axis.labels=TRUE
file.name <- "PPMC Analyses.R"
source(file.path(R.folder, file.name))## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs produced
## Warning in log(det(mod.imp.cov.matrix)): NaNs producedThis code doesn’t compute a Bayes Factor unless there is a second model’s LPML available. It just saves the LPML, which could be compared later if you run a second model (e.g., 2-factor CFA) and compute its LPML in the same way.
############################################################################################
# Create a subfolder if one isn't there
############################################################################################
Bayes.factor.folder <- file.path(main.folder, "Bayes Factor")
dir.create(Bayes.factor.folder, showWarnings = FALSE, recursive = TRUE)###########################################################################
# Read in the draws
# Using the read.coda function as the index name is different
###########################################################################
# Read cleaned chain manually
chain1_data <- read.table(file.path(coda.folder, "codaforCPO1.txt"), header = FALSE)
chain2_data <- read.table(file.path(coda.folder, "codaforCPO2.txt"), header = FALSE)
chain3_data <- read.table(file.path(coda.folder, "codaforCPO3.txt"), header = FALSE)
# Extract variable names from index file
var_names <- readLines(file.path(coda.folder, "codaforCPOIndex.txt"))
var_names <- sub(" .*", "", var_names) # strip off indices
colnames(chain1_data) <- var_names
colnames(chain2_data) <- var_names
colnames(chain3_data) <- var_names
# Convert to MCMC objects
chain1 <- mcmc(chain1_data)
chain2 <- mcmc(chain2_data)
chain3 <- mcmc(chain3_data)
# Combine into mcmc.list
draws.cpo <- mcmc.list(chain1, chain2, chain3)
# Convert the draws to matrix format
draws.cpo.as.matrix <- as.matrix(draws.cpo)
# Declare burn-in and thin if needed
n.burnin <- 500
# Trim the draws to exclude burn-in
draws.to.analyze <- window(draws.cpo, start = n.burnin + 1)
draws.to.analyze.as.matrix <- as.matrix(draws.to.analyze)
dim(draws.to.analyze.as.matrix)## [1] 15000 5Code from Summarize Posterior.R is pasted below no need to source
#############################################################################################
# Combine chains for summaries
#############################################################################################
coda.options(combine.stats=TRUE, combine.plots=TRUE)
#############################################################################################
# Extract the summary statistics
#############################################################################################
summary.stats <- summary(draws.to.analyze)
summary.stats$statistics## Mean SD Naive SE Time-series SE
## inv.psi[1] 1.719715 0.1334875 0.001089921 0.001602461
## inv.psi[2] 1.915477 0.1402696 0.001145296 0.001636833
## inv.psi[3] 2.082781 0.1657126 0.001353038 0.002175898
## inv.psi[4] 1.979501 0.1477880 0.001206684 0.001744413
## inv.psi[5] 1.852429 0.1472060 0.001201932 0.001886410summary.stats$quantiles## 2.5% 25% 50% 75% 97.5%
## inv.psi[1] 1.469772 1.627095 1.714782 1.806763 1.996843
## inv.psi[2] 1.656461 1.818285 1.909655 2.005496 2.206027
## inv.psi[3] 1.777551 1.967165 2.075257 2.189518 2.429536
## inv.psi[4] 1.708979 1.877965 1.973412 2.074387 2.285917
## inv.psi[5] 1.579080 1.750080 1.847248 1.946388 2.160508probability.for.HPD=.95
summary.statistics <- cbind(
summary.stats$statistics,
summary.stats$quantiles,
matrix(HPDinterval(draws.to.analyze, prob=probability.for.HPD)[[1]], ncol=2)
)
colnames(summary.statistics) <- c(
colnames(summary.stats$statistics),
colnames(summary.stats$quantiles),
c("95% HPD lower", "95% HPD Upper")
)This has been broguht here from source file ’Compute CPO and Pseudomarginal likelihood
#############################################################################################
# Compute CPO as inverse of the posterior mean (of the inverse likelihoods) from MCMC
#############################################################################################
CPO <- 1/summary.statistics[,1]
#############################################################################################
# Compute the log of the Psueduomarginal likelihood as the sum of the logged CPOs
#############################################################################################
log.PsML <- sum(log(CPO))
#############################################################################################
# Write out the CPOs
#############################################################################################
write.table(
x = matrix(CPO, ncol = 1),
file = file.path(Bayes.factor.folder, "CPO.out"),
row.names = FALSE,
col.names = FALSE
)
#############################################################################################
# Write out the log PsML
#############################################################################################
write(
x = log.PsML,
file = file.path(Bayes.factor.folder, "Log of PSML.out")
)# Reorient around the 2-factor model folder
main.folder <- "CFA-Two-Latent-Variables"
R.folder <- file.path(main.folder, "R")
data.folder <- file.path(main.folder, "Data")
coda.folder <- file.path(main.folder, "coda_files")
winbugs.folder <- file.path(main.folder, "WinBUGS")
rjags.folder <- file.path(main.folder, "Rjags")
# Read in the raw text file as a single string - this is in the WinBUGS format
data_text2 <- readLines(file.path(winbugs.folder, "Two Factor CFA data.txt"))
# Evaluate the text as R code to reconstruct the list
# This assigns values to n and J and all the x values in a way that R can deal with.
data_list2 <- eval(parse(text = data_text2)) inits1 <- list(
tau = rep(0.1, 5),
lambda = {
mat <- matrix(NA, 5, 2)
mat[2,1] <- 2.0
mat[3,1] <- 2.0
mat[5,2] <- 2.0
mat
},
inv.phi = matrix(c(1.0, 0.0, 0.0, 1.0), nrow = 2),
inv.psi = rep(1.0, 5)
)
inits2 <- list(
tau = rep(3.0, 5),
lambda = {
mat <- matrix(NA, 5, 2)
mat[2,1] <- 0.5
mat[3,1] <- 0.5
mat[5,2] <- 0.5
mat
},
inv.phi = matrix(c(1.33, -0.667, -0.667, 1.33), nrow = 2),
inv.psi = rep(2.0, 5)
)
inits3 <- list(
tau = rep(5.0, 5),
lambda = {
mat <- matrix(NA, 5, 2)
mat[2,1] <- 1.0
mat[3,1] <- 1.0
mat[5,2] <- 1.0
mat
},
inv.phi = matrix(c(1.96, -1.37, -1.37, 1.96), nrow = 2),
inv.psi = rep(0.5, 5)
)
# Bundle inits
inits_list2 <- list(inits1, inits2, inits3)model3 <- jags.model(file = file.path(rjags.folder, "Two Factor CFA.jags"), data = data_list2, inits = inits_list2, n.chains = 3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2500
## Unobserved stochastic nodes: 514
## Total graph size: 9035
##
## Initializing model# Draw samples from iteration 501 to 5500
samples3 <- coda.samples(model3, variable.names = c("lambda", "tau", "phi", "ksi", "psi"), n.iter = 5500)This has been modified to work in R from the original WinBUGS coda function
# Rebuild chain files without formatting issues
for (i in seq_along(samples3)) {
write.table(samples3[[i]],
file = file.path(coda.folder, paste0("coda", i, ".txt")),
sep = " ", # force space
eol = "\n", # force Unix line endings
row.names = FALSE,
col.names = FALSE,
quote = FALSE,
fileEncoding = "UTF-8") # clean encoding
}
# Variable names include the 500 individual ksi[n], the five factors' lambda[J] (loading), psi[J] (variance of factor 1), and tau[J] (intercept), and phi.
var_names <- varnames(samples3)
index_lines <- sapply(seq_along(var_names), function(i) {
paste(var_names[i], i, i)
})
writeLines(index_lines, file.path(coda.folder, "codaIndex.txt"), useBytes = TRUE)Initialize the model
model4 <- jags.model(file = file.path(rjags.folder, "Two Factor CFA with CPO.jags"), data = data_list2, inits = inits_list2, n.chains = 3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2500
## Unobserved stochastic nodes: 514
## Total graph size: 15561
##
## Initializing modelsamples4 <- coda.samples(model4, variable.names = c("inv.psi"), n.iter = 5500)# Write individual chain outputs
for (i in 1:length(samples4)) {
write.table(samples4[[i]],
file = file.path(coda.folder, paste0("codaforCPO", i, ".txt")),
sep = " ", # force space
eol = "\n", # force Unix line endings
row.names = FALSE,
col.names = FALSE,
quote = FALSE,
fileEncoding = "UTF-8") # clean encoding
}
# Create codaforCPOIndex.txt
varnames <- varnames(samples4)
index_lines <- character(length(varnames))
for (i in seq_along(varnames)) {
index_lines[i] <- paste(varnames[i], i, i)
}
writeLines(index_lines, file.path(coda.folder, "codaforCPOIndex.txt"))########################################################################################################################
# Read in the observed data
########################################################################################################################
file.name <- "IIS.dat"
x <- as.matrix(read.table(file.path(data.folder, file.name), header=TRUE))
############################################################################################################################
# Extract features
# n is the number of examinees
# J is the number of observables
############################################################################################################################
n <- nrow(x)
J <- ncol(x)
########################################################################################################################
# Conduct Bayes Factor analyses
# Create a subfolder if one isn't there
########################################################################################################################
Bayes.factor.folder <- file.path(main.folder, "Bayes Factor")
dir.create(Bayes.factor.folder, showWarnings = FALSE, recursive = TRUE)
# Read cleaned chain manually
c.1 <- read.table(file.path(coda.folder, "codaforCPO1.txt"), header = FALSE)
c.2 <- read.table(file.path(coda.folder, "codaforCPO2.txt"), header = FALSE)
c.3 <- read.table(file.path(coda.folder, "codaforCPO3.txt"), header = FALSE)
# Extract variable names from index file
var_names <- readLines(file.path(coda.folder, "codaforCPOIndex.txt"))
var_names <- sub(" .*", "", var_names) # strip off indices
colnames(c.1) <- var_names
colnames(c.2) <- var_names
colnames(c.3) <- var_names
# Convert to MCMC objects
chain1 <- mcmc(c.1)
chain2 <- mcmc(c.2)
chain3 <- mcmc(c.3)
# Combine into mcmc.list
draws.cpo <- mcmc.list(chain1, chain2, chain3)
# Convert the draws to matrix format
draws.cpo.as.matrix <- as.matrix(draws.cpo)
# Declare burn-in and thin if needed
n.burnin <- 500
# Trim the draws to exclude burn-in
draws.to.analyze <- window(draws.cpo, start = n.burnin + 1)
draws.to.analyze.as.matrix <- as.matrix(draws.to.analyze)
dim(draws.to.analyze.as.matrix)## [1] 15000 5Code from Summarize Posterior.R is pasted below no need to source
#############################################################################################
# Combine chains for summaries
#############################################################################################
coda.options(combine.stats=TRUE, combine.plots=TRUE)
#############################################################################################
# Extract the summary statistics
#############################################################################################
summary.stats <- summary(draws.to.analyze)
summary.stats$statistics## Mean SD Naive SE Time-series SE
## inv.psi[1] 1.697288 0.1335982 0.001090825 0.002257104
## inv.psi[2] 1.961008 0.1503670 0.001227741 0.001940108
## inv.psi[3] 2.147196 0.1813686 0.001480868 0.002700056
## inv.psi[4] 2.120295 0.1812764 0.001480116 0.003451615
## inv.psi[5] 1.944745 0.1777444 0.001451277 0.003072972summary.stats$quantiles## 2.5% 25% 50% 75% 97.5%
## inv.psi[1] 1.453832 1.604943 1.692590 1.783870 1.973574
## inv.psi[2] 1.681310 1.857467 1.954695 2.057121 2.274558
## inv.psi[3] 1.818913 2.020899 2.136870 2.263899 2.527855
## inv.psi[4] 1.791924 1.994716 2.109366 2.237516 2.500667
## inv.psi[5] 1.628816 1.821612 1.933874 2.058081 2.328866probability.for.HPD=.95
summary.statistics <- cbind(
summary.stats$statistics,
summary.stats$quantiles,
matrix(HPDinterval(draws.to.analyze, prob=probability.for.HPD)[[1]], ncol=2)
)
colnames(summary.statistics) <- c(
colnames(summary.stats$statistics),
colnames(summary.stats$quantiles),
c("95% HPD lower", "95% HPD Upper")
)This has been broguht here from source file ’Compute CPO and Pseudomarginal likelihood
#############################################################################################
# Compute CPO as inverse of the posterior mean (of the inverse likelihoods) from MCMC
#############################################################################################
CPO <- 1/summary.statistics[,1]
#############################################################################################
# Compute the log of the Psueduomarginal likelihood as the sum of the logged CPOs
#############################################################################################
log.PsML <- sum(log(CPO))
#############################################################################################
# Write out the CPOs
#############################################################################################
write.table(
x = matrix(CPO, ncol = 1),
file = file.path(Bayes.factor.folder, "CPO.out"),
row.names = FALSE,
col.names = FALSE
)
#############################################################################################
# Write out the log PsML
#############################################################################################
write(
x = log.PsML,
file = file.path(Bayes.factor.folder, "Log of PSML.out")
)View summary of CFA One Latent Var model so far
# Just the model parameters, not the 500 subject-level summary
# Step 1: Define the parameter names you want
model_pars <- c(
paste0("lambda[", 1:5, "]"),
"phi",
paste0("psi[", 1:5, "]"),
paste0("tau[", 1:5, "]")
)
# Step 2: Subset the samples before summarizing model parameters
samples_model1 <- samples1[, model_pars]
# Step 3: Summarize the subset
summary_model1 <- summary(samples_model1)
# Step 4: Combine stats and quantiles into one clean table
summary_df1 <- cbind(
summary_model1$statistics,
summary_model1$quantiles
)
# Step 5 (optional): Round and print
print(round(summary_df1, 3))## Mean SD Naive SE Time-series SE 2.5% 25% 50% 75% 97.5%
## lambda[1] 1.000 0.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000
## lambda[2] 0.788 0.110 0.001 0.003 0.646 0.735 0.784 0.833 0.936
## lambda[3] 0.903 0.114 0.001 0.003 0.753 0.848 0.898 0.951 1.058
## lambda[4] 0.827 0.108 0.001 0.003 0.686 0.773 0.822 0.873 0.973
## lambda[5] 0.928 0.119 0.001 0.003 0.774 0.871 0.924 0.979 1.088
## phi 0.445 0.071 0.001 0.002 0.353 0.410 0.442 0.477 0.548
## psi[1] 0.587 0.046 0.000 0.001 0.502 0.555 0.584 0.616 0.682
## psi[2] 0.524 0.039 0.000 0.000 0.453 0.497 0.523 0.549 0.606
## psi[3] 0.484 0.039 0.000 0.000 0.412 0.457 0.482 0.509 0.561
## psi[4] 0.508 0.038 0.000 0.000 0.437 0.482 0.507 0.533 0.586
## psi[5] 0.543 0.043 0.000 0.001 0.465 0.514 0.541 0.571 0.632
## tau[1] 3.766 0.053 0.000 0.001 3.675 3.736 3.766 3.796 3.856
## tau[2] 3.722 0.052 0.000 0.001 3.642 3.694 3.721 3.749 3.800
## tau[3] 3.530 0.053 0.000 0.001 3.448 3.502 3.529 3.557 3.611
## tau[4] 3.845 0.052 0.000 0.001 3.764 3.817 3.843 3.871 3.924
## tau[5] 3.704 0.056 0.000 0.001 3.619 3.674 3.704 3.733 3.787lambda are the one factor model’s loadings on each of five constructs tau are observable variable intercepts psi are covariance of delta errors phi is the covariance \(\Phi\) of observables? means are represented as \(\kappa\) in the model
\(x_{i} = \tau + {\Lambda \xi_{i}} + \delta_{i'}\) # Now for the withCPO model monitoring inv.psi
summary(samples2)##
## Iterations = 1:5500
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 5500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## inv.psi[1] 1.718 0.1339 0.001042 0.001624
## inv.psi[2] 1.915 0.1403 0.001093 0.001567
## inv.psi[3] 2.083 0.1665 0.001296 0.002143
## inv.psi[4] 1.979 0.1479 0.001152 0.001659
## inv.psi[5] 1.852 0.1475 0.001148 0.001793
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## inv.psi[1] 1.467 1.626 1.714 1.805 1.994
## inv.psi[2] 1.656 1.818 1.909 2.005 2.205
## inv.psi[3] 1.776 1.968 2.076 2.190 2.430
## inv.psi[4] 1.708 1.878 1.973 2.075 2.285
## inv.psi[5] 1.579 1.750 1.847 1.946 2.159Reconsider the CFA model with one latent variable for the data from the Institutional Integration Scale described in Section 9.3.a. Recreate the Bayesian analysis reported there using WinBUGS and:
dic_result2 <- dic.samples(model2, n.iter = 5000, type = "pD") # or type = "popt"
print(dic_result2)## Mean deviance: 5349
## penalty 402.6
## Penalized deviance: 5752LR.fit.ppp.value.out is empty with only an NA, so describe here:
# Realized LR is from :
summary(realized.LR.fit)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 3190 3805 3929 4062 4019 1885432 384# Posterior Predictive LR:
summary(postpred.LR.fit)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 4 63 143 2689 352 22615563 384# Compute the posterior predictive p-value
mean(postpred.LR.fit > realized.LR.fit, na.rm = TRUE)## [1] 0.02647783This ppp p-value should be between .025 and .975 This is very close to the edge of the 95% credible interval at .025, suggesting that the posterior predictive LR is centered much lower than the realized LR.
# Describe the realized SRMR distribution
summary(realized.SRMR)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.251 18.592 19.554 19.140 20.101 27.378No connection between realized discrepancy measure and the actual data
# Describe the posterior predictive SRMR distribution
summary(postpred.SRMR)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.006976 0.023529 0.031056 0.033084 0.040559 0.110677Good fit of model to predicted values to the actual data. The model can produce simulated data that are a plausible fit to the observed data.
# Compute the posterior predictive p-value
mean(postpred.SRMR > realized.SRMR, na.rm = TRUE)## [1] 0Standardized Root Mean Square Residuals are incredibly high for realized SRMR, showing a drastic misfit to the observed data based on fixed parameters.
Question: Why does realized discrepancy matter? They are more like a frequentist SRMR, in that they are generated from a single point estimate theta selected from the distribution of thetas. The Bayesian posterior predictive SRMR smooths out by integrating over the posterior.
With SRMR of the correlations between five observed variables drastically different from the SRMR, that means that the model fails to explain the patterns in the observed data.
Interpret these results in terms of the adequacy of model-data fit. What feature of the model is responsible for this (in)adequacy of fit?
summary2 <- as.data.frame(data_list)[3:7]
obs_means <- colMeans(summary2)
pred_means <- colMeans(x.postpred)
colMeans(t(pred_means) > obs_means)## V1 V2 V3 V4 V5
## 0 0 0 0 0Not a single posterior predicted value is above the observed mean for each category.
# Convert the matrix to a dataframe with variable names
x_df <- as.data.frame(data_list$x)
colnames(x_df) <- paste0("x", 1:5) # or use meaningful names if available
# Define the one-factor model
model <- '
F1 =~ x1 + x2 + x3 + x4 + x5
'
# Fit the model using maximum likelihood
fit <- cfa(model, data = x_df, estimator = "ML")
# Summary with standardized estimates and fit measures
summary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 23 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 10
##
## Number of observations 500
##
## Model Test User Model:
##
## Test statistic 6.655
## Degrees of freedom 5
## P-value (Chi-square) 0.248
##
## Model Test Baseline Model:
##
## Test statistic 529.098
## Degrees of freedom 10
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.997
## Tucker-Lewis Index (TLI) 0.994
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3010.716
## Loglikelihood unrestricted model (H1) -3007.389
##
## Akaike (AIC) 6041.433
## Bayesian (BIC) 6083.579
## Sample-size adjusted Bayesian (SABIC) 6051.838
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.026
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.071
## P-value H_0: RMSEA <= 0.050 0.766
## P-value H_0: RMSEA >= 0.080 0.020
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.018
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## x1 1.000 0.563 0.596
## x2 0.906 0.094 9.653 0.000 0.510 0.592
## x3 1.050 0.101 10.385 0.000 0.591 0.670
## x4 0.948 0.096 9.902 0.000 0.534 0.616
## x5 1.068 0.105 10.210 0.000 0.601 0.649
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.574 0.044 13.054 0.000 0.574 0.644
## .x2 0.482 0.037 13.119 0.000 0.482 0.650
## .x3 0.430 0.037 11.719 0.000 0.430 0.551
## .x4 0.466 0.037 12.749 0.000 0.466 0.620
## .x5 0.497 0.041 12.159 0.000 0.497 0.579
## F1 0.317 0.050 6.374 0.000 1.000 1.000The SRMR from the ML estimation model is 0.018 which is below 0.08 and considered a good fit. This compared to the realized SRMR of the Bayesian CFA model, which is 19.5, showing misspecification of the model.
The CFI and TLI are both above .95 which are evidence that the one-factor model is an excellent fit to compared to the baseline independence model.
# Convert the matrix to a dataframe with variable names
x_df2 <- as.data.frame(data_list$x)
colnames(x_df2) <- paste0("x", 1:5) # or use meaningful names if available
# Define the one-factor model
model <- '
F1 =~ x2 + x3 + x4
F2 =~ x1 + x4
'
# Fit the model using maximum likelihood
fit2 <- cfa(model, data = x_df2, estimator = "ML")## Warning: lavaan->lav_model_vcov():
## Could not compute standard errors! The information matrix could not be
## inverted. This may be a symptom that the model is not identified.# Summary with standardized estimates and fit measures
summary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 23 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 10
##
## Number of observations 500
##
## Model Test User Model:
##
## Test statistic 6.655
## Degrees of freedom 5
## P-value (Chi-square) 0.248
##
## Model Test Baseline Model:
##
## Test statistic 529.098
## Degrees of freedom 10
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.997
## Tucker-Lewis Index (TLI) 0.994
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3010.716
## Loglikelihood unrestricted model (H1) -3007.389
##
## Akaike (AIC) 6041.433
## Bayesian (BIC) 6083.579
## Sample-size adjusted Bayesian (SABIC) 6051.838
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.026
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.071
## P-value H_0: RMSEA <= 0.050 0.766
## P-value H_0: RMSEA >= 0.080 0.020
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.018
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## x1 1.000 0.563 0.596
## x2 0.906 0.094 9.653 0.000 0.510 0.592
## x3 1.050 0.101 10.385 0.000 0.591 0.670
## x4 0.948 0.096 9.902 0.000 0.534 0.616
## x5 1.068 0.105 10.210 0.000 0.601 0.649
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.574 0.044 13.054 0.000 0.574 0.644
## .x2 0.482 0.037 13.119 0.000 0.482 0.650
## .x3 0.430 0.037 11.719 0.000 0.430 0.551
## .x4 0.466 0.037 12.749 0.000 0.466 0.620
## .x5 0.497 0.041 12.159 0.000 0.497 0.579
## F1 0.317 0.050 6.374 0.000 1.000 1.000