R/lav_sam_step2.R

In lavaan: Latent Variable Analysis

# SAM step 2: estimate structural part

lav_sam_step2 <- function(STEP1 = NULL, FIT = NULL,
                          sam.method = "local", struc.args = list()) {

    lavoptions <- FIT@Options
    lavpta     <- FIT@pta
    nlevels    <- lavpta$nlevels
    PT         <- STEP1$PT

    LV.names <- unique(unlist(FIT@pta$vnames$lv.regular))

    # adjust options
    lavoptions.PA <- lavoptions
    if(lavoptions.PA$se == "naive") {
        lavoptions.PA$se <- "standard"
    } else {
        # twostep or none -> none
        lavoptions.PA$se <- "none"
    }
    #lavoptions.PA$fixed.x <- TRUE # may be false if indicator is predictor
    lavoptions.PA$fixed.x <- FALSE # until we fix this...
    lavoptions.PA$verbose <- FALSE # must be in struc.args
    lavoptions.PA$categorical <- FALSE
    lavoptions.PA$.categorical <- FALSE
    lavoptions.PA <- modifyList(lavoptions.PA, struc.args)

    # override, no matter what
    lavoptions.PA$do.fit <- TRUE

    if(sam.method %in% c("local", "fsr")) {
        lavoptions.PA$missing <- "listwise"
        lavoptions.PA$sample.cov.rescale <- FALSE
        #lavoptions.PA$baseline <- FALSE
        lavoptions.PA$h1 <- FALSE
        #lavoptions.PA$implied <- FALSE
        lavoptions.PA$loglik <- FALSE
    } else {
        lavoptions.PA$h1 <- FALSE
        #lavoptions.PA$implied <- FALSE
        lavoptions.PA$loglik <- FALSE
    }


    # construct PTS
    if(sam.method %in% c("local", "fsr")) {
        # extract structural part
        PTS <- lav_partable_subset_structural_model(PT, lavpta = lavpta,
                   add.idx = TRUE,
                   fixed.x = lavoptions.PA$fixed.x,
                   add.exo.cov = TRUE)
        # remove est/se/start columns
        PTS$est   <- NULL
        PTS$se    <- NULL
        PTS$start <- NULL

        if(nlevels > 1L) {
            PTS$level <- NULL
            PTS$group <- NULL
            PTS$group <- PTS$block
            NOBS <- FIT@Data@Lp[[1]]$nclusters
        } else {
            NOBS <- FIT@Data@nobs
        }
        # if meanstructure, 'free' user=0 intercepts?
        if(lavoptions.PA$meanstructure) {
            extra.int.idx <- which(PTS$op == "~1" & PTS$user == 0L &
                                   PTS$free == 0L &
                                   PTS$exo == 0L) # needed?
            if(length(extra.int.idx) > 0L) {
                PTS$free[  extra.int.idx ] <- 1L
                PTS$ustart[extra.int.idx ] <- as.numeric(NA)
                PTS$free[ PTS$free > 0L ] <-
                    seq_len( length(PTS$free[ PTS$free > 0L ]) )
                PTS$user[ extra.int.idx ] <- 3L
            }
        } else {
            extra.int.idx <- integer(0L)
        }
        reg.idx <- attr(PTS, "idx"); attr(PTS, "idx") <- NULL
    } else {

        # global SAM

        # the measurement model parameters now become fixed ustart values
        PT$ustart[PT$free > 0] <- PT$est[PT$free > 0]

        reg.idx <- lav_partable_subset_structural_model(PT = PT,
                          lavpta = lavpta, idx.only = TRUE)

        # remove 'exogenous' factor variances (if any) from reg.idx
        lv.names.x <- LV.names[ LV.names %in% unlist(lavpta$vnames$eqs.x)  &
                               !LV.names %in% unlist(lavpta$vnames$eqs.y) ]
        if(lavoptions.PA$fixed.x && length(lv.names.x) > 0L) {
            var.idx <- which(PT$lhs %in% lv.names.x &
                             PT$op == "~~" &
                             PT$lhs == PT$rhs)
            rm.idx <- which(reg.idx %in% var.idx)
            if(length(rm.idx) > 0L) {
                reg.idx <- reg.idx[ -rm.idx ]
            }
        }

        # adapt parameter table for structural part
        PTS <- PT

        # remove constraints we don't need
        con.idx <- which(PTS$op %in% c("==","<",">",":="))
        if(length(con.idx) > 0L) {
            needed.idx <- which(con.idx %in% reg.idx)
            if(length(needed.idx) > 0L) {
                con.idx <- con.idx[-needed.idx]
            }
            if(length(con.idx) > 0L) {
                PTS <- as.data.frame(PTS, stringsAsFactors = FALSE)
                PTS <- PTS[-con.idx, ]
            }
        }
        PTS$est <- NULL
        PTS$se <- NULL

        PTS$free[ !seq_len(length(PTS$lhs)) %in% reg.idx & PTS$free > 0L ] <- 0L
        PTS$free[ PTS$free > 0L ] <- seq_len( sum(PTS$free > 0L) )

        # set 'ustart' values for free FIT.PA parameter to NA
        PTS$ustart[ PTS$free > 0L ] <- as.numeric(NA)

        extra.int.idx <- integer(0L)
    } # global

   # fit structural model
    if(lavoptions$verbose) {
        cat("Fitting the structural part ... \n")
    }
    if(sam.method %in% c("local", "fsr")) {
        FIT.PA <- lavaan::lavaan(PTS,
                                 sample.cov = STEP1$VETA,
                                 sample.mean = STEP1$EETA,
                                 sample.nobs = NOBS,
                                 slotOptions = lavoptions.PA)

    } else {
        FIT.PA <- lavaan::lavaan(model = PTS,
                                 slotData = FIT@Data,
                                 slotSampleStats = FIT@SampleStats,
                                 slotOptions = lavoptions.PA)
    }
    if(lavoptions$verbose) {
        cat("Fitting the structural part ... done.\n")
    }

    # which parameters from PTS do we wish to fill in:
    # - all 'free' parameters
    # - :=, <, > (if any)
    # - and NOT element with user=3 (add.exo.cov = TRUE, extra.int.idx)
    pts.idx <- which( (PTS$free > 0L | (PTS$op %in% c(":=", "<", ">"))) &
                      !PTS$user == 3L )

    # find corresponding rows in PT
    PTS2 <- as.data.frame(PTS, stringsAsFactors = FALSE)
    pt.idx <- lav_partable_map_id_p1_in_p2(PTS2[pts.idx,], PT,
                                           exclude.nonpar = FALSE)
    # fill in
    PT$est[ pt.idx ] <- FIT.PA@ParTable$est[ pts.idx ]

    # create step2.free.idx
    p2.idx <- seq_len(length(PT$lhs)) %in% pt.idx & PT$free > 0 # no def!
    step2.free.idx <- STEP1$PT.free[ p2.idx ]

    # add 'step' column in PT
    PT$step <- rep(1L, length(PT$lhs))
    PT$step[seq_len(length(PT$lhs)) %in% reg.idx] <- 2L

    STEP2 <- list(FIT.PA = FIT.PA, PT = PT, reg.idx = reg.idx,
                  step2.free.idx = step2.free.idx, extra.int.idx = extra.int.idx)

    STEP2
}