7 Growth factors predicting a distal outcome

7.1 Why this model

In Chapter 5 the growth factors were the outcomes. A time-invariant covariate was used to explain why participants differ in their initial level and in their rate of change. In this chapter the direction of the arrow is reversed. The intercept factor and the slope factor are used as predictors of a variable that is measured later and only once. Such a variable is called a distal outcome.

The question is whether the way in which loneliness develops over the five waves is related to something measured at the end of the study. Here, the distal outcome is substance use at wave 5 (subs5). A model that only uses the loneliness score at wave 5 as a predictor would ignore where participants started and how fast they moved. Here, the intercept and slope factors stand for the whole trajectory and are both used as predictors.

7.2 The idea in plain language

The linear growth model of Chapter 2 describes every participant by an initial level (the intercept factor i) and a rate of change (the slope factor s). Both factors are latent variables, so their values are never observed directly. They are estimated from the five loneliness scores.

In the distal outcome model, one observed variable is added and regressed on i and s. The regression coefficient of i tells us how much higher or lower the outcome is expected to be for participants whose initial level of loneliness is one unit higher, given the same rate of change. The coefficient of s tells us the same for the rate of change, given the same initial level. Thus, the two coefficients are partial effects, as in any multiple regression.

The advantage over a regression on observed scores is that i and s are free of the wave-specific residuals of the loneliness scores.

7.3 A concrete example with the dataset

The dataset has five waves of loneliness (lone1 to lone5) and one measurement of substance use at wave 5 (subs5). We use the time scores 0, 1, 3, 4 and 6 as before, so the intercept refers to the first wave (2015) and the slope is the change per year.

Before fitting the model, we look at the outcome.

Code

table(data$subs5, useNA = "ifany")


   0    1    2    3    4    5    6    7    8   13 <NA> 
 193   77   44   31   17   14    2    2    1    1   60

The variable takes the values 0 to 13 and is observed for 382 of the 442 participants. Of these, 193 have a score of zero, and only 4 have a score above 6. The distribution is therefore strongly right skewed (skewness 2.14, excess kurtosis 7.53). Note that the model treats the outcome as continuous and normal.

7.4 The formal model

The measurement part is the linear growth model. For participant $i$ at wave $t$,

\[ y_{ti} = \eta_{0i} + \lambda_t \, \eta_{1i} + \varepsilon_{ti}, \]

where $\eta_{0i}$ is the intercept factor, $\eta_{1i}$ is the slope factor, and the time scores are $\lambda_t = 0, 1, 3, 4, 6$. The growth factors vary around their means,

\[ \eta_{0i} = \alpha_0 + \zeta_{0i}, \qquad \eta_{1i} = \alpha_1 + \zeta_{1i}, \]

and the two factors are allowed to covary. The distal outcome $z_i$ (here subs5) is added as a regression on both factors,

\[ z_i = \alpha_z + \gamma_0 \, \eta_{0i} + \gamma_1 \, \eta_{1i} + \zeta_{zi}. \]

The parameter $\gamma_0$ is the effect of the initial level and $\gamma_1$ is the effect of the rate of change. The term $\alpha_z$ is the expected outcome when both growth factors are zero, and the variance of $\zeta_{zi}$ is the part of the outcome not accounted for by the growth factors. No covariance between $\zeta_{zi}$ and the residuals $\varepsilon_{ti}$ is estimated in this specification. The diagram in Figure 7.1 shows the model as it was presented in the workshop.

Path diagram of a linear growth model with five repeated measures and a regression of a substance use outcome on the intercept and slope factors — Figure 7.1: The linear growth model with a regression of substance use on the intercept factor I and the slope factor S. The slope loadings 0, 1, 3, 4 and 6 are the time scores, and the two arrows from I and S to the box on the right are the new part. Source: workshop slides (Anyan, 2023).

7.5 R implementation

The model uses the same growth specification as the linear growth model in Chapter 2. The only addition is the regression subs5 ~ i + s. The line subs5 ~ 1 frees the intercept of subs5. Without it, lavaan fixes this intercept to zero, which adds a restriction to the model.

Code

distal <- ' i =~ 1*lone1 + 1*lone2 + 1*lone3 + 1*lone4 + 1*lone5
            s =~ 0*lone1 + 1*lone2 + 3*lone3 + 4*lone4 + 6*lone5

# Time invariant covariates
subs5 ~ i + s

#Estimate the intercept of subs
subs5 ~ 1'
fit_distal <- growth(distal, data = data)
summary(fit_distal, fit.measures = TRUE)

lavaan 0.6.17 ended normally after 62 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        14

                                                  Used       Total
  Number of observations                           290         442

Model Test User Model:
                                                      
  Test statistic                                74.672
  Degrees of freedom                                13
  P-value (Chi-square)                           0.000

Model Test Baseline Model:

  Test statistic                               438.031
  Degrees of freedom                                15
  P-value                                        0.000

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.854
  Tucker-Lewis Index (TLI)                       0.832

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)              -1240.424
  Loglikelihood unrestricted model (H1)      -1203.087
                                                      
  Akaike (AIC)                                2508.847
  Bayesian (BIC)                              2560.226
  Sample-size adjusted Bayesian (SABIC)       2515.829

Root Mean Square Error of Approximation:

  RMSEA                                          0.128
  90 Percent confidence interval - lower         0.101
  90 Percent confidence interval - upper         0.157
  P-value H_0: RMSEA <= 0.050                    0.000
  P-value H_0: RMSEA >= 0.080                    0.998

Standardized Root Mean Square Residual:

  SRMR                                           0.086

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  i =~                                                
    lone1             1.000                           
    lone2             1.000                           
    lone3             1.000                           
    lone4             1.000                           
    lone5             1.000                           
  s =~                                                
    lone1             0.000                           
    lone2             1.000                           
    lone3             3.000                           
    lone4             4.000                           
    lone5             6.000                           

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  subs5 ~                                             
    i                 0.896    0.398    2.251    0.024
    s                 5.889    3.222    1.828    0.068

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  i ~~                                                
    s                -0.015    0.003   -5.793    0.000

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .subs5            -0.006    0.539   -0.011    0.991
    i                 1.466    0.025   58.598    0.000
    s                -0.023    0.004   -5.114    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .lone1             0.055    0.011    5.231    0.000
   .lone2             0.118    0.012    9.727    0.000
   .lone3             0.171    0.015   11.052    0.000
   .lone4             0.122    0.012   10.569    0.000
   .lone5             0.052    0.011    4.929    0.000
   .subs5             2.365    0.200   11.825    0.000
    i                 0.142    0.016    8.925    0.000
    s                 0.003    0.001    5.308    0.000

7.6 Understanding the output

7.6.1 Sample size and missing data

Notice the first lines of the output. The model uses 290 of the 442 rows. As explained in Chapter 1, the workshop script uses listwise deletion. The loneliness scores are complete for 321 participants, and subs5 is missing for another 31 of these, which leaves 290. In total, 152 rows (34%) do not contribute to this model. This is a larger loss than for the linear growth model alone (321 rows).

7.6.2 Fit

The chi-square test of the model is 74.672 on 13 degrees of freedom and is rejected. The fit indices are CFI = .854, TLI = .832, RMSEA = .128 and SRMR = .086. Following the usual combination rules (Hu & Bentler, 1999), these values do not indicate a good fit. The RMSEA is well above the customary limit, and its 90% confidence interval (0.101 to 0.157) does not include values that are conventionally regarded as acceptable. The fit statistics are the same as the ones on the workshop slide in Figure 7.2.

Path diagram and results table for the growth model with a distal outcome, including fit indices and parameter estimates — Figure 7.2: Results of the distal outcome model from the workshop. The slide shows Mplus estimates. The fit statistics and the parameter estimates match the lavaan output in this chapter, and the numbers used in the text come from the R output. Source: workshop slides (Anyan, 2023).

On the slide, “Intercept” and “Slope” under “Distal outcome” are the lines subs5 ~ i and subs5 ~ s in the lavaan output. The intercept of subs5 itself is a different parameter and is not on the slide.

The misfit is mostly inherited from the growth part of the model. We can check this with a short addition that was not part of the workshop. The linear growth model without the distal outcome is fitted to the same 290 participants, and its fit is compared with the fit above.

Code

# Addition, not part of the workshop: linear growth model without subs5,
# fitted to the same participants as the distal outcome model
linear_only <- ' i =~ 1*lone1 + 1*lone2 + 1*lone3 + 1*lone4 + 1*lone5
                 s =~ 0*lone1 + 1*lone2 + 3*lone3 + 4*lone4 + 6*lone5'
fit_lin_same <- growth(linear_only, data = data[cc_all, ])
round(fitMeasures(fit_lin_same, c("chisq", "df", "cfi", "rmsea")), 3)

 chisq     df    cfi  rmsea 
67.954 10.000  0.861  0.141

The linear model alone has a chi-square of 67.954 on 10 degrees of freedom, CFI = .861 and RMSEA = .141. Adding the outcome and the regression adds 3 degrees of freedom and increases the chi-square by about 6.7. Thus, most of the lack of fit comes from the way in which a straight line describes the five loneliness scores, and little of it comes from the regression on the outcome. Remember that the growth factors come from a model that does not fit well, so the coefficients are tentative.

7.6.3 Estimates

The regression coefficients are the main result.

Code

reg[, c("lhs", "op", "rhs", "est", "se", "z", "pvalue", "ci.lower", "ci.upper")]

     lhs op rhs   est    se     z pvalue ci.lower ci.upper
11 subs5  ~   i 0.896 0.398 2.251  0.024    0.116    1.676
12 subs5  ~   s 5.889 3.222 1.828  0.068   -0.426   12.203

The other parameters of interest are the means, variances and covariance of the growth factors, which can be read in the summary above. The mean of the intercept factor is 1.466 and the mean of the slope factor is -0.023 per year. The variances are 0.142 for the intercept factor and 0.003 for the slope factor, and the covariance of the two factors is -0.015. These match the values on the slide, and they are close to those of the linear growth model of Chapter 2 (an intercept of 1.466 and a slope of -0.025 with 321 participants).

7.7 Interpretation

7.7.1 The effect of the initial level

The coefficient of i is 0.896 (standard error 0.398, p = .024). This means that, holding the rate of change constant, a participant whose initial level of loneliness is one unit higher is expected to have a subs5 score that is about 0.9 units higher. The estimate is statistically different from zero at the .05 level, although the 95% confidence interval is wide (from 0.12 to 1.68).

7.7.2 The effect of the rate of change

The coefficient of s is 5.889 (standard error 3.222, p = .068). The estimate is large in absolute size, but the standard error is also large, and the 95% confidence interval runs from -0.43 to 12.20. It includes zero, and the estimate is not statistically different from zero at the .05 level. The data are consistent with a positive association between a steeper increase in loneliness and higher substance use at wave 5, but they do not establish it. The result for the slope also changes with the handling of missing data, as the next section shows.

7.7.3 The scale of the slope coefficient

The size of the slope coefficient depends on the scale of the factor

The slope factor is measured in units of loneliness per year, and its variance is tiny (0.003, which is a standard deviation of 0.056). A difference of one unit in s is therefore far outside the range of values that occur in the data, and the coefficient of 5.889 is large only because it is expressed per unit of a factor that hardly varies. The size of the coefficient alone does not say that the slope is a stronger predictor than the intercept.

To compare the two predictors, we multiply each coefficient by the standard deviation of its factor. This gives the expected difference in subs5 for a participant who is one standard deviation higher on the factor.

Code

# Expected difference in subs5 for a one standard deviation difference in each factor
c(i = gi * sd_lv["i"], s = gs * sd_lv["s"])

      i.i       s.s 
0.3380402 0.3316519

On this scale the two effects are about equal, roughly 0.34 units for the intercept factor and 0.33 units for the slope factor. The standard deviation of subs5 among the 290 participants is 1.56, so each of these differences is about one fifth of a standard deviation of the outcome. The two growth factors are also strongly and negatively correlated (-0.70). Participants who start higher tend to increase less, and this overlap makes it difficult to separate the two coefficients, which is one reason for the large standard error of the slope coefficient.

The model explains only a small part of the variance of the outcome. The R-squared of subs5 is 0.028, which is about 3%. The estimated intercept of subs5 is -0.006. It is the expected outcome for a participant with i and s both equal to zero, which is a combination far from what is observed, and we do not interpret it.

7.8 Sensitivity to the missing data handling

Because only 290 of 442 rows are used, we refit the model with FIML (see Chapter 1). The code below is an addition that was not part of the workshop.

Code

# Addition, not part of the workshop: same model, FIML instead of listwise deletion
fit_distal_fiml <- growth(distal, data = data, missing = "fiml")
fitMeasures(fit_distal_fiml, c("chisq", "df", "cfi", "rmsea", "srmr"))

 chisq     df    cfi  rmsea   srmr 
76.779 13.000  0.888  0.106  0.074

Code

parameterEstimates(fit_distal_fiml) |>
  filter(op == "~") |>
  select(lhs, op, rhs, est, se, z, pvalue)

    lhs op rhs   est    se     z pvalue
1 subs5  ~   i 0.833 0.331 2.514  0.012
2 subs5  ~   s 5.865 2.853 2.056  0.040

lavaan reports that some cases are empty and ignored. These are rows without any observed value on the variables of the model, and they leave 438 rows in the analysis. The coefficient of i is 0.833 (standard error 0.331) and the coefficient of s is 5.865 (standard error 2.853). The point estimates are close to the listwise ones, and the standard errors are smaller. The slope coefficient now has p = .040 and the coefficient of the intercept factor has p = .012. Thus, whether the slope is statistically different from zero depends on the handling of missing data, and a cautious reading is that the evidence for the slope is weak. The fit improves somewhat (CFI = .888, RMSEA = .106) but is still not good.

FIML was found to perform better than listwise deletion in the conditions examined by Enders & Bandalos (2001). It rests on the missing at random assumption, and subs5 is missing for 60 participants, so a relation between missingness and the unobserved level of substance use cannot be ruled out.

7.9 The outcome is treated as continuous and normal

The model is estimated with maximum likelihood, which assumes that the outcome, given the growth factors, is normally distributed. The variable subs5 has a lower bound at zero, about half of the observed values are zero (193 of 382), and the distribution is skewed. The loneliness scores are also skewed, with skewness between 1.52 and 3.05 across the five waves (see Chapter 1). This has two consequences. The ML chi-square tends to be inflated with non-normal data (Curran et al., 1996), so the fit values above may look somewhat worse than they would with a correction. The standard errors also rest on normality, and the predicted values of a linear regression can fall below zero.

A correction is available in lavaan with the estimator "MLR" (Yuan & Bentler, 2000), which gives robust standard errors and a scaled test statistic.

Code

# Addition, not part of the workshop: robust ML with FIML
fit_distal_mlr <- growth(distal, data = data, estimator = "MLR", missing = "fiml")
fitMeasures(fit_distal_mlr, c("chisq.scaled", "df", "cfi.robust", "rmsea.robust"))

chisq.scaled           df   cfi.robust rmsea.robust 
      46.605       13.000        0.896        0.112

With the scaled test statistic, the chi-square is 46.605, CFI = .896 and RMSEA = .112, so the fit remains below the usual standards. A different model for the outcome, such as a count model, is not part of the workshop and is not fitted here.

7.10 What changed compared with the previous models

In Chapter 5 the covariate was a predictor of the growth factors, and the arrows pointed from the observed variable to i and s. In this chapter, the arrows point from i and s to an observed variable. In terms of the lavaan syntax, the growth factors move to the right side of a ~ line. The growth specification, the time scores and the interpretation of the means and variances of the factors are unchanged. What does change is the sample (290 rows instead of 321) and the number of estimated parameters (14), because the outcome adds a regression coefficient for each factor, an intercept and a residual variance.

7.11 What this model answers and what it cannot tell

The model answers the question of whether the initial level and the rate of change of loneliness, estimated from five waves, are statistically associated with substance use at wave 5. In this dataset the association with the initial level is statistically different from zero under listwise deletion, and the association with the slope is uncertain and depends on the handling of missing data. The model also has a poor fit.

Prediction here is statistical, not causal

The outcome subs5 is measured at wave 5, which is also the wave of the last loneliness score lone5. The slope factor is defined by all five waves, so s already contains information from the time at which subs5 is measured. The growth process and the outcome are not separated in time, and the temporal ordering that the word “predict” suggests is not clean. The regression coefficients describe an association, and they do not show that the trajectory of loneliness leads to substance use. The reverse direction, or a shared influence on both variables, is equally compatible with the data.

Several further limits apply. The model assumes a linear trajectory, and the fit shows that a straight line is only a rough description. No covariance between the residual of subs5 and the residual of lone5 is estimated, although both are measured at wave 5, and this restriction could contribute to the misfit. The estimates also come from 290 of 442 participants. Finally, the model says nothing about whether the relation is the same for different groups, which is the topic of Chapter 8.

7.12 Researcher’s perspective

Four choices matter when a distal outcome is added to a growth model: the timing of the outcome, the number of waves, the scale of the factors, and missing data.

An outcome measured after the last wave of the growth process allows a clearer ordering than one measured at the same wave, and the design of the data collection sets that limit. If the outcome is at wave 5, a growth model of waves 1 to 4 would give a cleaner temporal ordering at the cost of less information on the slope. This alternative is not fitted here.

Report the standardized or one-standard-deviation effects next to the raw coefficients, as we did above, so that the reader can judge size. Report the number of rows used, say that listwise deletion was applied, and show a FIML result as a check. When the conclusion changes between the two, as for the slope coefficient here, say so and avoid a strong claim.

Report the fit as it is. The model has CFI = .854 and RMSEA = .128, and the results should be described as tentative. The standard errors and the tests are less trustworthy than they are for a well-fitting model. Groups are compared in Chapter 8.

# Growth factors predicting a distal outcome {#sec-distal} ```{r} #| label: setup #| include: false library(lavaan) library(dplyr) library(tidyr) library(ggplot2) library(psych) data <- read.csv("data/data.csv", fileEncoding = "UTF-8-BOM") theme_set( theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold", color = "black"), axis.title = element_text(color = "black"), axis.text = element_text(color = "black"), strip.text = element_text(color = "black"), legend.text = element_text(color = "black")) ) ``` ## Why this model In @sec-tic the growth factors were the outcomes. A time-invariant covariate was used to explain why participants differ in their initial level and in their rate of change. In this chapter the direction of the arrow is reversed. The intercept factor and the slope factor are used as predictors of a variable that is measured later and only once. Such a variable is called a *distal outcome*. The question is whether the way in which loneliness develops over the five waves is related to something measured at the end of the study. Here, the distal outcome is substance use at wave 5 (`subs5`). A model that only uses the loneliness score at wave 5 as a predictor would ignore where participants started and how fast they moved. Here, the intercept and slope factors stand for the whole trajectory and are both used as predictors. ## The idea in plain language The *linear growth model* of @sec-basic describes every participant by an initial level (the intercept factor `i`) and a rate of change (the slope factor `s`). Both factors are latent variables, so their values are never observed directly. They are estimated from the five loneliness scores. In the distal outcome model, one observed variable is added and regressed on `i` and `s`. The regression coefficient of `i` tells us how much higher or lower the outcome is expected to be for participants whose initial level of loneliness is one unit higher, given the same rate of change. The coefficient of `s` tells us the same for the rate of change, given the same initial level. Thus, the two coefficients are partial effects, as in any multiple regression. The advantage over a regression on observed scores is that `i` and `s` are free of the wave-specific residuals of the loneliness scores. ## A concrete example with the dataset The dataset has five waves of loneliness (`lone1` to `lone5`) and one measurement of substance use at wave 5 (`subs5`). We use the time scores 0, 1, 3, 4 and 6 as before, so the intercept refers to the first wave (2015) and the slope is the change per year. Before fitting the model, we look at the outcome. ```{r} #| label: subs5-desc table(data$subs5, useNA = "ifany") ``` ```{r} #| label: subs5-numbers #| include: false n_total <- nrow(data) n_subs <- sum(!is.na(data$subs5)) n_subs_na <- sum(is.na(data$subs5)) n_zero <- sum(data$subs5 == 0, na.rm = TRUE) n_lone <- sum(complete.cases(data[, paste0("lone", 1:5)])) cc_all <- complete.cases(data[, c(paste0("lone", 1:5), "subs5")]) n_both <- sum(cc_all) sk_subs <- psych::skew(data$subs5, na.rm = TRUE) ku_subs <- psych::kurtosi(data$subs5, na.rm = TRUE) n_gt6 <- sum(data$subs5 > 6, na.rm = TRUE) sk_lone <- sapply(paste0("lone", 1:5), function(v) psych::skew(data[[v]], na.rm = TRUE)) ``` The variable takes the values 0 to 13 and is observed for `r n_subs` of the `r n_total` participants. Of these, `r n_zero` have a score of zero, and only `r n_gt6` have a score above 6. The distribution is therefore strongly right skewed (skewness `r sprintf("%.2f", sk_subs)`, excess kurtosis `r sprintf("%.2f", ku_subs)`). Note that the model treats the outcome as continuous and normal. ## The formal model The measurement part is the linear growth model. For participant $i$ at wave $t$, $$ y_{ti} = \eta_{0i} + \lambda_t \, \eta_{1i} + \varepsilon_{ti}, $$ where $\eta_{0i}$ is the intercept factor, $\eta_{1i}$ is the slope factor, and the time scores are $\lambda_t = 0, 1, 3, 4, 6$. The growth factors vary around their means, $$ \eta_{0i} = \alpha_0 + \zeta_{0i}, \qquad \eta_{1i} = \alpha_1 + \zeta_{1i}, $$ and the two factors are allowed to covary. The distal outcome $z_i$ (here `subs5`) is added as a regression on both factors, $$ z_i = \alpha_z + \gamma_0 \, \eta_{0i} + \gamma_1 \, \eta_{1i} + \zeta_{zi}. $$ The parameter $\gamma_0$ is the effect of the initial level and $\gamma_1$ is the effect of the rate of change. The term $\alpha_z$ is the expected outcome when both growth factors are zero, and the variance of $\zeta_{zi}$ is the part of the outcome not accounted for by the growth factors. No covariance between $\zeta_{zi}$ and the residuals $\varepsilon_{ti}$ is estimated in this specification. The diagram in @fig-distal-diagram shows the model as it was presented in the workshop. ![The linear growth model with a regression of substance use on the intercept factor I and the slope factor S. The slope loadings 0, 1, 3, 4 and 6 are the time scores, and the two arrows from I and S to the box on the right are the new part. Source: workshop slides [@anyan2023workshop].](../images/slides/slide_55.png){#fig-distal-diagram fig-alt="Path diagram of a linear growth model with five repeated measures and a regression of a substance use outcome on the intercept and slope factors" width="90%"} ## R implementation The model uses the same growth specification as the *linear growth model* in @sec-basic. The only addition is the regression `subs5 ~ i + s`. The line `subs5 ~ 1` frees the intercept of `subs5`. Without it, lavaan fixes this intercept to zero, which adds a restriction to the model. ```{r} #| label: distal distal <- ' i =~ 1*lone1 + 1*lone2 + 1*lone3 + 1*lone4 + 1*lone5 s =~ 0*lone1 + 1*lone2 + 3*lone3 + 4*lone4 + 6*lone5 # Time invariant covariates subs5 ~ i + s #Estimate the intercept of subs subs5 ~ 1' fit_distal <- growth(distal, data = data) summary(fit_distal, fit.measures = TRUE) ``` ```{r} #| label: distal-numbers #| include: false fm <- fitMeasures(fit_distal, c("chisq", "df", "cfi", "tli", "rmsea", "srmr")) pe <- parameterEstimates(fit_distal) reg <- pe[pe$op == "~", ] n_used <- lavInspect(fit_distal, "nobs") r2 <- lavInspect(fit_distal, "r2")["subs5"] sd_lv <- sqrt(diag(lavInspect(fit_distal, "cov.lv"))) cor_is <- lavInspect(fit_distal, "cor.lv")["i", "s"] gi <- reg$est[reg$rhs == "i"]; gs <- reg$est[reg$rhs == "s"] fp <- function(p) sub("^0", "", sprintf("%.3f", p)) sd_subs <- sd(data$subs5[cc_all]) rm_ci <- fitMeasures(fit_distal, c("rmsea.ci.lower", "rmsea.ci.upper")) npar <- fitMeasures(fit_distal, "npar") fit_lin_all <- growth(' i =~ 1*lone1 + 1*lone2 + 1*lone3 + 1*lone4 + 1*lone5 s =~ 0*lone1 + 1*lone2 + 3*lone3 + 4*lone4 + 6*lone5', data = data) pl <- parameterEstimates(fit_lin_all) n_lin <- lavInspect(fit_lin_all, "nobs") ``` ## Understanding the output ### Sample size and missing data Notice the first lines of the output. The model uses `r n_used` of the `r n_total` rows. As explained in @sec-prelims, the workshop script uses listwise deletion. The loneliness scores are complete for `r n_lone` participants, and `subs5` is missing for another `r n_lone - n_both` of these, which leaves `r n_both`. In total, `r n_total - n_both` rows (`r sprintf("%.0f", 100 * (n_total - n_both) / n_total)`%) do not contribute to this model. This is a larger loss than for the linear growth model alone (`r n_lin` rows). ### Fit The chi-square test of the model is `r sprintf("%.3f", fm["chisq"])` on `r fm["df"]` degrees of freedom and is rejected. The fit indices are CFI = `r fp(fm["cfi"])`, TLI = `r fp(fm["tli"])`, RMSEA = `r fp(fm["rmsea"])` and SRMR = `r fp(fm["srmr"])`. Following the usual combination rules [@hu1999cutoff], these values do not indicate a good fit. The RMSEA is well above the customary limit, and its 90% confidence interval (`r sprintf("%.3f", rm_ci[1])` to `r sprintf("%.3f", rm_ci[2])`) does not include values that are conventionally regarded as acceptable. The fit statistics are the same as the ones on the workshop slide in @fig-distal-slide. ![Results of the distal outcome model from the workshop. The slide shows Mplus estimates. The fit statistics and the parameter estimates match the lavaan output in this chapter, and the numbers used in the text come from the R output. Source: workshop slides [@anyan2023workshop].](../images/slides/slide_57.png){#fig-distal-slide fig-alt="Path diagram and results table for the growth model with a distal outcome, including fit indices and parameter estimates" width="95%"} On the slide, "Intercept" and "Slope" under "Distal outcome" are the lines `subs5 ~ i` and `subs5 ~ s` in the lavaan output. The intercept of `subs5` itself is a different parameter and is not on the slide. The misfit is mostly inherited from the growth part of the model. We can check this with a short addition that was not part of the workshop. The linear growth model without the distal outcome is fitted to the same `r n_both` participants, and its fit is compared with the fit above. ```{r} #| label: lin-same-n # Addition, not part of the workshop: linear growth model without subs5, # fitted to the same participants as the distal outcome model linear_only <- ' i =~ 1*lone1 + 1*lone2 + 1*lone3 + 1*lone4 + 1*lone5 s =~ 0*lone1 + 1*lone2 + 3*lone3 + 4*lone4 + 6*lone5' fit_lin_same <- growth(linear_only, data = data[cc_all, ]) round(fitMeasures(fit_lin_same, c("chisq", "df", "cfi", "rmsea")), 3) ``` ```{r} #| label: lin-same-numbers #| include: false fl <- fitMeasures(fit_lin_same, c("chisq", "df", "cfi", "rmsea")) ``` The linear model alone has a chi-square of `r sprintf("%.3f", fl["chisq"])` on `r fl["df"]` degrees of freedom, CFI = `r fp(fl["cfi"])` and RMSEA = `r fp(fl["rmsea"])`. Adding the outcome and the regression adds `r fm["df"] - fl["df"]` degrees of freedom and increases the chi-square by about `r sprintf("%.1f", fm["chisq"] - fl["chisq"])`. Thus, most of the lack of fit comes from the way in which a straight line describes the five loneliness scores, and little of it comes from the regression on the outcome. Remember that the growth factors come from a model that does not fit well, so the coefficients are tentative. ### Estimates The regression coefficients are the main result. ```{r} #| label: distal-reg reg[, c("lhs", "op", "rhs", "est", "se", "z", "pvalue", "ci.lower", "ci.upper")] ``` The other parameters of interest are the means, variances and covariance of the growth factors, which can be read in the summary above. The mean of the intercept factor is `r sprintf("%.3f", pe$est[pe$lhs == "i" & pe$op == "~1"])` and the mean of the slope factor is `r sprintf("%.3f", pe$est[pe$lhs == "s" & pe$op == "~1"])` per year. The variances are `r sprintf("%.3f", pe$est[pe$lhs == "i" & pe$op == "~~" & pe$rhs == "i"])` for the intercept factor and `r sprintf("%.3f", pe$est[pe$lhs == "s" & pe$op == "~~" & pe$rhs == "s"])` for the slope factor, and the covariance of the two factors is `r sprintf("%.3f", pe$est[pe$lhs == "i" & pe$op == "~~" & pe$rhs == "s"])`. These match the values on the slide, and they are close to those of the linear growth model of @sec-basic (an intercept of `r sprintf("%.3f", pl$est[pl$lhs == "i" & pl$op == "~1"])` and a slope of `r sprintf("%.3f", pl$est[pl$lhs == "s" & pl$op == "~1"])` with `r n_lin` participants). ## Interpretation ### The effect of the initial level The coefficient of `i` is `r sprintf("%.3f", gi)` (standard error `r sprintf("%.3f", reg$se[reg$rhs == "i"])`, *p* = `r fp(reg$pvalue[reg$rhs == "i"])`). This means that, holding the rate of change constant, a participant whose initial level of loneliness is one unit higher is expected to have a `subs5` score that is about `r sprintf("%.1f", gi)` units higher. The estimate is statistically different from zero at the .05 level, although the 95% confidence interval is wide (from `r sprintf("%.2f", reg$ci.lower[reg$rhs == "i"])` to `r sprintf("%.2f", reg$ci.upper[reg$rhs == "i"])`). ### The effect of the rate of change The coefficient of `s` is `r sprintf("%.3f", gs)` (standard error `r sprintf("%.3f", reg$se[reg$rhs == "s"])`, *p* = `r fp(reg$pvalue[reg$rhs == "s"])`). The estimate is large in absolute size, but the standard error is also large, and the 95% confidence interval runs from `r sprintf("%.2f", reg$ci.lower[reg$rhs == "s"])` to `r sprintf("%.2f", reg$ci.upper[reg$rhs == "s"])`. It includes zero, and the estimate is not statistically different from zero at the .05 level. The data are consistent with a positive association between a steeper increase in loneliness and higher substance use at wave 5, but they do not establish it. The result for the slope also changes with the handling of missing data, as the next section shows. ### The scale of the slope coefficient ::: {.callout-warning title="The size of the slope coefficient depends on the scale of the factor"} The slope factor is measured in units of loneliness per year, and its variance is tiny (`r sprintf("%.3f", pe$est[pe$lhs == "s" & pe$op == "~~" & pe$rhs == "s"])`, which is a standard deviation of `r sprintf("%.3f", sd_lv["s"])`). A difference of one unit in `s` is therefore far outside the range of values that occur in the data, and the coefficient of `r sprintf("%.3f", gs)` is large only because it is expressed per unit of a factor that hardly varies. The size of the coefficient alone does not say that the slope is a stronger predictor than the intercept. ::: To compare the two predictors, we multiply each coefficient by the standard deviation of its factor. This gives the expected difference in `subs5` for a participant who is one standard deviation higher on the factor. ```{r} #| label: scaled-effects # Expected difference in subs5 for a one standard deviation difference in each factor c(i = gi * sd_lv["i"], s = gs * sd_lv["s"]) ``` On this scale the two effects are about equal, roughly `r sprintf("%.2f", gi * sd_lv["i"])` units for the intercept factor and `r sprintf("%.2f", gs * sd_lv["s"])` units for the slope factor. The standard deviation of `subs5` among the `r n_both` participants is `r sprintf("%.2f", sd_subs)`, so each of these differences is about one fifth of a standard deviation of the outcome. The two growth factors are also strongly and negatively correlated (`r sprintf("%.2f", cor_is)`). Participants who start higher tend to increase less, and this overlap makes it difficult to separate the two coefficients, which is one reason for the large standard error of the slope coefficient. The model explains only a small part of the variance of the outcome. The *R*-squared of `subs5` is `r sprintf("%.3f", r2)`, which is about `r sprintf("%.0f", 100 * r2)`%. The estimated intercept of `subs5` is `r sprintf("%.3f", pe$est[pe$lhs == "subs5" & pe$op == "~1"])`. It is the expected outcome for a participant with `i` and `s` both equal to zero, which is a combination far from what is observed, and we do not interpret it. ## Sensitivity to the missing data handling Because only `r n_both` of `r n_total` rows are used, we refit the model with FIML (see @sec-prelims). The code below is an addition that was not part of the workshop. ```{r} #| label: distal-fiml # Addition, not part of the workshop: same model, FIML instead of listwise deletion fit_distal_fiml <- growth(distal, data = data, missing = "fiml") fitMeasures(fit_distal_fiml, c("chisq", "df", "cfi", "rmsea", "srmr")) parameterEstimates(fit_distal_fiml) |> filter(op == "~") |> select(lhs, op, rhs, est, se, z, pvalue) ``` ```{r} #| label: fiml-numbers #| include: false fmf <- fitMeasures(fit_distal_fiml, c("chisq", "cfi", "rmsea")) pef <- parameterEstimates(fit_distal_fiml) regf <- pef[pef$op == "~", ] n_fiml <- lavInspect(fit_distal_fiml, "nobs") ``` lavaan reports that some cases are empty and ignored. These are rows without any observed value on the variables of the model, and they leave `r n_fiml` rows in the analysis. The coefficient of `i` is `r sprintf("%.3f", regf$est[regf$rhs == "i"])` (standard error `r sprintf("%.3f", regf$se[regf$rhs == "i"])`) and the coefficient of `s` is `r sprintf("%.3f", regf$est[regf$rhs == "s"])` (standard error `r sprintf("%.3f", regf$se[regf$rhs == "s"])`). The point estimates are close to the listwise ones, and the standard errors are smaller. The slope coefficient now has *p* = `r fp(regf$pvalue[regf$rhs == "s"])` and the coefficient of the intercept factor has *p* = `r fp(regf$pvalue[regf$rhs == "i"])`. Thus, whether the slope is statistically different from zero depends on the handling of missing data, and a cautious reading is that the evidence for the slope is weak. The fit improves somewhat (CFI = `r fp(fmf["cfi"])`, RMSEA = `r fp(fmf["rmsea"])`) but is still not good. FIML was found to perform better than listwise deletion in the conditions examined by @enders2001relative. It rests on the missing at random assumption, and `subs5` is missing for `r n_subs_na` participants, so a relation between missingness and the unobserved level of substance use cannot be ruled out. ## The outcome is treated as continuous and normal The model is estimated with maximum likelihood, which assumes that the outcome, given the growth factors, is normally distributed. The variable `subs5` has a lower bound at zero, about half of the observed values are zero (`r n_zero` of `r n_subs`), and the distribution is skewed. The loneliness scores are also skewed, with skewness between `r sprintf("%.2f", min(sk_lone))` and `r sprintf("%.2f", max(sk_lone))` across the five waves (see @sec-prelims). This has two consequences. The ML chi-square tends to be inflated with non-normal data [@curran1996robustness], so the fit values above may look somewhat worse than they would with a correction. The standard errors also rest on normality, and the predicted values of a linear regression can fall below zero. A correction is available in lavaan with the estimator `"MLR"` [@yuan2000three], which gives robust standard errors and a scaled test statistic. ```{r} #| label: distal-mlr # Addition, not part of the workshop: robust ML with FIML fit_distal_mlr <- growth(distal, data = data, estimator = "MLR", missing = "fiml") fitMeasures(fit_distal_mlr, c("chisq.scaled", "df", "cfi.robust", "rmsea.robust")) ``` ```{r} #| label: mlr-numbers #| include: false fmr <- fitMeasures(fit_distal_mlr, c("chisq.scaled", "cfi.robust", "rmsea.robust")) ``` With the scaled test statistic, the chi-square is `r sprintf("%.3f", fmr["chisq.scaled"])`, CFI = `r fp(fmr["cfi.robust"])` and RMSEA = `r fp(fmr["rmsea.robust"])`, so the fit remains below the usual standards. A different model for the outcome, such as a count model, is not part of the workshop and is not fitted here. ## What changed compared with the previous models In @sec-tic the covariate was a predictor of the growth factors, and the arrows pointed from the observed variable to `i` and `s`. In this chapter, the arrows point from `i` and `s` to an observed variable. In terms of the lavaan syntax, the growth factors move to the right side of a `~` line. The growth specification, the time scores and the interpretation of the means and variances of the factors are unchanged. What does change is the sample (`r n_used` rows instead of `r n_lin`) and the number of estimated parameters (`r npar`), because the outcome adds a regression coefficient for each factor, an intercept and a residual variance. ## What this model answers and what it cannot tell The model answers the question of whether the initial level and the rate of change of loneliness, estimated from five waves, are statistically associated with substance use at wave 5. In this dataset the association with the initial level is statistically different from zero under listwise deletion, and the association with the slope is uncertain and depends on the handling of missing data. The model also has a poor fit. ::: {.callout-important title="Prediction here is statistical, not causal"} The outcome `subs5` is measured at wave 5, which is also the wave of the last loneliness score `lone5`. The slope factor is defined by all five waves, so `s` already contains information from the time at which `subs5` is measured. The growth process and the outcome are not separated in time, and the temporal ordering that the word "predict" suggests is not clean. The regression coefficients describe an association, and they do not show that the trajectory of loneliness leads to substance use. The reverse direction, or a shared influence on both variables, is equally compatible with the data. ::: Several further limits apply. The model assumes a linear trajectory, and the fit shows that a straight line is only a rough description. No covariance between the residual of `subs5` and the residual of `lone5` is estimated, although both are measured at wave 5, and this restriction could contribute to the misfit. The estimates also come from `r n_both` of `r n_total` participants. Finally, the model says nothing about whether the relation is the same for different groups, which is the topic of @sec-multigroup. ## Researcher's perspective Four choices matter when a distal outcome is added to a growth model: the timing of the outcome, the number of waves, the scale of the factors, and missing data. An outcome measured after the last wave of the growth process allows a clearer ordering than one measured at the same wave, and the design of the data collection sets that limit. If the outcome is at wave 5, a growth model of waves 1 to 4 would give a cleaner temporal ordering at the cost of less information on the slope. This alternative is not fitted here. Report the standardized or one-standard-deviation effects next to the raw coefficients, as we did above, so that the reader can judge size. Report the number of rows used, say that listwise deletion was applied, and show a FIML result as a check. When the conclusion changes between the two, as for the slope coefficient here, say so and avoid a strong claim. Report the fit as it is. The model has CFI = `r fp(fm["cfi"])` and RMSEA = `r fp(fm["rmsea"])`, and the results should be described as tentative. The standard errors and the tests are less trustworthy than they are for a well-fitting model. Groups are compared in @sec-multigroup.