# Examples of Longitudinal Models with Time-varying Covariates

## Load nlpsem package, dependent packages and set CSOLNP as the optimizer

    library(nlpsem)
    mxOption(model = NULL, key = "Default optimizer", "CSOLNP", reset = FALSE)

## Load pre-computed models

    load(system.file("extdata", "getTVCmodel_examples.RData", package = "nlpsem"))

## Load example data and preprocess data

    # Load ECLS-K (2011) data
    data("RMS_dat")
    RMS_dat0 <- RMS_dat
    # Re-baseline the data so that the estimated initial status is for the
    # starting point of the study
    baseT <- RMS_dat0$T1
    RMS_dat0$T1 <- (RMS_dat0$T1 - baseT)/12
    RMS_dat0$T2 <- (RMS_dat0$T2 - baseT)/12
    RMS_dat0$T3 <- (RMS_dat0$T3 - baseT)/12
    RMS_dat0$T4 <- (RMS_dat0$T4 - baseT)/12
    RMS_dat0$T5 <- (RMS_dat0$T5 - baseT)/12
    RMS_dat0$T6 <- (RMS_dat0$T6 - baseT)/12
    RMS_dat0$T7 <- (RMS_dat0$T7 - baseT)/12
    RMS_dat0$T8 <- (RMS_dat0$T8 - baseT)/12
    RMS_dat0$T9 <- (RMS_dat0$T9 - baseT)/12
    # Standardize time-invariant covariates (TICs)
    ## ex1 is standardized growth TIC in models
    RMS_dat0$ex1 <- scale(RMS_dat0$Approach_to_Learning)
    # Standardize time-varying covariate (TVC)
    BL_mean <- mean(RMS_dat0[, "R1"])
    BL_var <- var(RMS_dat0[, "R1"])
    RMS_dat0$Rs1 <- (RMS_dat0$R1 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs2 <- (RMS_dat0$R2 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs3 <- (RMS_dat0$R3 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs4 <- (RMS_dat0$R4 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs5 <- (RMS_dat0$R5 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs6 <- (RMS_dat0$R6 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs7 <- (RMS_dat0$R7 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs8 <- (RMS_dat0$R8 - BL_mean)/sqrt(BL_var)
    RMS_dat0$Rs9 <- (RMS_dat0$R9 - BL_mean)/sqrt(BL_var)
    xstarts <- mean(baseT)

## Example 1: This example includes two models. Model 1 is a full bilinear spline LGCM with a TVC to examine the influence of baseline teacher-reported approach to learning and the development in reading ability on the development of mathematics ability. It also includes a visualization showcasing the growth status of mathematics ability. Model 2 is a full bilinear spline LGCM with a decomposed TVC (interval-specific slopes) to examine the influence of baseline teacher-reported approach to learning and the development in reading ability on the development of mathematics ability. P values and Wald confidence intervals of all parameters are provided. It also includes a visualization showcasing the growth status of mathematics ability.

    set.seed(20191029)
    Math_TVC_BLS_f <- getTVCmodel(
      dat = RMS_dat0, t_var = "T", y_var = "M", curveFun = "BLS", intrinsic = TRUE,
      records = 1:9, y_model = "LGCM", TVC = "Rs", decompose = 0, growth_TIC = "ex1",
      res_scale = 0.1, tries = 10
      )
    paraBLS_TVC.f <- c(
      "Y_alpha0", "Y_alpha1", "Y_alpha2", "Y_alphag",
      paste0("Y_psi", c("00", "01", "02", "0g", "11", "12", "1g", "22", "2g", "gg")),
      "Y_residuals", "X_mueta0", "X_mueta1", paste0("X_psi", c("00", "01", "11")),
      paste0("X_rel_rate", 2:8), paste0("X_abs_rate", 1:8), "X_residuals",
      paste0("betaTIC", c(0:2, "g")), paste0("betaTVC", c(0:2, "g")), "muTIC", "phiTIC",
      "Y_mueta0", "Y_mueta1", "Y_mueta2", "Y_mu_knot", "covBL", "kappa", "Cov_XYres"
      )
    set.seed(20191029)
    Math_TVCslp_BLS_f <- getTVCmodel(
      dat = RMS_dat0, t_var = "T", y_var = "M", curveFun = "BLS", intrinsic = TRUE,
      records = 1:9, y_model = "LGCM", TVC = "Rs", decompose = 1, growth_TIC = "ex1",
      res_scale = c(0.1, 0.1), res_cor = 0.3, tries = 10, paramOut = TRUE,
      names = paraBLS_TVC.f
      )

    getEstimateStats(est_in = Math_TVCslp_BLS_f@Estimates, CI_type = "Wald")
    #> An object of class "StatsOutput"
    #> Slot "wald":
    #>             Estimate     SE p.value wald_lbound wald_ubound
    #> Y_alpha0     37.1158 0.4057 <0.0001     36.3206     37.9110
    #> Y_alpha1     20.4875 0.2061 <0.0001     20.0836     20.8914
    #> Y_alpha2      8.2065 0.1880 <0.0001      7.8380      8.5750
    #> Y_alphag      3.1600 0.0419 <0.0001      3.0779      3.2421
    #> Y_psi00      47.2135 4.1854 <0.0001     39.0103     55.4167
    #> Y_psi01      -1.2051 1.6321  0.4603     -4.4040      1.9938
    #> Y_psi02      -2.2946 1.4974  0.1254     -5.2295      0.6403
    #> Y_psi0g      -0.3718 0.3036  0.2207     -0.9668      0.2232
    #> Y_psi11      12.1669 1.2113 <0.0001      9.7928     14.5410
    #> Y_psi12      -0.2187 0.8047  0.7858     -1.7959      1.3585
    #> Y_psi1g      -0.7959 0.1918 <0.0001     -1.1718     -0.4200
    #> Y_psi22       3.3249 1.1547   0.004      1.0617      5.5881
    #> Y_psi2g      -0.3249 0.1806   0.072     -0.6789      0.0291
    #> Y_psigg       0.2939 0.0484 <0.0001      0.1990      0.3888
    #> Y_residuals  28.3922 0.8050 <0.0001     26.8144     29.9700
    #> X_mueta0      0.0093 0.0556  0.8672     -0.0997      0.1183
    #> X_mueta1      2.5036 0.0701 <0.0001      2.3662      2.6410
    #> X_psi00       1.2124 0.0857 <0.0001      1.0444      1.3804
    #> X_psi01      -0.0776 0.0203   1e-04     -0.1174     -0.0378
    #> X_psi11       0.0980 0.0106 <0.0001      0.0772      0.1188
    #> X_rel_rate2   0.6947 0.0384 <0.0001      0.6194      0.7700
    #> X_rel_rate3   1.2478 0.0431 <0.0001      1.1633      1.3323
    #> X_rel_rate4   0.4903 0.0303 <0.0001      0.4309      0.5497
    #> X_rel_rate5   0.7110 0.0330 <0.0001      0.6463      0.7757
    #> X_rel_rate6   0.2823 0.0169 <0.0001      0.2492      0.3154
    #> X_rel_rate7   0.3019 0.0168 <0.0001      0.2690      0.3348
    #> X_rel_rate8   0.2568 0.0165 <0.0001      0.2245      0.2891
    #> X_abs_rate1   2.5036 0.0701 <0.0001      2.3662      2.6410
    #> X_abs_rate2   1.7393 0.0664 <0.0001      1.6092      1.8694
    #> X_abs_rate3   3.1242 0.0664 <0.0001      2.9941      3.2543
    #> X_abs_rate4   1.2275 0.0680 <0.0001      1.0942      1.3608
    #> X_abs_rate5   1.7801 0.0657 <0.0001      1.6513      1.9089
    #> X_abs_rate6   0.7069 0.0380 <0.0001      0.6324      0.7814
    #> X_abs_rate7   0.7559 0.0367 <0.0001      0.6840      0.8278
    #> X_abs_rate8   0.6429 0.0375 <0.0001      0.5694      0.7164
    #> X_residuals   0.3587 0.0088 <0.0001      0.3415      0.3759
    #> betaTIC0      0.9665 0.3921  0.0137      0.1980      1.7350
    #> betaTIC1      0.1190 0.2145   0.579     -0.3014      0.5394
    #> betaTIC2     -0.2104 0.2004  0.2938     -0.6032      0.1824
    #> betaTICg      0.0049 0.0405  0.9037     -0.0745      0.0843
    #> betaTVC0      7.6230 0.3801 <0.0001      6.8780      8.3680
    #> betaTVC1      0.8498 0.2045 <0.0001      0.4490      1.2506
    #> betaTVC2     -0.3467 0.1912  0.0698     -0.7214      0.0280
    #> betaTVCg     -0.1272 0.0394  0.0012     -0.2044     -0.0500
    #> muTIC         0.0000 0.0447 >0.9999     -0.0876      0.0876
    #> phiTIC        0.9980 0.0631 <0.0001      0.8743      1.1217
    #> Y_mueta0     37.1867 0.5451 <0.0001     36.1183     38.2551
    #> Y_mueta1     20.4954 0.2062 <0.0001     20.0913     20.8995
    #> Y_mueta2      8.2032 0.1890 <0.0001      7.8328      8.5736
    #> Y_mu_knot     3.1588 0.0418 <0.0001      3.0769      3.2407
    #> covBL         0.4000 0.0523 <0.0001      0.2975      0.5025
    #> kappa         1.7789 0.1059 <0.0001      1.5713      1.9865
    #> Cov_XYres     0.5510 0.0619 <0.0001      0.4297      0.6723
    #>
    #> Slot "likelihood":
    #> data frame with 0 columns and 0 rows
    #>
    #> Slot "bootstrap":
    #> data frame with 0 columns and 0 rows
    Figure1 <- getFigure(
      model = Math_TVC_BLS_f@mxOutput, sub_Model = "TVC", y_var = "M", curveFun = "BLS",
      y_model = "LGCM", t_var = "T", records = 1:9, xstarts = xstarts, xlab = "Year",
      outcome = "Mathematics"
    )
    #> Treating first argument as an object that stores a character
    show(Figure1)
    #> figOutput Object
    #> --------------------
    #> Trajectories: 1
    #> Figure 1:
    #> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

    Figure2 <- getFigure(
      model = Math_TVCslp_BLS_f@mxOutput, sub_Model = "TVC", y_var = "M", curveFun = "BLS",
      y_model = "LGCM", t_var = "T", records = 1:9, xstarts = xstarts, xlab = "Year",
      outcome = "Mathematics"
    )
    #> Treating first argument as an object that stores a character
    show(Figure2)
    #> figOutput Object
    #> --------------------
    #> Trajectories: 1
    #> Figure 1:
    #> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

A comparison between Figure 1 and Figure 2 demonstrates that incorporating a TVC directly results in underestimation of the growth factor means.

## Example 2: Fit reduced bilinear spline LGCMs with a decomposed TVC (interval-specific slopes, interval-specific changes, and change from baseline) to examine the influence of baseline teacher-reported approach to learning and the development in reading ability on the development of mathematics ability. It also includes a visualization showcasing the growth status of mathematics ability.

    paraBLS_TVC.r <- c(
      "Y_alpha0", "Y_alpha1", "Y_alpha2", "Y_knot",
      paste0("Y_psi", c("00", "01", "02", "11", "12", "22")), "Y_residuals",
      "X_mueta0", "X_mueta1", paste0("X_psi", c("00", "01", "11")),
      paste0("X_rel_rate", 2:8), paste0("X_abs_rate", 1:8), "X_residuals",
      paste0("betaTIC", 0:2), paste0("betaTVC", 0:2), "muTIC", "phiTIC",
      "Y_mueta0", "Y_mueta1", "Y_mueta2", "covBL", "kappa", "Cov_XYres"
      )
    set.seed(20191029)
    Math_TVCslp_BLS_r <- getTVCmodel(
      dat = RMS_dat0, t_var = "T", y_var = "M", curveFun = "BLS", intrinsic = FALSE,
      records = 1:9, y_model = "LGCM", TVC = "R", decompose = 1, growth_TIC = "ex1",
      res_scale = c(0.1, 0.1), res_cor = 0.3, tries = 10,  paramOut = TRUE,
      names = paraBLS_TVC.r)
    set.seed(20191029)
    Math_TVCchg_BLS_r <- getTVCmodel(
      dat = RMS_dat0, t_var = "T", y_var = "M", curveFun = "BLS", intrinsic = FALSE,
      records = 1:9, y_model = "LGCM", TVC = "R", decompose = 2, growth_TIC = "ex1",
      res_scale = c(0.1, 0.1), res_cor = 0.3, tries = 10,  paramOut = TRUE,
      names = paraBLS_TVC.r)
    set.seed(20191029)
    Math_TVCchgBL_BLS_r <- getTVCmodel(
      dat = RMS_dat0, t_var = "T", y_var = "M", curveFun = "BLS", intrinsic = FALSE,
      records = 1:9, y_model = "LGCM", TVC = "R", decompose = 3, growth_TIC = "ex1",
      res_scale = c(0.1, 0.1), res_cor = 0.3, tries = 10,  paramOut = TRUE,
      names = paraBLS_TVC.r)

    Figure3 <- getFigure(
      model = Math_TVCslp_BLS_r@mxOutput, sub_Model = "TVC", y_var = "M", curveFun = "BLS",
      y_model = "LGCM", t_var = "T", records = 1:9, xstarts = xstarts, xlab = "Year",
      outcome = "Mathematics"
    )
    #> Treating first argument as an object that stores a character
    #> Treating first argument as an object that stores a character
    #> Treating first argument as an object that stores a character
    show(Figure3)
    #> figOutput Object
    #> --------------------
    #> Trajectories: 1
    #> Figure 1:
    #> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

    Figure4 <- getFigure(
      model = Math_TVCchg_BLS_r@mxOutput, sub_Model = "TVC", y_var = "M", curveFun = "BLS",
      y_model = "LGCM", t_var = "T", records = 1:9, xstarts = xstarts, xlab = "Year",
      outcome = "Mathematics"
    )
    #> Treating first argument as an object that stores a character
    #> Treating first argument as an object that stores a character
    #> Treating first argument as an object that stores a character
    show(Figure4)
    #> figOutput Object
    #> --------------------
    #> Trajectories: 1
    #> Figure 1:
    #> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

    Figure5 <- getFigure(
      model = Math_TVCchgBL_BLS_r@mxOutput, sub_Model = "TVC", y_var = "M", curveFun = "BLS",
      y_model = "LGCM", t_var = "T", records = 1:9, xstarts = xstarts, xlab = "Year",
      outcome = "Mathematics"
    )
    #> Treating first argument as an object that stores a character
    #> Treating first argument as an object that stores a character
    #> Treating first argument as an object that stores a character
    show(Figure5)
    #> figOutput Object
    #> --------------------
    #> Trajectories: 1
    #> Figure 1:
    #> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'