PLmixed: A package for estimating GLMMs with factor structures.

Description

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.

References

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.

Simulated multilevel IRT dataset.

Description

A simulated dataset that replicates the dataset from CITO.

Usage

IRTsim

Format

A data frame with 2500 rows and 4 variables:

sid

Student ID

school

School ID

item

Item ID

y

Response

Simulated Multi-rater Multi-response dataset.

Description

A simulated dataset that replicates the dataset from a multi-rater mult-reponse study where teachers and students provided responses about two student traits.

Usage

JUDGEsim

Format

A data frame with 54462 rows and 7 variables:

item

Item ID

method

1 = teacher response, 2 = student response

trait

1 = trait 1, 2 = trait 2

stu

Student ID

class

Classroom ID

tch

Teacher ID

response

Item response

Simulated KYPS item-level dataset.

Description

A simulated dataset that replicates the dataset item-level data from KYPS.

Usage

KYPSitemsim

Format

A data frame with 66947 rows and 6 variables:

id

Student ID

time

Time Identifier

item

Item ID

mid

Middle School ID

hid

High School ID

response

Item Response

Simulated KYPS dataset.

Description

A simulated dataset that replicates the dataset from KYPS.

Usage

KYPSsim

Format

A data frame with 11494 rows and 5 variables:

mid

Middle School ID

hid

High School ID

sid

Student ID

time

Time Identifier

esteem

Self Esteem

Fit GLMM with Factor Structure

Description

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.

Usage

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
)

Arguments

Details

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.

Value

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.

References

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.

See Also

lme4

glmer

lmer

Examples

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)

coef.PLmod

Description

Obtain coefficients for a model of class PLmod.

Usage

## S3 method for class 'PLmod'
coef(object, ...)

Arguments

Value

sum of the random and fixed effects coefficients for each explanatory variable for each level of the grouping factor.

fitted.PLmod

Description

Obtain fitted values for a model of class PLmod.

Usage

## S3 method for class 'PLmod'
fitted(object, ...)

Arguments

fixef.PLmod

Description

Obtain fixed effect estimates for a model of class PLmod.

Usage

## S3 method for class 'PLmod'
fixef(object, ...)

Arguments

iterPlot

Description

Plot parameter estimates at each optim iteration.

Usage

iterPlot(object)

Arguments

plot.PLmod

Description

Diagnostic plots for a model of class PLmod.

Usage

## S3 method for class 'PLmod'
plot(x, ...)

Arguments

predict.PLmod

Description

Predict response values from a model of class PLmod.

Usage

## S3 method for class 'PLmod'
predict(object, newdata = NULL, ...)

Arguments

print.PLmod

Description

Print the fitted PLmixed model object of class PLmod.

Usage

## S3 method for class 'PLmod'
print(x, digits = 4, ...)

Arguments

print.summary.PLmod

Description

Print the output for a PLmixed model object of class PLmod.

Usage

## S3 method for class 'summary.PLmod'
print(x, digits = 4, ...)

Arguments

ranef.PLmod

Description

Obtain conditional modes of the random effects for a model of class PLmod.

Usage

## S3 method for class 'PLmod'
ranef(object, ...)

Arguments

residuals.PLmod

Description

Obtain residuals for a model of class PLmod.

Usage

## S3 method for class 'PLmod'
residuals(object, ...)

Arguments

simulate.PLmod

Description

Simulate responses from a model of class PLmod.

Usage

## S3 method for class 'PLmod'
simulate(object, ...)

Arguments

summary.PLmod

Description

Obtain key output for a fitted PLmixed model object of class PLmod.

Usage

## S3 method for class 'PLmod'
summary(object, ...)

Arguments

Value

An object containing all parameter estimates and model characteristics.