PLmixed: Fit GLMM with Factor Structure

Description Usage Arguments Details Value See Also Examples

View source: R/PLmixed_main_function.R

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

1
2
3
4
5
6
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)

Arguments

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 factor argument.

data

A data frame containing the variables used in the model (but not factor names).

family

A GLM family, see glm and family.

load.var

A variable in the dataframe identifying what the factors load onto. Each unique element in load.var will have a unique factor loading. All rows in the dataset with the same value for load.var will have the same factor loading.

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 = 1) or vector of initial lambda values. If a scalar, the value is applied to all lambda parameters. If a vector, the values apply in row by column by matrix order.

nlp

A character vector containing the names of additional nonlinear parameters that are in the model formula.

init.nlp

A scalar (default = 1) or vector of initial nlp values. If a scalar, the value is applied to all nlp parameters. If a vector, the values apply in the order listed.

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 1, corresponding to the Laplace approximation. See glmer.

method

The optim optimization method. Defaults to L-BFGS-B.

lower

Lower bound on lambda parameters if method = "L-BFGS-B".

upper

Upper bound on lambda parameters if method = "L-BFGS-B".

lme4.optimizer

The lme4 optimization method.

lme4.start

Start values used for lme4.

lme4.optCtrl

A list controlling the lme4 optimization. See lmerControl or glmerControl

opt.control

Controls for the optim optimization.

REML

Use REML if model is linear? Defaults to FALSE.

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.

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.

See Also

lme4

glmer

lmer

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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)

PLmixed documentation built on May 3, 2018, 1:04 a.m.