| Title: | Estimate (Generalized) Linear Mixed Models with Factor Structures |
|---|---|
| Description: | Utilizes the 'lme4' and 'optimx' packages (previously the optim() function from 'stats') to estimate (generalized) linear mixed models (GLMM) with factor structures using a profile likelihood approach, as outlined in Jeon and Rabe-Hesketh (2012) <doi:10.3102/1076998611417628> and Rockwood and Jeon (2019) <doi:10.1080/00273171.2018.1516541>. Factor analysis and item response models can be extended to allow for an arbitrary number of nested and crossed random effects, making it useful for multilevel and cross-classified models. |
| Authors: | Minjeong Jeon [aut], Nicholas Rockwood [aut, cre] |
| Maintainer: | Nicholas Rockwood <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.7 |
| Built: | 2025-02-14 04:59:51 UTC |
| Source: | https://github.com/cran/PLmixed |
The PLmixed package's main function is PLmixed, which estimates
the model through nested maximizations using the lme4 and optimx packages
(previously the optim function). This extends the capabilities of lme4
to allow for estimated factor structures, making it useful for estimating multilevel
factor analysis and item response theory models with an arbitrary number of hierarchical
levels or crossed random effects.
Rockwood, N. J., & Jeon, M. (2019). Estimating complex measurement and growth models using the R package PLmixed.Multivariate Behavioral Research, 54(2), 288-306.
Jeon, M., & Rabe-Hesketh, S. (2012). Profile-likelihood approach for estimating generalized linear mixed models with factor structures. Journal of Educational and Behavioral Statistics, 37(4), 518-542.
Obtain coefficients for a model of class PLmod.
## S3 method for class 'PLmod' coef(object, ...)## S3 method for class 'PLmod' coef(object, ...)
object |
an object of class PLmod |
... |
Additional arguments from |
sum of the random and fixed effects coefficients for each explanatory variable for each level of the grouping factor.
Obtain fitted values for a model of class PLmod.
## S3 method for class 'PLmod' fitted(object, ...)## S3 method for class 'PLmod' fitted(object, ...)
object |
an object of class PLmod |
... |
Additional arguments from |
Obtain fixed effect estimates for a model of class PLmod.
## S3 method for class 'PLmod' fixef(object, ...)## S3 method for class 'PLmod' fixef(object, ...)
object |
an object of class PLmod |
... |
Additional arguments from |
A simulated dataset that replicates the dataset from CITO.
IRTsimIRTsim
A data frame with 2500 rows and 4 variables:
Student ID
School ID
Item ID
Response
A simulated dataset that replicates the dataset from a multi-rater mult-reponse study where teachers and students provided responses about two student traits.
JUDGEsimJUDGEsim
A data frame with 54462 rows and 7 variables:
Item ID
1 = teacher response, 2 = student response
1 = trait 1, 2 = trait 2
Student ID
Classroom ID
Teacher ID
Item response
A simulated dataset that replicates the dataset item-level data from KYPS.
KYPSitemsimKYPSitemsim
A data frame with 66947 rows and 6 variables:
Student ID
Time Identifier
Item ID
Middle School ID
High School ID
Item Response
A simulated dataset that replicates the dataset from KYPS.
KYPSsimKYPSsim
A data frame with 11494 rows and 5 variables:
Middle School ID
High School ID
Student ID
Time Identifier
Self Esteem
Fit a (generalized) linear mixed effects model (GLMM) with factor structures. Utilizes both the
lme4 package and optim function for estimation using a profile-likelihood based
approach.
PLmixed( formula, data, family = gaussian, load.var = NULL, lambda = NULL, factor = NULL, init = 1, nlp = NULL, init.nlp = 1, nAGQ = 1, method = "L-BFGS-B", lower = -Inf, upper = Inf, lme4.optimizer = "bobyqa", lme4.start = NULL, lme4.optCtrl = list(), opt.control = NULL, REML = FALSE, SE = 1, ND.method = "simple", check = "stop", est = TRUE, iter.count = TRUE )PLmixed( formula, data, family = gaussian, load.var = NULL, lambda = NULL, factor = NULL, init = 1, nlp = NULL, init.nlp = 1, nAGQ = 1, method = "L-BFGS-B", lower = -Inf, upper = Inf, lme4.optimizer = "bobyqa", lme4.start = NULL, lme4.optCtrl = list(), opt.control = NULL, REML = FALSE, SE = 1, ND.method = "simple", check = "stop", est = TRUE, iter.count = TRUE )
formula |
A formula following that of lme4, with the addition that factors can be specified
as random effects. Factor names should not be names of variables in the data set, and are instead
defined with the |
data |
A data frame containing the variables used in the model (but not factor names). |
family |
|
load.var |
A variable in the dataframe identifying what the factors load onto. Each unique element in |
lambda |
A matrix or list of matrices corresponding to the loading matrices. A value of NA indicates the loading is freely estimated, while a numeric entry indicates a constraint. |
factor |
A list of factors corresponding to the loading matrices and factors specified in model. |
init |
A scalar (default = |
nlp |
A character vector containing the names of additional nonlinear parameters that are in the model formula. |
init.nlp |
A scalar (default = |
nAGQ |
If family is non-gaussian, the number of points per axis for evaluating the adaptive
Gauss-Hermite approximation to the log-likelihood. Defaults to |
method |
The |
lower |
Lower bound on lambda parameters if applicable. |
upper |
Upper bound on lambda parameters if applicable. |
lme4.optimizer |
The lme4 optimization method. |
lme4.start |
Start values used for lme4. |
lme4.optCtrl |
A list controlling the lme4 optimization. See |
opt.control |
Controls for the |
REML |
Use REML if model is linear? Defaults to |
SE |
Method of calculating standard errors for fixed effects. |
ND.method |
Method of calculating numerical derivatives. |
check |
Check number of observations vs. levels and number of observations vd. random effects. |
est |
Return parameter estimates. |
iter.count |
Print the iteration counter during optimization. |
Factors are listed within the formula in the same way that random effects are specified
in lme4. The grouping variable listed after | defines what the factor values randomly
vary over, just as | does for other random effects. The names of factors and other random
effect terms can be listed within the same set of parentheses, allowing the covariance between the
factor(s) and random effect(s) to be estimated. The same factor may be specified for multiple grouping
variables, allowing for multilevel or crossed effects.
The factor argument must list any factor that appears in the formula. The ordering will
depend on the ordering of the matrices listed within lambda. The matrices in lambda
specify the factor loading matrices. The number of matrices in lambda should equal the number
of character vectors in factor and the number of elements in load.var. The number of
rows in the kth matrix listed in lambda should correspond to the number of unique elements
in the dataset for the kth variable listed in load.var, and the number of columns in the kth
matrix should correspond to the number of factors listed in the kth character vector of factor.
Within the kth matrix, the (i, j) cell corresponds to the factor loading for the ith unique element
of the kth variable listed in load.var on the jth factor listed in the kth character vector
of factor. Each element of the matrix should be either a number or NA. If the element is a
number, the loading will be constrained to that value. If the element is an NA, the loading will
be freely estimated. For identification, it is necessary (but not sufficient) for at least one element in
each column to be constrained.
The nlp argument can be viewed as a special case of the factor argument, where the character vector
listed in nlp is automatically linked to 1 x p lambda matrix, where p is the number of elements in nlp.
The load.var for these parameters is viewed as a constant, so that the nlp parameters are equivalent for
all rows in the dataset. Thus, nlp simplifies the process of adding additional nonlinear parameters to the model
without having to specify corresponding lambda and load.var values.
An object of class PLmod, which contains an object of class merMod as one of its elements.
Some functions for class merMod have been adapted to work with class PLmod. Others can be utilized
using object$'lme4 Model', where object is an object of class PLmod.
Rockwood, N. J., & Jeon, M. (2019). Estimating complex measurement and growth models using the R package PLmixed.Multivariate Behavioral Research, 54(2), 288-306.
Jeon, M., & Rabe-Hesketh, S. (2012). Profile-likelihood approach for estimating generalized linear mixed models with factor structures. Journal of Educational and Behavioral Statistics, 37(4), 518-542.
data("IRTsim") # Load the IRTsim data IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3 set.seed(12345) IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item irt.lam = c(1, NA, NA) # Specify the lambda matrix # Below, the # in front of family = binomial can be removed to change the response distribution # to binomial, where the default link function is logit. irt.model <- PLmixed(y ~ 0 + as.factor(item) + (0 + abil.sid |sid) +(0 + abil.sid |school), data = IRTsub, load.var = c("item"), # family = binomial, factor = list(c("abil.sid")), lambda = list(irt.lam)) summary(irt.model) ## Not run: # A more time-consuming example. # ~ 5-10 minutes data("KYPSsim") # Load the KYPSsim data kyps.lam <- rbind(c( 1, 0), # Specify the lambda matrix c(NA, 0), c(NA, 1), c(NA, NA)) kyps.model <- PLmixed(esteem ~ as.factor(time) + (0 + hs | hid) + (0 + ms | mid) + (1 | sid), data = KYPSsim, factor = list(c("ms", "hs")), load.var = c("time"), lambda = list(kyps.lam)) summary(kyps.model) data("JUDGEsim") JUDGEsim <- JUDGEsim[order(JUDGEsim$item), ] # Order by item unique(JUDGEsim$item) # Specify Lambda matrix judge.lam <- rbind(c( 1, 0, 1, 0, 0, 0), c(NA, 0, NA, 0, 0, 0), c(NA, 0, NA, 0, 0, 0), c( 0, 1, 0, 1, 0, 0), c( 0, NA, 0, NA, 0, 0), c( 0, NA, 0, NA, 0, 0), c( 0, 0, 0, 0, 1, 0), c( 0, 0, 0, 0, NA, 0), c( 0, 0, 0, 0, NA, 0), c( 0, 0, 0, 0, 0, 1), c( 0, 0, 0, 0, 0, NA), c( 0, 0, 0, 0, 0, NA)) # Conduct analysis judge.example <- PLmixed(response ~ 0 + as.factor(item) + (1 | class) + (0 + trait1.t + trait2.t + trait1.s + trait2.s | stu) + (0 + teacher1 + teacher2 | tch), data = JUDGEsim, lambda = list(judge.lam), load.var = "item", factor = list(c("teacher1", "teacher2", "trait1.t", "trait2.t", "trait1.s", "trait2.s"))) summary(judge.example) data("KYPSitemsim") time.lam <- rbind(c( 1, 0), # Specify time lambda matrix c(NA, 0), c(NA, 1), c(NA, NA)) item.lam <- c(1, NA, NA, NA, NA, NA) # Specify item lambda matrix KYPSitemsim$time2 <- (KYPSitemsim$time == 2) * 1 KYPSitemsim$time3 <- (KYPSitemsim$time == 3) * 1 KYPSitemsim$time4 <- (KYPSitemsim$time == 4) * 1 kyps.item.model <- PLmixed(response ~ 0 + as.factor(item) + lat.var:time2 + lat.var:time3 + lat.var:time4 + (0 + hs:lat.var | hid) + (0 + ms:lat.var | mid) + (0 + lat.var:as.factor(time) | id), data = KYPSitemsim, lambda = list(time.lam, item.lam), factor = list(c("ms", "hs"), "lat.var"), load.var = c("time", "item")) summary(kyps.item.model) ## End(Not run)data("IRTsim") # Load the IRTsim data IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3 set.seed(12345) IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item irt.lam = c(1, NA, NA) # Specify the lambda matrix # Below, the # in front of family = binomial can be removed to change the response distribution # to binomial, where the default link function is logit. irt.model <- PLmixed(y ~ 0 + as.factor(item) + (0 + abil.sid |sid) +(0 + abil.sid |school), data = IRTsub, load.var = c("item"), # family = binomial, factor = list(c("abil.sid")), lambda = list(irt.lam)) summary(irt.model) ## Not run: # A more time-consuming example. # ~ 5-10 minutes data("KYPSsim") # Load the KYPSsim data kyps.lam <- rbind(c( 1, 0), # Specify the lambda matrix c(NA, 0), c(NA, 1), c(NA, NA)) kyps.model <- PLmixed(esteem ~ as.factor(time) + (0 + hs | hid) + (0 + ms | mid) + (1 | sid), data = KYPSsim, factor = list(c("ms", "hs")), load.var = c("time"), lambda = list(kyps.lam)) summary(kyps.model) data("JUDGEsim") JUDGEsim <- JUDGEsim[order(JUDGEsim$item), ] # Order by item unique(JUDGEsim$item) # Specify Lambda matrix judge.lam <- rbind(c( 1, 0, 1, 0, 0, 0), c(NA, 0, NA, 0, 0, 0), c(NA, 0, NA, 0, 0, 0), c( 0, 1, 0, 1, 0, 0), c( 0, NA, 0, NA, 0, 0), c( 0, NA, 0, NA, 0, 0), c( 0, 0, 0, 0, 1, 0), c( 0, 0, 0, 0, NA, 0), c( 0, 0, 0, 0, NA, 0), c( 0, 0, 0, 0, 0, 1), c( 0, 0, 0, 0, 0, NA), c( 0, 0, 0, 0, 0, NA)) # Conduct analysis judge.example <- PLmixed(response ~ 0 + as.factor(item) + (1 | class) + (0 + trait1.t + trait2.t + trait1.s + trait2.s | stu) + (0 + teacher1 + teacher2 | tch), data = JUDGEsim, lambda = list(judge.lam), load.var = "item", factor = list(c("teacher1", "teacher2", "trait1.t", "trait2.t", "trait1.s", "trait2.s"))) summary(judge.example) data("KYPSitemsim") time.lam <- rbind(c( 1, 0), # Specify time lambda matrix c(NA, 0), c(NA, 1), c(NA, NA)) item.lam <- c(1, NA, NA, NA, NA, NA) # Specify item lambda matrix KYPSitemsim$time2 <- (KYPSitemsim$time == 2) * 1 KYPSitemsim$time3 <- (KYPSitemsim$time == 3) * 1 KYPSitemsim$time4 <- (KYPSitemsim$time == 4) * 1 kyps.item.model <- PLmixed(response ~ 0 + as.factor(item) + lat.var:time2 + lat.var:time3 + lat.var:time4 + (0 + hs:lat.var | hid) + (0 + ms:lat.var | mid) + (0 + lat.var:as.factor(time) | id), data = KYPSitemsim, lambda = list(time.lam, item.lam), factor = list(c("ms", "hs"), "lat.var"), load.var = c("time", "item")) summary(kyps.item.model) ## End(Not run)
Diagnostic plots for a model of class PLmod.
## S3 method for class 'PLmod' plot(x, ...)## S3 method for class 'PLmod' plot(x, ...)
x |
an object of class PLmod |
... |
Additional arguments from |
Predict response values from a model of class PLmod.
## S3 method for class 'PLmod' predict(object, newdata = NULL, ...)## S3 method for class 'PLmod' predict(object, newdata = NULL, ...)
object |
an object of class PLmod |
newdata |
data frame to obtain predictions for |
... |
Additional arguments from |
Print the fitted PLmixed model object of class PLmod.
## S3 method for class 'PLmod' print(x, digits = 4, ...)## S3 method for class 'PLmod' print(x, digits = 4, ...)
x |
an object of class PLmod |
digits |
minimal number of significant digits, see |
... |
Additional arguments. |
Print the output for a PLmixed model object of class PLmod.
## S3 method for class 'summary.PLmod' print(x, digits = 4, ...)## S3 method for class 'summary.PLmod' print(x, digits = 4, ...)
x |
an object of class PLmod |
digits |
minimal number of significant digits, see |
... |
Additional arguments. |
Obtain conditional modes of the random effects for a model of class PLmod.
## S3 method for class 'PLmod' ranef(object, ...)## S3 method for class 'PLmod' ranef(object, ...)
object |
an object of class PLmod |
... |
Additional arguments from |
Obtain residuals for a model of class PLmod.
## S3 method for class 'PLmod' residuals(object, ...)## S3 method for class 'PLmod' residuals(object, ...)
object |
an object of class PLmod |
... |
Additional arguments from |
Simulate responses from a model of class PLmod.
## S3 method for class 'PLmod' simulate(object, ...)## S3 method for class 'PLmod' simulate(object, ...)
object |
an object of class PLmod |
... |
Additional arguments from |
Obtain key output for a fitted PLmixed model object of class PLmod.
## S3 method for class 'PLmod' summary(object, ...)## S3 method for class 'PLmod' summary(object, ...)
object |
an object of class PLmod |
... |
Additional arguments. |
An object containing all parameter estimates and model characteristics.