psindex: Parameter Stability Index

# export to Mplus syntax

lav2mplus <- function(lav, group.label=NULL) {
    lav <- lav2check(lav)
    header <- "  ! this model syntax is autogenerated by lavExport\n"
    footer <- "\n"

    lav <- as.data.frame(lav, stringsAsFactors=FALSE)
    ngroups <- lav_partable_ngroups(lav)

    lav_one_group <- function(lav) {

        # mplus does not like variable names with a 'dot'
        # replace them by an underscore '_'
        lav$lhs <- gsub("\\.", "_", lav$lhs)
        lav$rhs <- gsub("\\.", "_", lav$rhs)
 
        # remove contraints (:=, <, >, ==) here
        con.idx <- which(lav$op %in% c(":=", "<",">","=="))
        if(length(con.idx) > 0L) {
            lav <- lav[-con.idx,]
        }

        # remove exogenous variances/covariances/intercepts...
        exo.idx <- which(lav$exo == 1L)
        if(length(exo.idx)) {
            lav <- lav[-exo.idx,]
        }

        # remove intercepts for categorical variables
        ord.names <- unique(lav$lhs[ lav$op == "|" ])
        ord.int.idx <- which(lav$op == "~1" & lav$lhs %in% ord.names)
        if(length(ord.int.idx)) {
            lav <- lav[-ord.int.idx,]
        }

        # end of line
        lav$eol <- rep(";", length(lav$lhs))
        lav$ustart <- ifelse(is.na(lav$ustart), "", lav$ustart)
        lav$rhs2 <- ifelse(lav$free == 0L, 
                           paste("@",lav$ustart,sep=""),
                           paste("*",lav$ustart,sep=""))
        lav$plabel <- gsub("\\.", "", lav$plabel)
        lav$plabel <- ifelse(lav$plabel == "", lav$plabel,
                            paste(" (",lav$plabel,")",sep=""))

        # remove variances for ordered variables
        ov.names.ord <- vnames(lav, type="ov.ord")
        ord.idx <- which(lav$lhs %in% ov.names.ord &
                         lav$op == "~~" &
                         lav$free == 0L &
                         lav$lhs == lav$rhs)
        lav$lhs[ord.idx] <- paste("! ", lav$lhs[ord.idx], sep="")
        lav$op[ord.idx] <- ""
        lav$rhs[ord.idx] <- ""

        # variances
        var.idx <- which(lav$op == "~~" & lav$rhs == lav$lhs)
        lav$op[var.idx] <- ""
        lav$rhs[var.idx] <- ""

        # scaling factors
        scal.idx <- which(lav$op == "~*~") 
        lav$op[scal.idx] <- ""
        lav$rhs2[scal.idx] <- paste(lav$rhs2[scal.idx],"}",sep="")
        lav$lhs[scal.idx] <- "{"

        # intercepts - excluding categorical observed
        int.idx <- which(lav$op == "~1")
        lav$op[int.idx] <- ""
        lav$rhs2[int.idx] <- paste(lav$rhs2[int.idx],"]",sep="")
        lav$lhs[int.idx] <- paste("[", lav$lhs[int.idx],sep="")

        # thresholds
        th.idx <- which(lav$op == "|")
        lav$op[th.idx] <- "$"
        lav$rhs[th.idx] <- gsub("t", "", x=lav$rhs[th.idx])
        lav$rhs2[th.idx] <- paste(lav$rhs2[th.idx],"]",sep="")
        lav$lhs[th.idx] <- paste("[", lav$lhs[th.idx],sep="")

        # replace binary operators
        lav$op <- ifelse(lav$op == "=~", " BY ", lav$op)
        lav$op <- ifelse(lav$op == "~", " ON ", lav$op)
        lav$op <- ifelse(lav$op == "~~", " WITH ", lav$op)


        lav2 <- paste(lav$lhs, lav$op, lav$rhs, lav$rhs2,
                      lav$plabel, lav$eol, sep="")
                      
        body <- paste(" ", lav2, collapse="\n")

        body
    }

    if(ngroups == 1L) {
        body <- lav_one_group(lav)
    } else {
        # group 1
        body <- lav_one_group(lav[lav$group == 1,])

        if(is.null(group.label) || length(group.label) == 0L) {
            group.label <- paste(1:ngroups)
        }

        for(g in 2:ngroups) {
            body <- paste(body,
                          paste("\nMODEL ", group.label[g], ":\n", sep=""),
                          lav_one_group(lav[lav$group == g,]),
                          sep="")
        }
    }

    # constraints go to a 'MODEL CONSTRAINTS' block
    con.idx <- which(lav$op %in% c(":=", "<",">","=="))
    if(length(con.idx) > 0L) {
        ### FIXME: we need to convert the operator
        ### eg b^2 --> b**2, others??
        lav$lhs[con.idx] <- gsub("\\^","**",lav$lhs[con.idx])
        lav$rhs[con.idx] <- gsub("\\^","**",lav$rhs[con.idx])

        constraints <- "\nMODEL CONSTRAINT:\n"
        # define 'new' variables
        def.idx <- which(lav$op == ":=")
        if(length(def.idx) > 0L) {
            def <- paste(lav$lhs[def.idx], collapse= " ")
            constraints <- paste(constraints, "NEW (", def, ");")
            lav$op[def.idx] <- "="
        }
        # replace '==' by '='
        eq.idx <- which(lav$op == "==")
        if(length(eq.idx) > 0L) {
            lav$op[eq.idx] <- "="
        }
        con <- paste(gsub("\\.","",lav$lhs[con.idx]), " ", 
                     lav$op[con.idx], " ", 
                     gsub("\\.","",lav$rhs[con.idx]), ";", sep="")
        con2 <- paste("  ", con, collapse="\n")
        constraints <- paste(constraints, con2, sep="\n") 
    } else {
        constraints <- ""
    }
  
    out <- paste(header, body, constraints, footer, sep="")
    class(out) <- c("psindex.character", "character")
    out
}


# helper functions
lav_mplus_estimator <- function(object) {

    estimator <- object@Options$estimator
    if(estimator == "DWLS") {
        estimator <- "WLS"
    }

    if(estimator == "ML") {
        if(object@Options$test == "yuan.bentler") {
            estimator <- "MLR"
        } else if(object@Options$test == "satorra.bentler") {
            estimator <- "MLM"
        } else if(object@Options$test == "scaled.shifted") {
            estimator <- "MLMV"
        } else if(object@Options$se == "first.order") {
            estimator <- "MLF"
        }
    } else if(estimator %in% c("ULS","WLS")) {
        if(object@Options$test == "satorra.bentler") {
            estimator <- paste(estimator, "M", sep="")
        } else if(object@Options$test == "scaled.shifted") {
            estimator <- paste(estimator, "MV", sep="")
        }
    } else if(estimator == "MML") {
        estimator <- "ML"
    }

    estimator
}

lav_mplus_header <- function(data.file=NULL, group.label="", ov.names="",
                             ov.ord.names="", estimator="ML", 
                             data.type="full", nobs=NULL) {

    # replace '.' by '_' in all variable names
    ov.names     <- gsub("\\.", "_", ov.names)
    ov.ord.names <- gsub("\\.", "_", ov.ord.names)

    ### FIXME!!
    ### this is old code from psindex 0.3-1
    ### surely, this can be done better...

    # TITLE command
    c.TITLE <- "TITLE:\n"
    c.TITLE <- paste(c.TITLE, 
                     "  [This syntax is autogenerated by lavExport]\n")

    # DATA command
    c.DATA <- "DATA:\n"
    ngroups <- length(data.file)    
    if(ngroups == 1L) {
        c.DATA  <- paste(c.DATA,
                         "  file is ", data.file, ";\n", sep="")
    } else {
        for(g in 1:ngroups) {
            c.DATA <- paste(c.DATA,
                            "  file (", group.label[g] ,") is ",
                            data.file[g], ";\n", sep="")
        }
    }
    if(data.type == "full") {
        c.DATA <- paste(c.DATA, "  type is individual;\n", sep="")
    } else if(data.type == "moment") {
        c.DATA <- paste(c.DATA, "  type is fullcov;\n", sep="")
        c.DATA <- paste(c.DATA, "  nobservations are ", nobs, ";\n", sep="")
    } else {
        stop("psindex ERROR: data.type must be full or moment")
    }
    
    # VARIABLE command
    c.VARIABLE <- "VARIABLE:\n"
    c.VARIABLE <- paste(c.VARIABLE, "  names are", sep="")
    nvar <- length(ov.names); tmp <- 0
    for(i in 1:nvar) {
        if(tmp%%6 == 0) { c.VARIABLE <- paste(c.VARIABLE,"\n    ", sep="") }
        c.VARIABLE <- paste(c.VARIABLE, ov.names[i], sep=" ")
        tmp <- tmp+1
    }
    c.VARIABLE <- paste(c.VARIABLE, ";\n", sep="")
    # missing
    if(data.type == "full") {
        c.VARIABLE <- paste(c.VARIABLE,
                            "  missing are all (-999999);\n",sep="")
    }
    # categorical?
    if(length(ov.ord.names)) {
        c.VARIABLE <- paste(c.VARIABLE, "  categorical are", sep="")
        nvar <- length(ov.ord.names); tmp <- 0
        for(i in 1:nvar) {
            if(tmp%%6 == 0) { c.VARIABLE <- paste(c.VARIABLE,"\n    ", sep="") }
            c.VARIABLE <- paste(c.VARIABLE, ov.ord.names[i])
            tmp <- tmp+1
        }
        c.VARIABLE <- paste(c.VARIABLE,";\n",sep="")
    }

    # ANALYSIS command
    c.ANALYSIS <- paste("ANALYSIS:\n  type = general;\n", sep="")
    c.ANALYSIS <- paste(c.ANALYSIS, "  estimator = ", toupper(estimator),
                        ";\n", sep="")

    # MODEL command
    c.MODEL <- paste("MODEL:\n")

    # assemble pre-model header
    out <- paste(c.TITLE, c.DATA, c.VARIABLE, c.ANALYSIS, c.MODEL, sep="")

    out
}