R/run.sem.sim.R

In MonteCarloSEM: Monte Carlo Data Simulation Package

#' @title This function runs a model for simulated data by using the lavaan package.
#' @description Generated data sets (Generated by sim.skewed() or sim.normal() functions) will be fitted to pre-specified model. The lavaan package is used to fit the model.
#' After running the model, fit indices and parameters estimated with their standard errors will be printed to a Comma Separated Values (CSV) file.
#' The name of the output file is "All_Results.csv". Each line in the file represents the results of a simulation.
#' The columns are self-explanatory but the second column (named Notes) needs a more detailed explanation.
#' This column shows if the model convergency. If the model is converged without any problem the value will be "CONVERGE"
#' If it is not converged the value will be "NONCONVERGE" and all the values in the line will be "NA"
#' If there are some kind of warnings (such as negative variance) during the model run the value will be "WARNING" and based on the warning type some of the values might be "NA".
#' To run the simulation, previously generated (either with the package functions or any other software) data sets and the list of the data sets (i.e., "Data_List.dat") should be located in the same folder in the working directory.
#'
#' @author Fatih Orcan
#' @importFrom lavaan cfa parameterEstimates fitmeasures standardizedSolution lavTech lavInspect
#' @importFrom utils read.table write.csv
#' @param model Lavaan model
#' @param PEmethod Parameter estimation method. The default is ML.
#' @param dataList List of the names of data sets generated earlier either with the package functions or any other software.
#' @param f.loc File location. It indicates where the simulated data sets and "dataList" are located.
#' @param missing as in the lavaan package (See lavOptions)
#' @export
#' @examples
#'
#' #   Data needed to be generated at the first step.
#' fc<-fcors.value(nf=3, cors=c(1,.5,.6,.5,1,.4,.6,.4,1))
#' fl<-loading.value(nf=3, fl.loads=c(.5,.5,.5,0,0,0,0,0,0,0,0,.6,.6,.6,0,0,0,0,0,0,0,0,.4,.4))
#' sim.normal(nd=10, ss=100, fcors=fc, loading<-fl,  f.loc=tempdir())
#'
#' #  Then simulation should be run at the second step.
#' lavaanM<-'
#' #CFA Model
#' f1	=~ NA*x1 + x2 + x3
#' f2	=~ NA*x4 + x5 + x6
#' f3 =~ NA*x7 + x8
#' #Factor Correlations
#' f1	~~ f2
#' f1	~~ f3
#' f2	~~ f3
#' #Factor variance
#' f1	~~ 1*f1
#' f2	~~ 1*f2
#' f3	~~ 1*f3
#' '
#' dl<-"Data_List.dat"  # should be located in the working directory.
#'
#' # Please note that this function uses data sets and the list files which were generated previously.
#' # If there are no such a data sets and the list file, it will print an error message
#' #  saying "cannot open the connection"
#'
#' fit.simulation(model=lavaanM, PEmethod="MLR", dataList=dl, f.loc=tempdir())
#'

fit.simulation<-function(model, PEmethod="ML", dataList="Data_List.dat", f.loc,
                         missing=NULL){

  data.names<-read.table(paste(f.loc, "/", dataList,sep=""), header = FALSE)

  fit.names<-c("rep#","Notes", "chisq", "df", "pvalue","chisq.scaled", "df.scaled", "pvalue.scaled", "chisq.scaling.factor",
               "baseline.chisq", "baseline.df", "baseline.pvalue", "cfi","tli","srmr", "rmsea", "rmsea.ci.lower",
               "rmsea.ci.upper", "rmsea.pvalue", "logl","aic", "bic")
  par.names1<-c("est","se","pvalue")
  par.names2<-c("std.est","std.se","pvalue")

  veri<-read.table(paste(f.loc,"/", data.names[1,],sep=""))

  colnames(veri)<-c("ID", paste("x",seq(1:(dim(veri)[2]-1)),sep=""))
  veri<-veri[,-1]

  sonuc<-cfa(model,veri, estimator= PEmethod )

  #P.Est<-parameterEstimates(sonuc)
  #Sp.Est<-standardizedSolution(sonuc)
  tum.sonuc<-matrix(NA,dim(data.names)[1],(length(fit.names)+(dim(parameterEstimates(sonuc))[1])*6))

  for (i in 1:dim(data.names)[1]){  #Simulation replication start

    P.Est<-parameterEstimates(sonuc)
    Sp.Est<-standardizedSolution(sonuc)

    veri<-read.table(paste(f.loc,"/",data.names[i,], sep = ""))
    colnames(veri)<-c("ID", paste("x",seq(1:(dim(veri)[2]-1)),sep=""))
    veri<-veri[,-1]

    sonuc<-cfa(model,veri, estimator =PEmethod)
    tum.sonuc[i,1]<-i

    if(lavTech(sonuc, "converged")==TRUE){
      tum.sonuc[i,2:length(fit.names)]<-round(fitmeasures(sonuc)[fit.names[-1]],3)

      for(k in 1:(dim(P.Est)[1])){
        tum.sonuc[i,length(fit.names)+6*k-5]<-round(P.Est[k,"est"],3)
        tum.sonuc[i,length(fit.names)+6*k-4]<-round(P.Est[k,"se"],3)
        tum.sonuc[i,length(fit.names)+6*k-3]<-round(P.Est[k,"pvalue"],3)
        tum.sonuc[i,length(fit.names)+6*k-2]<-round(Sp.Est[k,"est.std"],3)
        tum.sonuc[i,length(fit.names)+6*k-1]<-round(Sp.Est[k,"se"],3)
        tum.sonuc[i,length(fit.names)+6*k]<-round(Sp.Est[k,"pvalue"],3)
      }

      print(paste(round(100*i/dim(data.names)[1],2),"% has completed...", sep=""))
      if(lavInspect(sonuc, "post.check")==TRUE){tum.sonuc[i,2]<-c("CONVERGED")}
      if(lavInspect(sonuc, "post.check")==FALSE){tum.sonuc[i,2]<-c("WARNING")}

    }

    if(lavTech(sonuc, "converged")==FALSE){
      for(k in 1:(dim(P.Est)[1])){
        tum.sonuc[i,]<-NA
        tum.sonuc[i,1]<-i
      }

      print(paste(round(100*i/dim(data.names)[1],2),"% has completed...", sep=""))
      tum.sonuc[i,2]<-c("NOT_CONVERGED")
    }
  } # replication end

  print("All Done !!!")
  colnames(tum.sonuc)<-c(paste("x",seq(1:(length(fit.names)+(dim(P.Est)[1])*6)),sep = ""))

  colnames(tum.sonuc)[1:length(fit.names)]<-c(fit.names)
  for(k in 1:(dim(P.Est)[1])){
    eft<-paste(P.Est[k,c("lhs","op","rhs")],sep = "", collapse="")
    colnames(tum.sonuc)[length(fit.names)+6*k-5]<-eft
    colnames(tum.sonuc)[length(fit.names)+6*k-4]<-par.names1[2]
    colnames(tum.sonuc)[length(fit.names)+6*k-3]<-par.names1[3]
    colnames(tum.sonuc)[length(fit.names)+6*k-2]<-par.names2[1]
    colnames(tum.sonuc)[length(fit.names)+6*k-1]<-par.names2[2]
    colnames(tum.sonuc)[length(fit.names)+6*k]<-par.names2[3]
  }
  write.csv(tum.sonuc, file= paste(f.loc,"/All_Results.csv", sep = ""), row.names = FALSE)
}