Nothing
# A function that takes the output of Mplus and reads it to an mlVAR object. Can be used to make changes to codes.
importMplus <- function(outfile){
# Read output:
foo <- capture.output(output <- suppressMessages(suppressWarnings(readModels(outfile,what = "all"))))
# Read results into R:
# saveInfo <- suppressWarnings(getSavedata_Fileinfo(outfile))
saveInfo <- output$savedata_info
# Read samples into R:
bayesSamples <- output$bparameters$valid_draw
#suppressWarnings(getSavedata_Bparams(outfile) )
if (!is(bayesSamples,"mcmc")){
bayesSamples <- do.call(rbind,bayesSamples)
}
# Read data into R:
Data <- read.table(gsub("\\.out",".dat",outfile),na.strings = ".")
names(Data) <- unlist(strsplit(output$input$variable$names,split = " "))
# Read mlVAR model into R:
mlvarFile <- gsub("\\.out",".mlVARmodel",outfile)
if (!file.exists(mlvarFile)){
stop(paste0("mlVAR output not found, expected a file ",mlvarFile))
}
mlVARinfo <- read.table(mlvarFile,header=TRUE,stringsAsFactors=FALSE)
# Names of headers in samples:
sampleNames <- saveInfo$bayesVarNames
sampleNames <- gsub("^Parameter\\.\\d*\\_","",sampleNames)
### Collect the results ####
Results <- list()
# Outcomes:
Outcomes <- unique(mlVARinfo$dep)
nOut <- length(Outcomes)
# Predictors:
Predictors <- unique(mlVARinfo$pred)
nPred <- length(Predictors)
# Number of lags:
Lags <- sort(unique(mlVARinfo$lag))
Lags <- Lags[Lags>0]
nLag <- length(Lags)
# Number of samples:
nSample <- nrow(bayesSamples)
if (!identical(Outcomes,Predictors)){
stop("Mplus estimation with different predictors than outcomes is not yet supported.")
}
### BETWEEN SUBJECTS ###
# Always random
# Gather the variance--covariance matrices:
Omega_mu_cov_samples <- array(,c(nOut, nOut, nSample))
for (i in 1:nOut){
for (j in i:nOut){
# Name of column:
if (i == j){
colName <- paste0("%BETWEEN%:.",toupper(Outcomes[i]))
colID <- which(sampleNames == colName)
Omega_mu_cov_samples[i,j,] <- bayesSamples[,colID]
} else {
colName <- paste0("%BETWEEN%:.",toupper(Outcomes[j]),".WITH.",toupper(Outcomes[i]))
colID <- which(sampleNames == colName)
Omega_mu_cov_samples[i,j,] <- Omega_mu_cov_samples[j,i,] <- bayesSamples[,colID]
}
}
}
# Store to output:
Results[["Omega_mu"]] <- covSamples_Mplus(Omega_mu_cov_samples)
### TEMPORAL NETWORK ###
tempModel <- unique(mlVARinfo[mlVARinfo$model == "temporal","random"])
if (length(tempModel) > 1){
stop("Temporal model is a mixture of fixed/orthogonal/correlated")
}
subModel <- mlVARinfo %>% dplyr::filter(.data[['model']] == "temporal")
# Gather the parameter samples:
BetaFixed <- array(,c(nOut, nPred, nLag,nSample))
BetaSD <- array(0,c(nOut, nPred, nLag,nSample))
if (tempModel != "fixed"){
for (i in 1:nOut){
for (j in 1:nPred){
for (l in 1:nLag){
# Name of parameter:
parName <- mlVARinfo$par[mlVARinfo$dep == Outcomes[i] & mlVARinfo$pred == Predictors[j] & mlVARinfo$lag == Lags[l]]
if (length(parName) == 1){
# Fixed effect:
sampleParName <- paste0("%BETWEEN%:.MEAN.",toupper(parName))
sampleParID <- which(sampleNames == sampleParName)
BetaFixed[i,j,l,] <- bayesSamples[,sampleParID]
# Standard deviation:
sampleParName <- paste0("%BETWEEN%:.",toupper(parName))
sampleParID <- which(sampleNames == sampleParName)
BetaSD[i,j,l,] <- bayesSamples[,sampleParID]
}
}
}
}
} else {
### FIXED EFFECTS ###
for (i in 1:nOut){
for (j in 1:nPred){
for (l in 1:nLag){
# Name of column:
colName <- paste0("%WITHIN%:.",toupper(Outcomes[1]),".ON.",toupper(Predictors[j]),"&",Lags[l])
# Fixed effect:
sampleParID <- which(sampleNames == colName)
BetaFixed[i,j,l,] <- bayesSamples[,sampleParID]
}
}
}
}
# Store to output:
Results[["Beta"]] <- parSamples_Mplus(BetaFixed,BetaSD)
### CONTEMPORANEOUS EFFECTS ###
conModel <- unique(mlVARinfo[grepl("cont",mlVARinfo$model),"random"])
# Factor names:
submodel <- mlVARinfo[mlVARinfo$model == "cont cov",]
# For every cov, create a dummy factors with equal factor loadings:
dummyFactors <- paste0("Fac",seq_len(nrow(submodel)))
if (length(conModel) > 1){
stop("Contemporaneous model is a mixture of fixed/orthogonal/correlated")
}
if (conModel != "fixed"){
# First collect the matrix as is:
Theta_samples <- array(,c(nOut, nOut, nSample))
for (i in 1:nOut){
for (j in i:nOut){
# Name of column:
# Name of parameter:
parName <- mlVARinfo$par[mlVARinfo$dep == Outcomes[i] & mlVARinfo$pred == Predictors[j] & mlVARinfo$lag == 0]
if (length(parName) == 1){
# Fixed effect:
sampleParName <- paste0("%BETWEEN%:.MEAN.",toupper(parName))
sampleParID <- which(sampleNames == sampleParName)
Theta_samples[i,j,] <- Theta_samples[j,i,] <- exp(bayesSamples[,sampleParID])
}
}
}
# Read signs into R:
signsFile <- gsub("\\.out",".mlVARsigns",outfile)
if (!file.exists(signsFile)){
stop(paste0("mlVAR signs file not found, expected a file ",mlvarFile))
}
signs <- as.matrix(read.table(signsFile,header=FALSE,stringsAsFactors=FALSE))
# Now add sum of off-diagonals to each diagonal!
for (t in 1:dim(Theta_samples)[3]){
diag(Theta_samples[,,t]) <- colSums(Theta_samples[,,t])
# Multiply with signs:
Theta_samples[,,t] <- Theta_samples[,,t] * signs
}
} else {
# First collect the matrix as is:
Theta_samples <- array(,c(nOut, nOut, nSample))
for (i in 1:nOut){
for (j in i:nOut){
# Name of column:
if (i == j){
colName <- paste0("%WITHIN%:.",toupper(Outcomes[i]))
colID <- which(sampleNames == colName)
Theta_samples[i,j,] <- bayesSamples[,colID]
} else {
# Which dummy?
colName <- paste0("%WITHIN%:.",toupper(Outcomes[j]),".WITH.",toupper(Outcomes[i]))
colID <- which(sampleNames == colName)
Theta_samples[i,j,] <- Theta_samples[j,i,] <- bayesSamples[,colID]
}
}
}
}
# # Now add sum of off-diagonals to each diagonal!
# for (t in 1:dim(Theta_samples)[3]){
# diag(Theta_samples[,,t]) <- colSums(Theta_samples[,,t])
#
# # Multiply with signs:
# Theta_samples[,,t] <- Theta_samples[,,t] * signs
# }
# Store to output:
Results[["Theta"]] <- covSamples_Mplus(Theta_samples)
# Obtain DIC:
if (is.null(output$summaries$DIC)){
warning("No DIC obtained. Model might not have converged.")
DIC <- NA
pD <- NA
} else {
DIC <- output$summaries$DIC
pD <- output$summaries$pD
}
# Combine results:
Results <- list(
results = Results,
output = output,
fit = NULL,
data = Data,
model = mlVARinfo,
DIC = DIC,
pD = pD,
input = list(
vars = Outcomes,
lags = Lags,
temporal = tempModel,
contemporaneous = conModel,
estimator = sprintf("Mplus (version %s)",output$summaries$Mplus.version)
)
)
class(Results) <- "mlVAR"
return(Results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.