Covariates, outcomes and staged estimation
Written by Mohammed Saqr and Sonsoles López-Pernas
Source:vignettes/articles/covariates.Rmd
covariates.RmdVariables outside the measurement model can be related to the latent classes in two ways. The three-step approach fits the profiles first and then relates them to an external variable, correcting for classification error, so the profiles stay fixed. The one-step approach adds the variable to the model, so it can also change the profiles. Staged estimation applies the same idea of fixed profiles to the group classes.
Data
course_engagement has one row per course enrolment:
1,422 enrolments from 106 students. student identifies the
group, and five activity measures (browse,
lectures, forum_read, forum_post,
attendance) are the indicators. previous_grade
is the standardized grade from the course before each enrolment; it is
not an indicator, so it can serve as an external variable.
The measurement model
To fit the profiles without the external variable, we call
multilpa() with vars for the indicators,
id for the group, n_profiles,
n_group_classes, n_starts and
seed.
fit <- multilpa(course_engagement,
vars = c("browse", "lectures", "forum_read", "forum_post",
"attendance"),
id = "student", n_profiles = 2, n_group_classes = 2,
n_starts = 3, seed = 1)
fit
#> Two-level latent profile analysis: 2 profiles, 2 group classes
#> 1422 individuals in 106 groups; varying diagonal residual covariance (VVI)
#> Log likelihood: -8439.816424 | AIC: 16925.633 | BIC (groups): 16986.892
#> Converged: TRUE | iterations: 9 | best start: 1/3
#>
#> profile browse lectures forum_read forum_post attendance count proportion
#> 1 0.540 0.427 0.620 0.530 0.649 834 0.587
#> 2 -0.766 -0.606 -0.879 -0.753 -0.921 588 0.413
#>
#> Variances and standard errors: get_results(x, "profiles").
#> Every other table: get_results(x, what = ), or get_results(x, "all").Classification error
To obtain the classification-error matrix, we call
get_results() with "classification_errors" and
level for the observations or the groups. Each row gives
the probability of each assigned class for cases in a true class; the
three-step methods use it to correct for misclassification.
get_results(fit, "classification_errors", level = "individuals")
#> level true_class assigned_class probability
#> 1 individuals 1 1 0.9872
#> 2 individuals 1 2 0.0128
#> 3 individuals 2 1 0.0214
#> 4 individuals 2 2 0.9786Class means of an external variable
To estimate the mean of an external variable in each profile, we call
three_step() with data for the data frame and
outcome for the variable. The default method, BCH, corrects
for classification error; standard errors are clustered on groups.
three_step(fit, data = course_engagement, outcome = "previous_grade")
#> Distal outcome by class (bch, 95% CI)
#>
#> Class Mean 95% CI SE Effective n
#> ----- ----- -------------- ----- -----------
#> 1 0.22 [ 0.15, 0.29] 0.035 810.27
#> 2 -0.31 [-0.39, -0.24] 0.036 564.36To compare the profiles, we call three_step() with
contrast = "pairs", which returns each difference between
profile means with a Benjamini-Hochberg adjusted p-value.
three_step(fit, data = course_engagement, outcome = "previous_grade",
contrast = "pairs")
#> Distal outcome: class differences (bch, 95% CI)
#>
#> Class vs Difference 95% CI p p (BH)
#> ----- -- ---------- -------------- ----- ------
#> 2 1 -0.54 [-0.64, -0.43] <.001 <.001To see the effect of ignoring classification error, we call
three_step() with method = "modal", which
treats each case’s most probable profile as known.
three_step(fit, data = course_engagement, outcome = "previous_grade",
method = "modal")
#> Distal outcome by class (modal, 95% CI)
#>
#> Class Mean 95% CI SE Effective n
#> ----- ----- -------------- ----- -----------
#> 1 0.21 [ 0.15, 0.28] 0.034 836
#> 2 -0.30 [-0.37, -0.23] 0.036 586Predicting class membership
To regress profile membership on a covariate, we call
r3step() with data for the data frame and
covariates for the predictors. It returns multinomial logit
coefficients, corrected for classification error.
vcov_type = "robust" clusters the standard errors on
groups.
r3step(fit, data = course_engagement, covariates = "previous_grade",
vcov_type = "robust")
#> level outcome term estimate standard_error statistic p_value p_adjusted
#> 1 individuals profile_1 (Intercept) 0.376 0.1303 2.89 3.90e-03 NA
#> 2 individuals profile_1 previous_grade 0.575 0.0593 9.69 3.35e-22 3.35e-22
#> conf_low conf_high
#> 1 0.121 0.631
#> 2 0.459 0.691One-step estimation
To estimate the covariate together with the profiles, we call
multilpa() with profile_covariates for
predictors of profile membership. group_covariates does the
same for group-class membership and needs a variable that is constant
within groups. Covariate models require complete data.
one_step <- multilpa(course_engagement,
vars = c("browse", "lectures", "forum_read",
"forum_post", "attendance"),
id = "student", n_profiles = 2, n_group_classes = 2,
profile_covariates = "previous_grade", n_starts = 3,
seed = 1)
one_step
#> Multilevel LPA with covariates: 2 profiles, 2 group classes
#> Log likelihood -8422.837649; AIC 16893.675; BIC (groups) 16957.598; converged TRUE
#>
#> profile browse lectures forum_read forum_post attendance count proportion
#> 1 0.540 0.429 0.619 0.530 0.648 834 0.587
#> 2 -0.766 -0.609 -0.879 -0.753 -0.920 588 0.413
#>
#> Variances and standard errors: get_results(x, "profiles").
#> Every other table: get_results(x, what = ), or get_results(x, "all").To obtain the membership coefficients, we call
get_results() with "coefficients", and for
their standard errors and intervals,
parameter_inference().
get_results(one_step, "coefficients")
#> level outcome term parameter estimate
#> 1 profile profile_1 group_class_1 logit -1.539
#> 2 profile profile_1 group_class_2 logit 1.267
#> 3 profile profile_1 previous_grade logit 0.422
#> 4 group group_class_1 (Intercept) logit -0.763
parameter_inference(one_step)
#> level outcome term parameter estimate standard_error statistic p_value
#> 1 measurement profile_1 browse mean 0.540 0.0271 19.88 6.21e-88
#> 2 measurement profile_1 lectures mean 0.429 0.0324 13.25 4.74e-40
#> 3 measurement profile_1 forum_read mean 0.619 0.0247 25.02 3.74e-138
#> 4 measurement profile_1 forum_post mean 0.530 0.0294 18.02 1.38e-72
#> 5 measurement profile_1 attendance mean 0.648 0.0223 29.12 1.87e-186
#> 6 measurement profile_2 browse mean -0.766 0.0311 -24.65 4.13e-134
#> 7 measurement profile_2 lectures mean -0.609 0.0309 -19.71 1.68e-86
#> 8 measurement profile_2 forum_read mean -0.879 0.0261 -33.65 3.58e-248
#> 9 measurement profile_2 forum_post mean -0.753 0.0277 -27.14 3.78e-162
#> 10 measurement profile_2 attendance mean -0.920 0.0264 -34.87 2.34e-266
#> 11 measurement profile_1 browse variance 0.591 0.0296 NA NA
#> 12 measurement profile_1 lectures variance 0.842 0.0419 NA NA
#> 13 measurement profile_1 forum_read variance 0.483 0.0245 NA NA
#> 14 measurement profile_1 forum_post variance 0.688 0.0347 NA NA
#> 15 measurement profile_1 attendance variance 0.385 0.0196 NA NA
#> 16 measurement profile_2 browse variance 0.526 0.0321 NA NA
#> 17 measurement profile_2 lectures variance 0.538 0.0320 NA NA
#> 18 measurement profile_2 forum_read variance 0.363 0.0227 NA NA
#> 19 measurement profile_2 forum_post variance 0.422 0.0262 NA NA
#> 20 measurement profile_2 attendance variance 0.375 0.0228 NA NA
#> 21 profile profile_1 group_class_1 coefficient -1.539 0.1677 -9.18 4.35e-20
#> 22 profile profile_1 group_class_2 coefficient 1.267 0.0943 13.44 3.62e-41
#> 23 profile profile_1 previous_grade coefficient 0.422 0.0736 5.74 9.73e-09
#> 24 group group_class_1 (Intercept) coefficient -0.763 0.2213 -3.45 5.64e-04
#> p_adjusted conf_low conf_high
#> 1 6.21e-88 0.486 0.593
#> 2 4.74e-40 0.365 0.492
#> 3 3.74e-138 0.571 0.668
#> 4 1.38e-72 0.473 0.588
#> 5 1.87e-186 0.605 0.692
#> 6 4.13e-134 -0.827 -0.705
#> 7 1.68e-86 -0.669 -0.548
#> 8 3.58e-248 -0.930 -0.828
#> 9 3.78e-162 -0.807 -0.698
#> 10 2.34e-266 -0.972 -0.868
#> 11 NA 0.536 0.652
#> 12 NA 0.764 0.928
#> 13 NA 0.437 0.534
#> 14 NA 0.623 0.759
#> 15 NA 0.348 0.425
#> 16 NA 0.467 0.593
#> 17 NA 0.478 0.604
#> 18 NA 0.321 0.410
#> 19 NA 0.374 0.477
#> 20 NA 0.333 0.422
#> 21 4.35e-20 -1.868 -1.211
#> 22 3.62e-41 1.082 1.452
#> 23 9.73e-09 0.278 0.567
#> 24 5.64e-04 -1.197 -0.329Staged estimation
To hold the profiles fixed while the group classes are estimated, we
call fit_staged() with the same arguments as
multilpa(). It first fits the profiles with one group
class, then estimates the group classes with the profile means and
variances held fixed.
staged <- fit_staged(course_engagement,
vars = c("browse", "lectures", "forum_read",
"forum_post", "attendance"),
id = "student", n_profiles = 2, n_group_classes = 2,
n_starts = 3, seed = 1)
staged
#> Two-level latent profile analysis: 2 profiles, 2 group classes
#> 1422 individuals in 106 groups; varying diagonal residual covariance (VVI)
#> Log likelihood: -8440.080258 | AIC: 16886.161 | BIC (groups): 16894.151
#> Converged: TRUE | iterations: 9 | best start: 2/3
#> Held fixed, not estimated here: means, variances (first stage)
#> Parameters estimated here: 3; with the held measurement: 23
#>
#> profile browse lectures forum_read forum_post attendance count proportion
#> 1 0.537 0.430 0.618 0.533 0.653 834 0.586
#> 2 -0.759 -0.608 -0.873 -0.753 -0.923 588 0.414
#>
#> Variances and standard errors: get_results(x, "profiles").
#> Every other table: get_results(x, what = ), or get_results(x, "all").To see what each stage estimated, we call get_results()
with "stages". parameters counts the
parameters estimated at a stage, and
parameters_with_measurement adds those held fixed.
get_results(staged, "stages")
#> stage group_classes fixed log_likelihood parameters parameters_with_measurement
#> 1 measurement 1 <NA> -8615 21 21
#> 2 membership 2 means, variances -8440 3 23
#> converged
#> 1 TRUE
#> 2 TRUEReference
The calls below list other functions for external variables and fixed parameters. Each comment states what the call returns.
get_results(fit, "classification_errors", level = "groups") # group-level errors
get_results(fit, "bch_weights", level = "individuals") # BCH weights per case
three_step(fit, data = course_engagement, outcome = "previous_grade",
method = "proportional") # weighting by posterior probabilities
three_step(fit, data = course_engagement, outcome = "previous_grade",
vcov_type = "independent") # standard errors without clustering
r3step(fit, data = course_engagement,
covariates = "previous_grade") # observed-information standard errors
measurement <- multilpa(course_engagement,
vars = c("browse", "lectures", "forum_read",
"forum_post", "attendance"),
id = "student", n_profiles = 2,
n_group_classes = 1, seed = 1) # a one-class measurement fit
fit_staged(course_engagement,
vars = c("browse", "lectures", "forum_read", "forum_post",
"attendance"),
id = "student", n_profiles = 2, n_group_classes = 2,
measurement = measurement) # stage two from that fit
multilpa(course_engagement,
vars = c("browse", "lectures", "forum_read", "forum_post",
"attendance"),
id = "student", n_profiles = 2, n_group_classes = 2,
start = starting_values(fit, what = "measurement"),
fixed = "variances") # hold chosen parameters fixed
parameter_inference(staged) # standard errors of the estimated stageAssumptions
- Three-step methods assume classification errors are independent of the external variable given the true class, and treat the error matrix as known.
- One-step coefficients are conditional on the group class and can change the profiles, so they need not equal the three-step coefficients.
- Staged estimation does not carry first-stage uncertainty into the second stage, and its likelihood is not the joint maximum.
- All associations are observational.