Presenting results effectively is a crucial step in longitudinal data analysis. Traditionally this was done separately from the analysis by copying and pasting outputs. In this section, we will explore several tools and techniques to present results from longitudinal data analysis. These tools can be used together with dynamic documents to create a fully reproducible workflow.

Presenting results from regression based models

To show how this can work, let’s begin by fitting several linear mixed models using the lme4 package. These models examine how the mental component score (“sf12mcs”) changes over time and how it varies by sex:

library(lme4)

m0 <- lmer(data = usl, sf12mcs ~ 1 + (1 | pidp))

m1 <- lmer(data = usl, sf12mcs ~ 1 + wave0 + (1 | pidp)) 

m3 <- lmer(data = usl, sf12mcs ~ 1 + wave0 + gndr.fct + 
             (1 + wave0 | pidp)) 

m4 <- lmer(data = usl, sf12mcs ~ 1 + wave0 + gndr.fct + gndr.fct:wave0 + 
             (1 + wave0 | pidp))

These models incrementally build complexity, starting with a random intercept model (m0) and progressing to include time (“wave0”), gender (“gndr.fct”), and their interaction (m4).

The stargazer package is an excellent tool for creating publication-quality tables summarizing regression results. By default, stargazer produces concise summaries of model coefficients and key statistics:

library(stargazer)


stargazer(list(m0, m1, m3, m4), 
          type = "text", 
          title = "Regression Results")

## 
## Regression Results
## ========================================================================
##                                      Dependent variable:                
##                      ---------------------------------------------------
##                                            sf12mcs                      
##                          (1)          (2)          (3)          (4)     
## ------------------------------------------------------------------------
## wave0                              -0.419***    -0.419***    -0.551***  
##                                     (0.020)      (0.022)      (0.033)   
##                                                                         
## gndr.fctFemale                                  -1.373***    -1.635***  
##                                                  (0.071)      (0.086)   
##                                                                         
## wave0:gndr.fctFemale                                          0.236***  
##                                                               (0.044)   
##                                                                         
## Constant              49.813***    50.285***    51.054***    51.201***  
##                        (0.035)      (0.042)      (0.059)      (0.065)   
##                                                                         
## ------------------------------------------------------------------------
## Observations           127,989      127,989      127,989      127,989   
## Log Likelihood       -468,160.100 -467,954.300 -467,570.900 -467,558.600
## Akaike Inf. Crit.    936,326.300  935,916.700  935,155.800  935,133.300 
## Bayesian Inf. Crit.  936,355.600  935,955.700  935,224.200  935,211.400 
## ========================================================================
## Note:                                        *p<0.1; **p<0.05; ***p<0.01

This command generates a basic summary of the models, including fixed effects, random effects, and standard errors. The option type can be adjusted for pdf (type = "latex") or HTML (type = "html") and can be combined with dynamic documents, such as RMarkdown. While informative, the default table format can be further improved, with most of the table being editable.

For example, to make the table more accessible, we can customize variable names, add column labels, and format output:

stargazer(
  m0, m1, m3, m4,
  type = "html",
  dep.var.labels = "Mental Component Score (MCS)",
  column.labels = c("Intercept Only", 
                    "Wave Effect", 
                    "Wave + Gender", 
                    "Wave x Gender Interaction"),
  covariate.labels = c("(Intercept)", "Wave", "Female", "Wave x Female"),
  omit.stat = c("LL", "ser", "f"),
  digits = 3,
  title = "Linear Mixed Models for MCS",
  align = TRUE
)

An alternative package to this is texreg. This package can also be used to create tables from regression models. For example, if we want to present the results from these models we can use the following code:

library(texreg)

htmlreg(list(m0, m1, m3, m4), 
       custom.model.names = c("Intercept Only", 
                              "Wave Effect", 
                              "Wave + Gender", 
                              "Wave x Gender Interaction"),
       custom.coef.names = c("(Intercept)", 
                             "Wave", 
                             "Female", 
                             "Wave x Female"))

Extracting results from regression based models

An alternative strategy to using such models is to extract the coefficients of interest from the models and create our own tables. A useful package for extracting coefficients from regression models are broom and broom.mixed (for mixture models).

library(broom.mixed)
library(knitr)

model_results <- tidy(m3, effects = "fixed")

Once we extract the information we can present it using a simple command such as kable().

kable(model_results, caption = "Fixed Effects Results")

This could be further improved by cleaning the table, adding a caption and restricting the decimal points:

model_results |> 
  select(-effect) |> 
  setNames(c("Predictor", "Estimate", "SE", "t_value", "p_value")) |>
  kable(caption = "Fixed Effects Results", digits = 2)

There are multiple packages that can be used to create tables from regression models, such as gtsummary, flextable, and kableExtra. Each package offers different features and customization options, so it’s worth exploring to find the one that best suits your needs.

Presenting SEM-Based Models

We also discussed using Structural Equation Modelling framework to analyse longitudinal data. Presenting SEM results tends to be more complicated as models can be complex and often researchers need to make decisions regarding what to present. Most of the packages mentioned above do not work in a straightforward way with SEM models. As a result, users often need to extract the model coefficients.

Let’s look at an example. Imagine we run a Latent Growth Model and we want to present the results.

library(lavaan)

model <- 'i =~ 1*logincome_1 + 1*logincome_2 + 
                1*logincome_3 + 1*logincome_4
          s =~ 0*logincome_1 + 1*logincome_2 + 
                2*logincome_3 + 3*logincome_4'

fit1 <- growth(model, data = usw)

Extracting fit indices

Typically we want to extract both the goodness of fit indices as well as the key coefficients in the model. We can extract fit indices using the fitmeasures() command. Running the command leads to a large number of fit indices:

fitMeasures(fit1)

##                  npar                  fmin                 chisq 
##                 9.000                 0.021              1102.403 
##                    df                pvalue        baseline.chisq 
##                 5.000                 0.000             37441.623 
##           baseline.df       baseline.pvalue                   cfi 
##                 6.000                 0.000                 0.971 
##                   tli                  nnfi                   rfi 
##                 0.965                 0.965                 0.965 
##                   nfi                  pnfi                   ifi 
##                 0.971                 0.809                 0.971 
##                   rni                  logl     unrestricted.logl 
##                 0.971           -159108.351           -158557.150 
##                   aic                   bic                ntotal 
##            318234.702            318308.412             26635.000 
##                  bic2                 rmsea        rmsea.ci.lower 
##            318279.810                 0.091                 0.086 
##        rmsea.ci.upper        rmsea.ci.level          rmsea.pvalue 
##                 0.095                 0.900                 0.000 
##        rmsea.close.h0 rmsea.notclose.pvalue     rmsea.notclose.h0 
##                 0.050                 1.000                 0.080 
##                   rmr            rmr_nomean                  srmr 
##                 0.086                 0.100                 0.054 
##          srmr_bentler   srmr_bentler_nomean                  crmr 
##                 0.054                 0.062                 0.034 
##           crmr_nomean            srmr_mplus     srmr_mplus_nomean 
##                 0.041                 0.110                 0.052 
##                 cn_05                 cn_01                   gfi 
##               268.473               365.497                 0.999 
##                  agfi                  pgfi                   mfi 
##                 0.998                 0.357                 0.980 
##                  ecvi 
##                 0.042

We can extract a subset of these indices to create a more concise summary table:

fit_select <- c("chisq", "df", "pvalue", "cfi", "rmsea", "aic", "bic")
fitmeasures(fit1, fit_select)

##      chisq         df     pvalue        cfi      rmsea        aic        bic 
##   1102.403      5.000      0.000      0.971      0.091 318234.702 318308.412

We can generalize this process to multiple models by using the map_df() command which loops our function and gathers the fit measures into a clean table:

model2 <- 'i =~ 1*logincome_1 + 1*logincome_2 + 
                  1*logincome_3 + 1*logincome_4
           s =~ 0*logincome_1 + 1*logincome_2 + 
                  2*logincome_3 + 3*logincome_4
                  
           logincome_1 ~~ a*logincome_1
           logincome_2 ~~ a*logincome_2
           logincome_3 ~~ a*logincome_3
           logincome_4 ~~ a*logincome_4'

fit2 <- growth(model2, data = usw)

gof <- map_df(list(fit1, fit2), fitmeasures, fit_select) |> 
  mutate(model = c("Model 1", "Model 2")) |> 
  select(model, everything())

gof

## # A tibble: 2 × 8
##   model   chisq      df         pvalue    cfi       rmsea      aic      bic     
##   <chr>   <lvn.vctr> <lvn.vctr> <lvn.vct> <lvn.vct> <lvn.vctr> <lvn.vc> <lvn.vc>
## 1 Model 1 1102.403   5          0         0.9706856 0.09077613 318234.7 318308.4
## 2 Model 2 1420.270   8          0         0.9622747 0.08141184 318546.6 318595.7

This table could then be presented in a dynamic report using a command such as kable()

gof |> 
  setNames(str_to_title(names(gof))) |>
  kable(caption = "Fit Measures for Latent Growth Models",
        digits = 3)

Alternatively, we can use a package such as tidySEM to create a table with fit indices. This also works with multiple models:

library(tidySEM)
table_fit(list(fit1, fit2)) |> 
  kable(digits = 2)

Extracting model coefficients

The other type of information we often want to extract are model coefficients, such as regression slopes and intercepts. When we have a complex model such as the one below we need to decide what are the key coefficients to present.

model3 <- 'i =~ 1*logincome_1 + 1*logincome_2 + 
                  1*logincome_3 + 1*logincome_4
           s =~ 0*logincome_1 + 1*logincome_2 +     
                  2*logincome_3 + 3*logincome_4
                  
           i ~ gndr.fct + age
           s ~ gndr.fct + age'

fit3 <- growth(model3, data = usw)

Coefficients can be extracted using the parameterEstimates() function.

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            s =~ logincome_1   0.000 0.000      NA     NA    0.000    0.000
## 6            s =~ logincome_2   1.000 0.000      NA     NA    1.000    1.000
## 7            s =~ logincome_3   2.000 0.000      NA     NA    2.000    2.000
## 8            s =~ logincome_4   3.000 0.000      NA     NA    3.000    3.000
## 9            i  ~    gndr.fct  -0.380 0.017 -22.944  0.000   -0.412   -0.347
## 10           i  ~         age   0.010 0.000  20.761  0.000    0.009    0.011
## 11           s  ~    gndr.fct   0.001 0.006   0.166  0.868   -0.010    0.012
## 12           s  ~         age  -0.003 0.000 -18.949  0.000   -0.003   -0.003
## 13 logincome_1 ~~ logincome_1   0.505 0.011  45.998  0.000    0.483    0.526
## 14 logincome_2 ~~ logincome_2   0.732 0.008  88.993  0.000    0.716    0.748
## 15 logincome_3 ~~ logincome_3   0.703 0.008  92.653  0.000    0.688    0.718
## 16 logincome_4 ~~ logincome_4   0.573 0.010  58.180  0.000    0.554    0.593
## 17           i ~~           i   1.412 0.017  84.015  0.000    1.379    1.445
## 18           s ~~           s   0.090 0.002  38.548  0.000    0.085    0.094
## 19           i ~~           s  -0.245 0.005 -46.828  0.000   -0.255   -0.235
## 20    gndr.fct ~~    gndr.fct   0.246 0.000      NA     NA    0.246    0.246
## 21    gndr.fct ~~         age  -0.172 0.000      NA     NA   -0.172   -0.172
## 22         age ~~         age 297.013 0.000      NA     NA  297.013  297.013
## 23 logincome_1 ~1               0.000 0.000      NA     NA    0.000    0.000
## 24 logincome_2 ~1               0.000 0.000      NA     NA    0.000    0.000
## 25 logincome_3 ~1               0.000 0.000      NA     NA    0.000    0.000
## 26 logincome_4 ~1               0.000 0.000      NA     NA    0.000    0.000
## 27    gndr.fct ~1               1.559 0.000      NA     NA    1.559    1.559
## 28         age ~1              47.520 0.000      NA     NA   47.520   47.520
## 29           i ~1               6.916 0.036 194.030  0.000    6.846    6.986
## 30           s ~1               0.230 0.012  19.332  0.000    0.207    0.253

The variables “lhs”, “op” and “rhs” can be used to identify the type of coefficient presented. For example, values of “=~” indicate a regression coefficient, while values of “~1” indicate an intercept and “~” a regression coefficient. We can filter the table to extract only the coefficients of interest and present them in a clean table. We use the spread() command to create a wide table and the setNames() command to rename the columns.

parameterEstimates(fit3) |> 
  filter(op == "~" | op == "~1") |> 
  filter(lhs %in% c("i", "s")) |> 
  mutate(rhs = str_to_title(rhs), 
         rhs = ifelse(op == "~1", "Intercept", rhs),
         coef = str_c(round(est, 3), " (", round(se, 3), ")")) |> 
  select(lhs, rhs, coef) |> 
  spread(lhs, coef) |> 
  setNames(c("Predictor", "Intercept", "Slope")) |> 
  kable()

We can also use the table_results() function from the tidySEM package to create a table with the coefficients:

table_results(fit3) |> 
  kable()

And again, with some cleaning it be presented in a nicer way

table_results(fit3) |> 
  filter(str_detect(label, "Means.i|Means.s|ON")) |>
  kable()

Visualizing SEM models

An alternative way to present results from SEM models is to use visualizations. Packages such as, tidySEM, semPlot and lavaanPlot can be used to create diagrams that show the model structure and estimated coefficients. For example, the graph_sem() function from the tidySEM package can be used to create a path diagram of the model:

graph_sem(fit3)

This could be further customized as discussed here.

Another option is to use the semPlot package to create diagrams that show model structure and estimated coefficients:

library(semPlot)
semPaths(fit3)

We can further enhance the graph by adding labels and changing the layout:

semPlot::semPaths(
  fit3,
  layout = "tree",
  sizeMan = 8,
  residuals = FALSE,
  rotation = 1,
  whatLabels = "est",
  nCharNodes = 5,
  normalize = TRUE,
  width = 12,
  height = 6
)

You can also explore other packages, such as lavaanPlot to see which one fits best what you want to present.

Conclusions

Presenting results from longitudinal data analysis is a critical step in the research process. By using tools like stargazer, texreg, broom, and tidySEM, you can create clear, concise tables that summarize your findings. These tables can be further customized to include variable names, column labels, and formatting options. Additionally, packages like semPlot and lavaanPlot can help you create visualizations that illustrate the structure of your models and the relationships between variables.

By combining these approaches with dynamic documents (as discussed here) and integrating them in reproducible workflows (as presented here) you can create reports that are transparent, reproducible, and easy to share with others.