setwd("C:/Work Files/8190 MLM/Data Assignment Keys")Data Assignment 2 Answer Key
After setting the working directory, you need to load the requisite packages. I always load tidyverse which includes ggplot2 among several other data management packages. Because you are running both latent growth models (SEM) and MLM growth models you will need both lavaan and lme4. Because lme4 does not provide significance tests for whatever reason, you also need to load lmerTest.
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lavaan)This is lavaan 0.6-18
lavaan is FREE software! Please report any bugs.library(lme4)Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpacklibrary(lmerTest)
Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':
lmer
The following object is masked from 'package:stats':
stepNext load the data set for the assignment
Assn_2_Data <-read.csv("Assn_2_Data.csv")Once you have loaded the data set, you should look over the data. The head command will give you a quick look at the data in a spreadsheet format. I also usually run the names command to make sure I have a reference for the names of the variables and the columns in which they are located. Finally it is important to run the str command to see how the data is structured in terms of variable types.
head(Assn_2_Data) StudentID Time GPA Work Sex
1 1 1 2.3 2 1
2 1 2 2.1 2 1
3 1 3 3.0 2 1
4 1 4 3.0 2 1
5 1 5 3.0 2 1
6 1 6 3.3 2 1names(Assn_2_Data)[1] "StudentID" "Time" "GPA" "Work" "Sex" str(Assn_2_Data)'data.frame': 1200 obs. of 5 variables:
$ StudentID: int 1 1 1 1 1 1 2 2 2 2 ...
$ Time : int 1 2 3 4 5 6 1 2 3 4 ...
$ GPA : num 2.3 2.1 3 3 3 3.3 2.2 2.5 2.6 2.6 ...
$ Work : int 2 2 2 2 2 2 2 3 2 2 ...
$ Sex : int 1 1 1 1 1 1 0 0 0 0 ...The data looks like it loaded correctly and there are no issues with the variable names (e.g., sometimes the id gets loaded as ~…id when pulling data in from a csv or excel file.) One thing to note is that both Sex and Work are shown as integer data, however, they are actually factors. Since Sex is coded as binary technically you don’t have to change it, but you might not get data labels in your output if you leave it as an integer. Work on the other hand has 3 categories and lme4 will read it as an unordered nominal variable and will not handle it correctly in regression-based models. So before proceeding you will need to convert both Sex and Work to factors.
Assn_2_Data$Work <- factor(Assn_2_Data$Work,
levels=c(1,2,3),
labels = c("No Work","Part-time", "Full-time"))
Assn_2_Data$Sex <- factor(Assn_2_Data$Sex,
levels = c(0,1),
labels = c("Male","Female"))
str(Assn_2_Data$Sex) Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 1 1 1 1 ...str(Assn_2_Data$Work) Factor w/ 3 levels "No Work","Part-time",..: 2 2 2 2 2 2 2 3 2 2 ...table(Assn_2_Data$Work)
No Work Part-time Full-time
51 971 178 table(Assn_2_Data$Sex)
Male Female
408 408 After running the code to convert Sex and Work to factors I usually run the str command on those variables separately **(it wouldn’t make much difference if you did it for the whole data set in this case, because there are a small number of variables, but with a larger data set confirming that the structure was changed appropriately would be more time consuming.) **
I also usually run a simple table for each factor to make sure that the distributions are correct and that the labels are being applied correctly.
The next step is to plot the trajectories at the sample and individual levels. This serves the same purpose as looking at the univariate descriptive statistics and histograms before conducting non-nested analyses; it gives a quick insight into how the outcome is changing over time and how much variability in trajectories there is within the sample. First, create mean values and standard errors for the outcome of interest (GPA). This will make plotting the sample average trajectory more straightforward and it will allow for the inclusion of 95% CI error bars.
mean_data <- Assn_2_Data %>%
group_by(Time) %>%
summarise(mean_GPA = mean(GPA),
se_GPA = sd(GPA) / sqrt(n()))
mean_traj_plot <-ggplot(mean_data, aes(x = Time, y = mean_GPA)) +
geom_line(color = "blue", size = 1) +
geom_point(color = "blue", size = 2) +
geom_errorbar(aes(ymin = mean_GPA - se_GPA, ymax = mean_GPA + se_GPA), width = 0.2) +
labs(title = "Mean Trajectory Over Time",
x = "Time",
y = "GPA") +
theme_minimal()Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.mean_traj_plot
Looking at the plot for changes in GPA overtime, it is clear that there is a general upward trend from grade 11 to 4th year in university (11th grade to college senior in American English) . Also, the error bars look to be very similar so the time specific variance seems to be homoscedastic. Because ggplot is using the range to determine the scaling of the y-axis, this upward trajectory looks steeper than it actually is. To get a more realistic graph, the min and max for the y-axis can be set to range from 0-4 (full range of GPA).
mean_traj_plot +
ylim(0,4) 
The mean GPA trajectory suggests a linear growth process and the time-specific variance looks homoscedastic, but this doesn’t give much insight into the variability in intra-individual variance. To get a better sense of that we need to run spaghetti plots. To run spaghetti plots in ggplot2 the data has to be in long form (or at least it is easier to do in long form). The x-axis will be our index of time, the y-axis is the outcome (GPA) and group is set to StudentID so that ggplot will generate a separate plot for each participant by ID. The trajectories can also be color coded to assess potential group differences in trajectory by predictors of interest. In this example I am using participant sex.
ggplot(Assn_2_Data, aes(x = Time, y = GPA, group = StudentID, color = Sex)) +
geom_line(alpha = 0.5) +
geom_point(alpha = 0.5) +
labs(title = "Individual Trajectories (Spaghetti Plot)",
x = "Time",
y = "Outcome") +
theme_minimal()
In this spaghetti plot the general upward trend in GPA is still clear, but there is a good deal of variability. There are some individual nonlinearities, but the overall pattern is still linear and so a nonlinear growth model is not really warranted. The color coding by participant gender also suggests that there might be gender differences, at least in intercept, with female students having slightly higher GPA’s on average.
Before moving on to testing the model, we need to take a look at the missingness patterns. First run the descriptives to see how much missing data there is. I personally like the psych package for this.
library(psych)
Attaching package: 'psych'The following object is masked from 'package:lavaan':
cor2covThe following objects are masked from 'package:ggplot2':
%+%, alphanames(Assn_2_Data[,-1])[1] "Time" "GPA" "Work" "Sex" describe(Assn_2_Data[,-1]) vars n mean sd median trimmed mad min max range skew kurtosis se
Time 1 1200 3.50 1.71 3.5 3.50 2.22 1.0 6 5.0 0.00 -1.27 0.05
GPA 2 1200 2.87 0.39 2.8 2.86 0.44 1.7 4 2.3 0.18 -0.02 0.01
Work* 3 1200 2.11 0.42 2.0 2.06 0.00 1.0 3 2.0 0.62 1.90 0.01
Sex* 4 816 1.50 0.50 1.5 1.50 0.74 1.0 2 1.0 0.00 -2.00 0.02The data is mostly complete with the exception of participant sex, which is missing quite a bit. To look at this in more dept there are a few options. mice, VIM, nanair, missForest, and Amelia are the most comprehensive packages for visualizing and imputing missing data. For our purposes VIM will do.
library(VIM)Loading required package: colorspaceLoading required package: gridVIM is ready to use.Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
Attaching package: 'VIM'The following object is masked from 'package:datasets':
sleepaggr_plot <- aggr(Assn_2_Data[,-1], col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(Assn_2_Data[,-1]), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))
Variables sorted by number of missings:
Variable Count
Sex 0.32
Time 0.00
GPA 0.00
Work 0.00marginplot(Assn_2_Data[c(3,5)])
The first graph displays the percent of missing data by variable in a histogram format. We can see that sex is missing 32% of the values. You could run a Little’s MCAR test here with the nanair package. Unfortunately it is not available for the newest version of R - there are ways around that, but we don’t need to spend time with it now. Perhaps a better indicator anyway is to examine how the missingness is related to our outcome GPA. The best option for this data would be to only run it as a latent growth curve (LGC) rather than a MLM because LGC can use FIML. The other option is to impute the missing data, but I am not completely comfortable with imputing demographics, especially with a limited number of variables. But, this is an assignment so we will do both.
Modeling Growth using an MLM approach
It isn’t necessary, but is advisable to first conduct an unconditional growth model before adding in the predictors. This is useful in testing model assumptions (e.g., linear v. nonlinear) and assessing the degree of variance in both intercept and slope.
To specify this model we will use the lmer (linear mixed effects regression) command from lme4 with GPA as the DV, the tilda indicates a regression and the level 1 predictor is time. To specify the level 2 or random effects we will add to the simple regression a set of parentheses and a vertical pipe. The left side of the pipe specifies the random effects. A value of 1 specifies a random intercept - if we did not include a 1 or put in a zero the model would only estimate fixed effects for the intercept. We also want random effects for our level 1 predictor (time) which will give us random slopes as well. To the right of the vertical pipe we indicate the nesting variable. Since for growth curves, time is nested within person, our index of person (StudentID) is the nesting variable. Finally, we specify which data set to use in estimating this model.
UC_model <- lmer(GPA ~ Time +
(1 + Time | StudentID), data = Assn_2_Data)
# Summary of the model
summary(UC_model)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: GPA ~ Time + (1 + Time | StudentID)
Data: Assn_2_Data
REML criterion at convergence: 261
Scaled residuals:
Min 1Q Median 3Q Max
-3.2695 -0.5377 -0.0128 0.5326 3.1939
Random effects:
Groups Name Variance Std.Dev. Corr
StudentID (Intercept) 0.052503 0.22913
Time 0.004504 0.06711 -0.38
Residual 0.042388 0.20588
Number of obs: 1200, groups: StudentID, 200
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.493e+00 2.112e-02 1.990e+02 118.02 <2e-16 ***
Time 1.063e-01 5.885e-03 1.990e+02 18.07 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
Time -0.579Looking at the output we see that the fixed effect (true sample average) intercept, or GPA at time 1 is 2.49 and the fixed effect slope is .01. So for every year we add .01 GPA points to the sample average at baseline. The intercept and slope are correlated at -0.38 indicating that for every unit of GPA a student’s baseline GPA is higher than the mean we have to subtract 0.38 from the trajectory over time and vice versa for scores below the mean.
lme4 does not have a convenient way to check the significance of the random effects, but we can check it by comparing the model fit with and without a random effect specified. The random effect for the slope (on average the degree to which individual slopes differ from each other) is .004, but we can’t tell if that is significantly different from 0. So we will specify an new model with only a random intercept by removing the time variable from the parentheses. We need to use a new model name as well. Then we can see if adding the random effect significantly improved model fit.
UC_model_2 <- lmer(GPA ~ Time +
(1 | StudentID), data = Assn_2_Data)
anova(UC_model, UC_model_2)refitting model(s) with ML (instead of REML)Data: Assn_2_Data
Models:
UC_model_2: GPA ~ Time + (1 | StudentID)
UC_model: GPA ~ Time + (1 + Time | StudentID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
UC_model_2 4 401.65 422.01 -196.82 393.65
UC_model 6 258.23 288.77 -123.12 246.23 147.41 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The output indicates that the random variation in slopes is significant. The model AIC and BIC are both substantially lower for the random slope effects model and the chi-square difference test is significant. The key thing here is that because the intercepts and slopes have significant variation between participants we can try and model that variability with person level (level 2) predictors. In the conditional model below, participant sex, participant work level, and the interaction of sex by work are included as predictors of the initial GPA and the rate of change in GPA over time.
Cond_model <- lmer(GPA ~ Time + Sex + Work + Time:Sex + Time:Work +
(1 + Time | StudentID), data = Assn_2_Data)
# Summary of the model
summary(Cond_model)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: GPA ~ Time + Sex + Work + Time:Sex + Time:Work + (1 + Time |
StudentID)
Data: Assn_2_Data
REML criterion at convergence: 95.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.2310 -0.5353 0.0183 0.5296 2.8413
Random effects:
Groups Name Variance Std.Dev. Corr
StudentID (Intercept) 0.028691 0.1694
Time 0.002992 0.0547 -0.64
Residual 0.044290 0.2105
Number of obs: 816, groups: StudentID, 136
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.917670 0.247577 663.771899 11.785 <2e-16 ***
Time -0.008471 0.054129 690.942440 -0.157 0.8757
SexFemale 0.107903 0.044735 136.851572 2.412 0.0172 *
WorkPart-time -0.515520 0.246287 649.054112 -2.093 0.0367 *
WorkFull-time -0.634121 0.250559 657.980488 -2.531 0.0116 *
Time:SexFemale 0.027976 0.012860 136.504978 2.175 0.0313 *
Time:WorkPart-time 0.081672 0.053605 674.843050 1.524 0.1281
Time:WorkFull-time 0.086795 0.054932 686.010531 1.580 0.1146
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Time SexFml WrkPr- WrkFl- Tm:SxF Tm:WP-
Time -0.958
SexFemale -0.124 0.129
WorkPart-tm -0.991 0.949 0.030
WorkFull-tm -0.978 0.938 0.052 0.978
Time:SexFml 0.097 -0.158 -0.772 -0.024 -0.045
Tm:WrkPrt-t 0.954 -0.984 -0.033 -0.963 -0.941 0.034
Tm:WrkFll-t 0.935 -0.969 -0.059 -0.934 -0.960 0.063 0.969Here we see that the intercept, which represents the baseline GPA when our predictors (sex and work) are at 0, so this represents the average GPA for males (sex = 0) and not working (work = 0) = 2.92. The average rate of change in GPA over time is no longer significant. The fixed effect for sex is significant meaning that on average female students had a GPA .11 units higher or (2.92 + 0.11) 3.03 at baseline. Working part-time (b=-0.51) and full-time (b=-0.63) were both significant predictors of the intercept relative to not working. Again we start with a baseline sample average GPA of 2.92. For students working part-time that initial GPA drops by .51 GPA units or a baseline average GPA of 2.41. Working full-time reduces the baseline even further (b=-0.63) for a baseline average GPA of 2.29.
Now looking at the slope we start with the average rate of change at -0.008, but for female students the average rate of change increases by (b=0.028) relative to males. In other words the average rate of change in GPA each year is (-0.008 + 0.028) = .02. We can then say that for female students they have an average GPA of 3.03 and it increases by .02 units every year. Neither part-time work or full-time work had a significant impact on the rate of change in GPA over time.
Modeling growth using an SEM framework
Next we can assess growth using a latent growth model or SEM approach. Before we can specify the growth model we need to rotate the data from long format to wide format. To save code and time I did this in SPSS and am just importing that data file directly from SPSS using the haven package.
library(haven)
Assn_2_Data_Wide <-read_sav("Assn_2_Data_Wide.sav")Just like with the MLM, it is a good idea to start with an unconditional growth model. This is a highly constrained SEM with two latent variables, a latent intercept and a latent slope. The factor loadings for the intercept are all fixed to 1, because we need to borrow the scaling metric from the DV (GPA). These factor loadings cannot be changed because doing so would conflate real mean level changes overtime with changes in measurement scale over time. The latent slope factor loadings are where we include our index of time in the model. The factor loadings correspond to the number of time units from baseline for each wave. At baseline the factor loading is 0, because it is 0 time units from baseline, the second factor loading is 1 because it is one time unit (year) from baseline, the third factor loading is 2 because it is two time units from the baseline, etc. The model fit commands are the same as those used in SEM. Importantly, we can use FIML to handle missing data. The unconditional models will both have n=200 because sex is not included in either model and that is the only variable with missing data. However, there will be a sample size difference for the conditional models.
UC_LGC <-'
I=~ 1*GPA.1 + 1*GPA.2 + 1*GPA.3 + 1*GPA.4 + 1*GPA.5 + 1*GPA.6
S=~ 0*GPA.1 + 1*GPA.2 + 2*GPA.3 + 3*GPA.4 + 4*GPA.5 + 5*GPA.6
'
UC_LGC_fit <- growth(UC_LGC, estimator= "ML", data=Assn_2_Data_Wide, mimic = "Mplus", missing = "FIML")
summary(UC_LGC_fit, fit.measures = TRUE, standardized=TRUE)lavaan 0.6-18 ended normally after 73 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 11
Number of observations 200
Number of missing patterns 1
Model Test User Model:
Test statistic 43.945
Degrees of freedom 16
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 811.632
Degrees of freedom 15
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.965
Tucker-Lewis Index (TLI) 0.967
Robust Comparative Fit Index (CFI) 0.965
Robust Tucker-Lewis Index (TLI) 0.967
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -49.385
Loglikelihood unrestricted model (H1) -27.413
Akaike (AIC) 120.770
Bayesian (BIC) 157.052
Sample-size adjusted Bayesian (SABIC) 122.203
Root Mean Square Error of Approximation:
RMSEA 0.093
90 Percent confidence interval - lower 0.061
90 Percent confidence interval - upper 0.127
P-value H_0: RMSEA <= 0.050 0.016
P-value H_0: RMSEA >= 0.080 0.772
Robust RMSEA 0.093
90 Percent confidence interval - lower 0.061
90 Percent confidence interval - upper 0.127
P-value H_0: Robust RMSEA <= 0.050 0.016
P-value H_0: Robust RMSEA >= 0.080 0.772
Standardized Root Mean Square Residual:
SRMR 0.233
Parameter Estimates:
Standard errors Standard
Information Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
I =~
GPA.1 1.000 0.187 0.552
GPA.2 1.000 0.187 0.553
GPA.3 1.000 0.187 0.558
GPA.4 1.000 0.187 0.568
GPA.5 1.000 0.187 0.535
GPA.6 1.000 0.187 0.471
S =~
GPA.1 0.000 0.000 0.000
GPA.2 1.000 0.057 0.170
GPA.3 2.000 0.114 0.342
GPA.4 3.000 0.172 0.522
GPA.5 4.000 0.229 0.656
GPA.6 5.000 0.286 0.721
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
I ~~
S 0.002 0.002 1.565 0.118 0.232 0.232
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
I 2.598 0.018 141.886 0.000 13.923 13.923
S 0.106 0.005 20.317 0.000 1.862 1.862
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.GPA.1 0.080 0.010 8.049 0.000 0.080 0.696
.GPA.2 0.071 0.008 8.518 0.000 0.071 0.621
.GPA.3 0.054 0.006 9.020 0.000 0.054 0.483
.GPA.4 0.029 0.003 8.486 0.000 0.029 0.268
.GPA.5 0.015 0.003 5.589 0.000 0.015 0.121
.GPA.6 0.016 0.004 4.336 0.000 0.016 0.100
I 0.035 0.007 4.937 0.000 1.000 1.000
S 0.003 0.001 5.593 0.000 1.000 1.000The results of this model indicate that a linear growth model provides an adequate fit to the data [ \(\chi^2\) (16) = 43.95, p<.001, CFI = .96, TLI = .96, RMSEA = .09, SRMR = .233]. It is notable that the \(\chi^2\) , CFI, and TLI do not align with the RMSEA and SRMR, particularly the latter - one common reason for this is that the SRMR and to a lesser extent the RMSEA work well for estimating fit of the model to the variance/covariance matrix, but do not take into account the mean structure, which often causes them to produce erroneous misfit information (see also, Full article: The Problem with Having Two Watches: Assessment of Fit When RMSEA and CFI Disagree)
Similar to the MLM model we see that the sample average GPA at baseline is 2.59 and the sample average rate of change is 0.11 GPA units per year. These values are the same as the MLM unconditional model
What we are really looking for in this model is the fact that, just like the unconditional MLM, there is significant variability in the intercept (\(\sigma^2\) = .035) and slope ( \[\sigma^2\] = .003) factors. It makes logical sense then to add predictors to see if they account for that variability. We can specify this conditional model by adding code to regress sex and work onto the Intercept and Slope variables.
Cond_LGC <-'
I=~ 1*GPA.1 + 1*GPA.2 + 1*GPA.3 + 1*GPA.4 + 1*GPA.5 + 1*GPA.6
S=~ 0*GPA.1 + 1*GPA.2 + 2*GPA.3 + 3*GPA.4 + 4*GPA.5 + 5*GPA.6
I ~ Sex + Work
S ~ Sex + Work
'
Cond_LGC_fit <- growth(Cond_LGC, estimator= "ML", data=Assn_2_Data_Wide, mimic = "Mplus", missing = "FIML")
summary(Cond_LGC_fit, fit.measures = TRUE, standardized=TRUE)lavaan 0.6-18 ended normally after 79 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 15
Number of observations 200
Number of missing patterns 1
Model Test User Model:
Test statistic 64.422
Degrees of freedom 24
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 993.679
Degrees of freedom 27
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.958
Tucker-Lewis Index (TLI) 0.953
Robust Comparative Fit Index (CFI) 0.958
Robust Tucker-Lewis Index (TLI) 0.953
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) 31.400
Loglikelihood unrestricted model (H1) NA
Akaike (AIC) -32.800
Bayesian (BIC) 16.675
Sample-size adjusted Bayesian (SABIC) -30.846
Root Mean Square Error of Approximation:
RMSEA 0.092
90 Percent confidence interval - lower 0.065
90 Percent confidence interval - upper 0.119
P-value H_0: RMSEA <= 0.050 0.007
P-value H_0: RMSEA >= 0.080 0.781
Robust RMSEA 0.092
90 Percent confidence interval - lower 0.065
90 Percent confidence interval - upper 0.119
P-value H_0: Robust RMSEA <= 0.050 0.007
P-value H_0: Robust RMSEA >= 0.080 0.781
Standardized Root Mean Square Residual:
SRMR 0.160
Parameter Estimates:
Standard errors Standard
Information Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
I =~
GPA.1 1.000 0.189 0.576
GPA.2 1.000 0.189 0.556
GPA.3 1.000 0.189 0.565
GPA.4 1.000 0.189 0.575
GPA.5 1.000 0.189 0.541
GPA.6 1.000 0.189 0.483
S =~
GPA.1 0.000 0.000 0.000
GPA.2 1.000 0.059 0.172
GPA.3 2.000 0.117 0.350
GPA.4 3.000 0.176 0.535
GPA.5 4.000 0.234 0.670
GPA.6 5.000 0.293 0.749
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
I ~
Sex 0.148 0.019 7.626 0.000 0.781 0.634
Work -0.144 0.048 -3.020 0.003 -0.760 -0.251
S ~
Sex 0.037 0.006 6.215 0.000 0.627 0.509
Work 0.028 0.015 1.897 0.058 0.470 0.155
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I ~~
.S -0.001 0.001 -0.921 0.357 -0.193 -0.193
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I 2.610 0.113 23.158 0.000 13.806 13.806
.S -0.024 0.034 -0.711 0.477 -0.417 -0.417
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.GPA.1 0.072 0.009 7.961 0.000 0.072 0.669
.GPA.2 0.073 0.008 8.645 0.000 0.073 0.629
.GPA.3 0.055 0.006 9.097 0.000 0.055 0.492
.GPA.4 0.030 0.004 8.599 0.000 0.030 0.279
.GPA.5 0.017 0.003 6.110 0.000 0.017 0.136
.GPA.6 0.013 0.003 3.734 0.000 0.013 0.083
.I 0.018 0.005 3.374 0.001 0.502 0.502
.S 0.003 0.000 5.045 0.000 0.733 0.733The overall model fit is adequate [ \(\chi^2\) (24) = 64.42, p <.001; CFI = .96, TLI = .95, RMSEA = .09, SRMR = .16] with the caveat that the RMSEA and SRMR are inconsistent with the rest of the fit statistics. For this model the sample average intercept for males (sex = 0) and non-working students (work = 0) is 2.61 and the average rate of change for non-working male students is -0.024. These are a bit different than the MLM model, however, the MLM estimates were based on a sample size of 136 and the LGC model estimates are based on a sample size of 200. But they are still in the same ballpark. Female students have a significantly higher baseline GPA than male students (b = .148) with an average GPA of (2.61 + .148) = 2.76. Similarly work status is also significant for the intercept, reducing the baseline GPA by -0.144 for part-time workers (GPA = 2.46) and -.288 for full-time workers (GPA = 2.32). This is based on the work variable being coded 0,1,2 - to be more explicit you could dummy code the work status variable. Being female increases the rate of change in GPA overtime significantly (b=.037) so the female rate of change would be -0.024 + .037 = .013 GPA units per year. Work status does not influence the rate of change in GPA over time, just like in the MLM version.