In a previous post, we covered how to use Latent Growth Modeling in R to examine change in time. In that post, we assumed a simple linear model, which is often unrealistic. Here, I am going to show how we can free this assumption and find the best way to treat change in time.
We can use exploratory analysis and previous research to understand how to model change in time. Moreover, we can also compare models that treat change in time in different ways to find the best fit for our data.
Want to really learn? Follow along using the data and code.
Get them when you join the newsletter. Additionally, you will get my 63-point pre-submission checklist for your research. All for free!
If we look again at the log income in the Understanding Society data, we get this graph (see previous post for an explanation of the data and how it is structured):
set.seed(20260803)
plot_ids <- sample(unique(usl$pidp), 2000)
observed_means <- usl |>
group_by(wave) |>
summarise(logincome = mean(logincome), .groups = "drop")
usl |>
filter(pidp %in% plot_ids) |>
ggplot(aes(wave, logincome, group = pidp)) +
geom_line(alpha = 0.01) +
geom_line(
data = observed_means,
aes(wave, logincome),
inherit.aes = FALSE,
linewidth = 1.5,
colour = "red"
) +
theme_bw() +
labs(x = "Wave", y = "Log income")
The graph would indicate we have an average change that is overall linear with a slight downward bend.
To see the individual level change, I also sampled 20 individuals and plotted each one’s change in income.
set.seed(20260803) people <- sample(unique(usl$pidp), 20) usl |> filter(pidp %in% people) |> ggplot(aes(wave, logincome, group = 1)) + geom_line() + facet_wrap(~pidp) + theme_bw() + labs(x = "Wave", y = "Log income")
It appears that at the individual level, we have a more mixed bag, although linear change would not be a bad approximation for quite a few people.
Keeping this in mind, we can also decide on the best way to model change in time by comparing several different models and seeing which fits the data best. As a starting point, we can run the linear model, which will be our reference (see previous post for an explanation of the model and syntax):
model <- 'i =~ 1*logincome_1 + 1*logincome_2 + 1*logincome_3 +
1*logincome_4 + 1*logincome_5 + 1*logincome_6
s =~ 0*logincome_1 + 1*logincome_2 + 2*logincome_3 +
3*logincome_4 + 4*logincome_5 + 5*logincome_6'
fit1 <- growth(model, data = usw)
summary(fit1, standardized = TRUE)
## lavaan 0.6-19 ended normally after 40 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 11
##
## Number of observations 21178
##
## Model Test User Model:
##
## Test statistic 1122.274
## Degrees of freedom 16
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i =~
## logincome_1 1.000 1.121 0.791
## logincome_2 1.000 1.121 0.908
## logincome_3 1.000 1.121 0.955
## logincome_4 1.000 1.121 0.970
## logincome_5 1.000 1.121 0.981
## logincome_6 1.000 1.121 1.003
## s =~
## logincome_1 0.000 0.000 0.000
## logincome_2 1.000 0.166 0.135
## logincome_3 2.000 0.333 0.284
## logincome_4 3.000 0.499 0.432
## logincome_5 4.000 0.666 0.583
## logincome_6 5.000 0.832 0.745
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i ~~
## s -0.116 0.003 -44.994 0.000 -0.622 -0.622
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i 6.885 0.009 801.397 0.000 6.144 6.144
## s 0.056 0.002 33.268 0.000 0.334 0.334
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .logincome_1 0.751 0.010 74.715 0.000 0.751 0.374
## .logincome_2 0.473 0.006 74.687 0.000 0.473 0.310
## .logincome_3 0.473 0.006 83.529 0.000 0.473 0.344
## .logincome_4 0.527 0.006 87.255 0.000 0.527 0.394
## .logincome_5 0.533 0.006 83.182 0.000 0.533 0.409
## .logincome_6 0.460 0.007 64.709 0.000 0.460 0.368
## i 1.256 0.015 81.824 0.000 1.000 1.000
## s 0.028 0.001 44.340 0.000 1.000 1.000The model estimates an average initial log income of 6.885 (about £978) and an average linear increase of 0.056 per wave. The interpretation is the same as before, although the numerical values change with the synthetic data.
We can also visualize the change in time based on our model (again, check the previous post for explanations).
pred_lgm <- predict(fit1)
pred_lgm_long <- map(
0:5,
function(x) pred_lgm[, 1] + x * pred_lgm[, 2]
) |>
reduce(cbind) |>
as.data.frame() |>
setNames(str_c("Wave ", 1:6)) |>
mutate(id = row_number()) |>
pivot_longer(-id, names_to = "wave", values_to = "pred")
set.seed(20260803)
prediction_ids <- sample(unique(pred_lgm_long$id), 2000)
linear_means <- pred_lgm_long |>
group_by(wave) |>
summarise(pred = mean(pred), .groups = "drop")
pred_lgm_long |>
filter(id %in% prediction_ids) |>
ggplot(aes(wave, pred, group = id)) +
geom_line(alpha = 0.01) +
geom_line(
data = linear_means,
aes(wave, pred, group = 1),
inherit.aes = FALSE,
linewidth = 1.5,
colour = "red"
) +
theme_bw() +
labs(y = "Log income", x = "Wave")
There are two general ways to expand this model to include non-linear change. One is by including polynomials while the other is by looking at relative change in time. We will cover both below.
Estimating non-linear LGM using polynomials
Including polynomials to model nonlinear effects has a similar motivation to regression modelling. A polynomial (or interaction) allows the effect to change depending on the values of a predictor. In the case of LGM, this would mean that we allow the slope to be higher or lower as time passes. This, in effect, would bend the trend upwards or downwards. If we want to allow for multiple bends, then we need to include multiple polynomials. Below, we will include just the square effects modelled as a latent variable “q” (but the model can be easily expanded to include cubed effects and so on).
model <- 'i =~ 1*logincome_1 + 1*logincome_2 + 1*logincome_3 +
1*logincome_4 + 1*logincome_5 + 1*logincome_6
s =~ 0*logincome_1 + 1*logincome_2 + 2*logincome_3 +
3*logincome_4 + 4*logincome_5 + 5*logincome_6
q =~ 0*logincome_1 + 1*logincome_2 + 4*logincome_3 +
9*logincome_4 + 16*logincome_5 + 25*logincome_6'
fit2 <- growth(model, data = usw)
summary(fit2, standardized = TRUE)
## lavaan 0.6-19 ended normally after 75 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 15
##
## Number of observations 21178
##
## Model Test User Model:
##
## Test statistic 461.571
## Degrees of freedom 12
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i =~
## logincome_1 1.000 1.173 0.839
## logincome_2 1.000 1.173 0.941
## logincome_3 1.000 1.173 0.996
## logincome_4 1.000 1.173 1.011
## logincome_5 1.000 1.173 1.018
## logincome_6 1.000 1.173 1.064
## s =~
## logincome_1 0.000 0.000 0.000
## logincome_2 1.000 0.416 0.333
## logincome_3 2.000 0.831 0.706
## logincome_4 3.000 1.247 1.074
## logincome_5 4.000 1.663 1.442
## logincome_6 5.000 2.078 1.885
## q =~
## logincome_1 0.000 0.000 0.000
## logincome_2 1.000 0.069 0.055
## logincome_3 4.000 0.275 0.234
## logincome_4 9.000 0.619 0.533
## logincome_5 16.000 1.101 0.955
## logincome_6 25.000 1.720 1.560
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i ~~
## s -0.235 0.010 -24.013 0.000 -0.482 -0.482
## q 0.020 0.002 12.521 0.000 0.246 0.246
## s ~~
## q -0.026 0.001 -21.952 0.000 -0.905 -0.905
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i 6.851 0.009 736.245 0.000 5.839 5.839
## s 0.095 0.005 18.794 0.000 0.229 0.229
## q -0.007 0.001 -8.026 0.000 -0.106 -0.106
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .logincome_1 0.577 0.013 44.837 0.000 0.577 0.295
## .logincome_2 0.484 0.006 76.119 0.000 0.484 0.311
## .logincome_3 0.438 0.006 76.450 0.000 0.438 0.316
## .logincome_4 0.481 0.006 79.039 0.000 0.481 0.357
## .logincome_5 0.533 0.006 82.386 0.000 0.533 0.401
## .logincome_6 0.387 0.010 38.078 0.000 0.387 0.319
## i 1.377 0.020 70.478 0.000 1.000 1.000
## s 0.173 0.007 24.934 0.000 1.000 1.000
## q 0.005 0.000 21.692 0.000 1.000 1.000The quadratic model estimates an initial increase of 0.095 per wave and a negative quadratic term of -0.007. The average increase therefore becomes smaller in later waves. The variance of “q” is 0.005, showing that respondents differ in the non-linear part of their income trajectories. All three latent variances are positive.
Next, we plot the new estimates of change from the new model. We will use a procedure similar to the one above. The main change is to the formula. Now, we need to add a new term, which is time squared (x^2) multiplied by the coefficient for the square effect (pred_lgm2[, 3]). We also added the line from the linear model for comparison.
pred_lgm2 <- predict(fit2)
pred_lgm2_long <- map(
0:5,
function(x) {
pred_lgm2[, 1] + x * pred_lgm2[, 2] + x^2 * pred_lgm2[, 3]
}
) |>
reduce(cbind) |>
as.data.frame() |>
setNames(str_c("Wave ", 1:6)) |>
mutate(id = row_number()) |>
pivot_longer(-id, names_to = "wave", values_to = "pred")
quadratic_means_plot <- pred_lgm2_long |>
group_by(wave) |>
summarise(pred = mean(pred), .groups = "drop")
ggplot(quadratic_means_plot, aes(wave, pred, group = 1)) +
geom_line(linewidth = 1.5, colour = "blue") +
geom_line(
data = linear_means,
aes(wave, pred, group = 1),
linewidth = 1.5,
colour = "red",
alpha = 0.5
) +
theme_bw() +
labs(y = "Log income", x = "Wave")
The blue quadratic trajectory starts below the red linear trajectory, rises slightly above it in the middle waves and finishes below it. The difference is modest, but the fit statistics show that allowing this curvature improves the model.
Non-linear change in time using relative change
The alternative way to model non-linear change is to estimate relative change. This is similar in spirit to including dummy variables in a regression model. The only thing we need to do is to tweak the loadings for the slope latent variable. We will fix only the first and last loading to 0 and 1. The rest of the loadings will not be fixed and will be estimated. Now, the interpretation of the slope will be the total amount of change from the first wave to the last one. The newly estimated loadings will tell us the proportion of change from the start until that point out of the total change observed.
model <- 'i =~ 1*logincome_1 + 1*logincome_2 + 1*logincome_3 +
1*logincome_4 + 1*logincome_5 + 1*logincome_6
s =~ 0*logincome_1 + logincome_2 + logincome_3 +
logincome_4 + logincome_5 + 1*logincome_6'
fit3 <- growth(model, data = usw)
summary(fit3, standardized = TRUE)
## lavaan 0.6-19 ended normally after 75 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 15
##
## Number of observations 21178
##
## Model Test User Model:
##
## Test statistic 1004.879
## Degrees of freedom 12
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i =~
## logincome_1 1.000 1.122 0.792
## logincome_2 1.000 1.122 0.903
## logincome_3 1.000 1.122 0.962
## logincome_4 1.000 1.122 0.976
## logincome_5 1.000 1.122 0.985
## logincome_6 1.000 1.122 0.992
## s =~
## logincome_1 0.000 0.000 0.000
## logincome_2 0.156 0.016 9.732 0.000 0.124 0.100
## logincome_3 0.468 0.011 41.324 0.000 0.373 0.320
## logincome_4 0.681 0.011 60.655 0.000 0.543 0.473
## logincome_5 0.876 0.012 70.421 0.000 0.698 0.613
## logincome_6 1.000 0.797 0.705
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i ~~
## s -0.554 0.017 -32.126 0.000 -0.620 -0.620
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## i 6.888 0.009 766.311 0.000 6.142 6.142
## s 0.258 0.008 30.516 0.000 0.323 0.323
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .logincome_1 0.746 0.013 58.510 0.000 0.746 0.372
## .logincome_2 0.442 0.008 54.014 0.000 0.442 0.286
## .logincome_3 0.480 0.006 83.935 0.000 0.480 0.353
## .logincome_4 0.521 0.006 85.982 0.000 0.521 0.395
## .logincome_5 0.520 0.007 75.861 0.000 0.520 0.402
## .logincome_6 0.492 0.008 63.034 0.000 0.492 0.385
## i 1.258 0.019 65.513 0.000 1.000 1.000
## s 0.635 0.020 31.726 0.000 1.000 1.000From wave 1 to wave 6, the relative-change model estimates a total average increase of 0.258. The loadings are 0, 0.156, 0.468, 0.681, 0.876 and 1. About 16% of the total change occurs by wave 2 and 47% by wave 3. The largest increment, about 31%, occurs between waves 2 and 3, while about 12% occurs between waves 5 and 6. This pattern is not perfectly linear.
We need to extract the loadings to make a nice graph using the formula. We can use the parameterestimates() command to do that:
parameterEstimates(fit3) ## lhs op rhs est se z pvalue ci.lower ci.upper ## 1 i =~ logincome_1 1.000 0.000 NA NA 1.000 1.000 ## 2 i =~ logincome_2 1.000 0.000 NA NA 1.000 1.000 ## 3 i =~ logincome_3 1.000 0.000 NA NA 1.000 1.000 ## 4 i =~ logincome_4 1.000 0.000 NA NA 1.000 1.000 ## 5 i =~ logincome_5 1.000 0.000 NA NA 1.000 1.000 ## 6 i =~ logincome_6 1.000 0.000 NA NA 1.000 1.000 ## 7 s =~ logincome_1 0.000 0.000 NA NA 0.000 0.000 ## 8 s =~ logincome_2 0.156 0.016 9.732 0 0.125 0.187 ## 9 s =~ logincome_3 0.468 0.011 41.324 0 0.446 0.490 ## 10 s =~ logincome_4 0.681 0.011 60.655 0 0.659 0.704 ## 11 s =~ logincome_5 0.876 0.012 70.421 0 0.851 0.900 ## 12 s =~ logincome_6 1.000 0.000 NA NA 1.000 1.000 ## 13 logincome_1 ~~ logincome_1 0.746 0.013 58.510 0 0.721 0.771 ## 14 logincome_2 ~~ logincome_2 0.442 0.008 54.014 0 0.426 0.458 ## 15 logincome_3 ~~ logincome_3 0.480 0.006 83.935 0 0.469 0.491 ## 16 logincome_4 ~~ logincome_4 0.521 0.006 85.982 0 0.509 0.533 ## 17 logincome_5 ~~ logincome_5 0.520 0.007 75.861 0 0.507 0.534 ## 18 logincome_6 ~~ logincome_6 0.492 0.008 63.034 0 0.477 0.508 ## 19 i ~~ i 1.258 0.019 65.513 0 1.220 1.296 ## 20 s ~~ s 0.635 0.020 31.726 0 0.596 0.674 ## 21 i ~~ s -0.554 0.017 -32.126 0 -0.588 -0.520 ## 22 logincome_1 ~1 0.000 0.000 NA NA 0.000 0.000 ## 23 logincome_2 ~1 0.000 0.000 NA NA 0.000 0.000 ## 24 logincome_3 ~1 0.000 0.000 NA NA 0.000 0.000 ## 25 logincome_4 ~1 0.000 0.000 NA NA 0.000 0.000 ## 26 logincome_5 ~1 0.000 0.000 NA NA 0.000 0.000 ## 27 logincome_6 ~1 0.000 0.000 NA NA 0.000 0.000 ## 28 i ~1 6.888 0.009 766.311 0 6.871 6.906 ## 29 s ~1 0.258 0.008 30.516 0 0.241 0.274
With some manipulation, we can extract just what we want:
loadings <- parameterEstimates(fit3) |> filter(lhs == "s", op == "=~") |> pull(est) loadings ## [1] 0.0000000 0.1559671 0.4677954 0.6814987 0.8756568 1.0000000
We can follow a similar approach to the one before to create the long data with predicted scores from the LGM. The only difference is that we now loop over the loadings instead of the numbers 0 to 5:
pred_lgm3 <- predict(fit3)
pred_lgm3_long <- map(
loadings,
function(x) pred_lgm3[, 1] + x * pred_lgm3[, 2]
) |>
reduce(cbind) |>
as.data.frame() |>
setNames(str_c("Wave ", 1:6)) |>
mutate(id = row_number()) |>
pivot_longer(-id, names_to = "wave", values_to = "pred")
relative_means_plot <- pred_lgm3_long |>
group_by(wave) |>
summarise(pred = mean(pred), .groups = "drop")
ggplot(relative_means_plot, aes(wave, pred, group = 1)) +
geom_line(linewidth = 1.5, colour = "green") +
geom_line(
data = quadratic_means_plot,
aes(wave, pred, group = 1),
linewidth = 1.5,
colour = "blue",
alpha = 0.5
) +
geom_line(
data = linear_means,
aes(wave, pred, group = 1),
linewidth = 1.5,
colour = "red",
alpha = 0.5
) +
theme_bw() +
labs(y = "Log income", x = "Wave")
The relative-change trajectory shows the same broad pattern. Income increases most quickly in the earlier waves and the rate of increase becomes smaller towards wave 6.
Want to really learn? Follow along using the data and code.
Get them when you join the newsletter. Additionally, you will get my 63-point pre-submission checklist for your research. All for free!
Finding the best fit
We can use chi-square, AIC and BIC to compare the models. The linear and quadratic models are nested, but the quadratic and relative-change models have the same degrees of freedom and are not nested. Their chi-square difference is therefore descriptive rather than a likelihood-ratio test:
model_comparison <- anova(fit1, fit2, fit3) model_comparison ## ## Chi-Squared Difference Test ## ## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) ## fit2 12 341502 341621 461.57 ## fit3 12 342045 342164 1004.88 543.31 0.000000 0 ## fit1 16 342154 342242 1122.27 117.40 0.036587 4 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The quadratic model has the lowest chi-square, AIC and BIC. It improves on the linear model by 660.70 chi-square points with four additional parameters. The information criteria also favour it over the relative-change model, so we retain the quadratic specification.
Conclusions
Hopefully, that will give you an idea of how to estimate non-linear LGM, how to visualize this change, and how to interpret it. Visualizing these models is always helpful, as the interpretation can get quite tricky.
If you liked this, you could look at other blog posts, such as this introduction to multilevel modelling for longitudinal data or this one visualizing transition in time for categorical variables. You can also learn how to include time-constant and time-varying predictors in LGM models here and here.
Want to take your skills to the next level? Join our next live course to learn how to efficiently prepare and explore data as well as the main frameworks for analysing longitudinal data.