,
Douglas Bates [ctb],
Ben Bolker [ctb],
Martin Maechler [ctb],
Allan Leal [ctb],
Fabian Scheipl [ctb],
Steven Walker [ctb],
Simon Wood [ctb]| Title: | Generalized Additive Latent and Mixed Models |
|---|---|
| Description: | Estimates generalized additive latent and mixed models using maximum marginal likelihood, as defined in Sorensen et al. (2023) <doi:10.1007/s11336-023-09910-z>, which is an extension of Rabe-Hesketh and Skrondal (2004)'s unifying framework for multilevel latent variable modeling <doi:10.1007/BF02295939>. Efficient computation is done using sparse matrix methods, Laplace approximation, and automatic differentiation. The framework includes generalized multilevel models with heteroscedastic residuals, mixed response types, factor loadings, smoothing splines, crossed random effects, and combinations thereof. Syntax for model formulation is close to 'lme4' (Bates et al. (2015) <doi:10.18637/jss.v067.i01>) and 'PLmixed' (Rockwood and Jeon (2019) <doi:10.1080/00273171.2018.1516541>). |
| Authors: | Øystein Sørensen [aut, cre] ,
Douglas Bates [ctb],
Ben Bolker [ctb],
Martin Maechler [ctb],
Allan Leal [ctb],
Fabian Scheipl [ctb],
Steven Walker [ctb],
Simon Wood [ctb] |
| Maintainer: | Øystein Sørensen <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.2.1.9000 |
| Built: | 2024-11-21 05:37:10 UTC |
| Source: | https://github.com/LCBC-UiO/galamm |
Anova function for comparing different GALAMMs fitted on the same data.
## S3 method for class 'galamm' anova(object, ...)## S3 method for class 'galamm' anova(object, ...)
object |
An object of class |
... |
Other fitted models of class |
A table with model comparison metric.
Some of the source code for this function is adapted from
lme4:::anova.merMod, with authors Douglas M. Bates, Martin Maechler,
Ben Bolker, and Steve Walker.
Bates DM, Mächler M, Bolker B, Walker S (2015). “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software, 67(1), 1–48. ISSN 1548-7660, doi:10.18637/jss.v067.i01.
summary.galamm() for the summary method and anova() for the
generic function.
Other summary functions:
plot.galamm(),
plot_smooth.galamm(),
print.galamm(),
print.summary.galamm(),
summary.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Model without interaction count_mod0 <- galamm( formula = y ~ lbas + treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Model comparison anova(count_mod, count_mod0)# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Model without interaction count_mod0 <- galamm( formula = y ~ lbas + treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Model comparison anova(count_mod, count_mod0)
Currently, this function only returns the fixed effects.
## S3 method for class 'galamm' coef(object, ...)## S3 method for class 'galamm' coef(object, ...)
object |
An object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
A matrix with the requested coefficients.
fixef.galamm() for fixed effects, ranef.galamm() for random
effects, and coef() for the generic function.
Other details of model fit:
VarCorr(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract coefficients coef(count_mod)# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract coefficients coef(count_mod)
Simulated dataset mimicking the measurement of abilities in three cognitive
domains. The latent traits (cognitive ability in a given domain) are based on
the functions in mgcv::gamSim
(Wood 2017), and depend on the
explanatory variable x.
cognitioncognition
cognition A data frame with 14400 rows and 7 columns:Subject ID.
Factor variable denoting the cognitive domain.
Explanatory variable.
Factor variable denoting the timepoint.
Factor variable denoting the item within the tests of each cognitive domain.
Number of trials, if applicable.
Response variable. For domain 1 a real number, for domain 2 a binomially distributed variable based on a single trial, for domain 3 a real number.
Wood SN (2017). Generalized Additive Models: An Introduction with R, 2 edition. Chapman and Hall/CRC.
Other datasets:
diet,
epilep,
hsced,
latent_covariates,
latent_covariates_long,
lifespan,
mresp,
mresp_hsced
Confidence intervals for model parameters
## S3 method for class 'galamm' confint(object, parm, level = 0.95, method = "Wald", ...)## S3 method for class 'galamm' confint(object, parm, level = 0.95, method = "Wald", ...)
object |
An object of class |
parm |
Parameters for which to compute intervals. Use |
level |
Decimal number specifying the confidence level. Defaults to 0.95. |
method |
Character of length one specifying the type of confidence interval. Currently only "Wald" is available. The argument is case sensitive. |
... |
Other arguments passed on to other methods. Currently not used. |
A matrix with the requested confidence intervals.
fixef.galamm() for fixed effects, coef.galamm() for coefficients
more generally, and vcov.galamm() for the variance-covariance matrix.
confint() is the generic function.
Other details of model fit:
VarCorr(),
coef.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) confint(count_mod, parm = "beta", level = .99)# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) confint(count_mod, parm = "beta", level = .99)
Extract deviance of galamm object
## S3 method for class 'galamm' deviance(object, ...)## S3 method for class 'galamm' deviance(object, ...)
object |
Object of class |
... |
Other arguments passed on to other methods. Currently not used. |
A numeric value giving the deviance of the model fit.
logLik.galamm() for a function returning the log likelihood and
deviance() for the generic function.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract deviance deviance(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract deviance deviance(mod)
Longitudinal epilepsy data from Morris et al. (1977). This documentation is based on Chapter 14.2 of Skrondal and Rabe-Hesketh (2004), where the dataset is used. See also Rabe-Hesketh et al. (2003).
dietdiet
diet A data frame with 236 rows and 7 columns:Subject ID.
Age (standardized).
Dummy variable indicating whether the subject is a bus driver or banking staff.
Integer indicating whether the outcome is fiber intake at time 1 (item = 1), fiber intake at time 2 (item = 2), or coronary heart disease (item = 3).
Outcome.
Dummy variable indicating whether y is an indicator for coronary heart disease, coded as 0/1.
Dummy variable indicating whether y is a fiber measurement at either timepoint 1 or 2.
Dummy variable indicating whether y is a fiber measurement at timepoint 2.
http://www.gllamm.org/books/readme.html#14.2
Morris JN, Marr JW, Clayton DG (1977).
“Diet and Heart: A Postscript.”
Br Med J, 2(6098), 1307–1314.
ISSN 0007-1447, 1468-5833, doi:10.1136/bmj.2.6098.1307.
Rabe-Hesketh S, Pickles A, Skrondal A (2003).
“Correcting for Covariate Measurement Error in Logistic Regression Using Nonparametric Maximum Likelihood Estimation.”
Statistical Modelling, 3(3), 215–232.
ISSN 1471-082X, doi:10.1191/1471082X03st056oa.
Skrondal A, Rabe-Hesketh S (2004).
Generalized Latent Variable Modeling, Interdisciplinary Statistics Series.
Chapman and Hall/CRC, Boca Raton, Florida.
Other datasets:
cognition,
epilep,
hsced,
latent_covariates,
latent_covariates_long,
lifespan,
mresp,
mresp_hsced
Longitudinal epilepsy data from Leppik et al. (1987). This documentation is based on Chapter 11.3 of Skrondal and Rabe-Hesketh (2004), where the dataset is used.
epilepepilep
epilep A data frame with 236 rows and 7 columns:Subject ID.
Number of seizures.
Dummy variable for treatment group.
Time at visit.
Dummy for visit 4.
Logarithm of age.
Logarithm of a quarter of the number of seizures in the eight weeks preceding entry into the trial.
http://www.gllamm.org/books/readme.html#11.3
Leppik IE, Dreifuss FE, Porter R, Bowman T, Santilli N, Jacobs M, Crosby C, Cloyd J, Stackman J, Graves N, Sutula T, Welty T, Vickery J, Brundage R, Gates J, Gumnit RJ, Gutierrez A (1987).
“A Controlled Study of Progabide in Partial Seizures: Methodology and Results.”
Neurology, 37(6), 963–963.
ISSN 0028-3878, 1526-632X, doi:10.1212/WNL.37.6.963.
Skrondal A, Rabe-Hesketh S (2004).
Generalized Latent Variable Modeling, Interdisciplinary Statistics Series.
Chapman and Hall/CRC, Boca Raton, Florida.
Other datasets:
cognition,
diet,
hsced,
latent_covariates,
latent_covariates_long,
lifespan,
mresp,
mresp_hsced
This function extracts parameter values from a fitted model object in a form that can be directly provided as initial values for a new model fit.
## S3 method for class 'galamm' extract_optim_parameters(object)## S3 method for class 'galamm' extract_optim_parameters(object)
object |
Object of class |
A list object containing the following elements:
theta Numerical vector of variance components, i.e., entries of
the lower Cholesky form of the covariance matrix of random effects.
beta Fixed regression coefficients.
lambda Factor loadings.
weights Weights for heteroscedastic residuals.
Other optimization functions:
galamm_control()
# Fit linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract parameters start <- extract_optim_parameters(mod) # Fit again using the Nelder-Mead algorithm, using start as initial values: mod_nm <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced, start = start, control = galamm_control(method = "Nelder-Mead") )# Fit linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract parameters start <- extract_optim_parameters(mod) # Fit again using the Nelder-Mead algorithm, using start as initial values: mod_nm <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced, start = start, control = galamm_control(method = "Nelder-Mead") )
Extract factor loadings from galamm object
## S3 method for class 'galamm' factor_loadings(object)## S3 method for class 'galamm' factor_loadings(object)
object |
Object of class |
This function has been named factor_loadings rather than just
loadings to avoid conflict with stats::loadings.
A matrix containing the estimated factor loadings with corresponding standard deviations.
The example for this function comes from PLmixed, with authors
Nicholas Rockwood and Minjeong Jeon
(Rockwood and Jeon 2019).
fixef.galamm() for fixed regression coefficients,
confint.galamm() for confidence intervals, and coef.galamm() for
coefficients more generally.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Logistic mixed model with factor loadings, example from PLmixed data("IRTsim", package = "PLmixed") # Reduce data size for the example to run faster IRTsub <- IRTsim[IRTsim$item < 4, ] IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] IRTsub$item <- factor(IRTsub$item) # Fix loading for first item to 1, and estimate the two others freely loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # Estimate model mod <- galamm(y ~ item + (0 + ability | sid) + (0 + ability | school), data = IRTsub, family = binomial, load.var = "item", factor = "ability", lambda = loading_matrix ) # Show estimated factor loadings, with standard errors factor_loadings(mod)# Logistic mixed model with factor loadings, example from PLmixed data("IRTsim", package = "PLmixed") # Reduce data size for the example to run faster IRTsub <- IRTsim[IRTsim$item < 4, ] IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] IRTsub$item <- factor(IRTsub$item) # Fix loading for first item to 1, and estimate the two others freely loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # Estimate model mod <- galamm(y ~ item + (0 + ability | sid) + (0 + ability | school), data = IRTsub, family = binomial, load.var = "item", factor = "ability", lambda = loading_matrix ) # Show estimated factor loadings, with standard errors factor_loadings(mod)
This function returns a list of families for an object of class
galamm, returned from galamm.
## S3 method for class 'galamm' family(object, ...)## S3 method for class 'galamm' family(object, ...)
object |
An object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
A list of family objects.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Mixed response model loading_matrix <- matrix(c(1, NA), ncol = 1) families <- c(gaussian, binomial) family_mapping <- ifelse(mresp$itemgroup == "a", 1, 2) mixed_resp <- galamm( formula = y ~ x + (0 + level | id), data = mresp, family = families, family_mapping = family_mapping, load.var = "itemgroup", lambda = loading_matrix, factor = "level" ) # This model has two family objects family(mixed_resp)# Mixed response model loading_matrix <- matrix(c(1, NA), ncol = 1) families <- c(gaussian, binomial) family_mapping <- ifelse(mresp$itemgroup == "a", 1, 2) mixed_resp <- galamm( formula = y ~ x + (0 + level | id), data = mresp, family = families, family_mapping = family_mapping, load.var = "itemgroup", lambda = loading_matrix, factor = "level" ) # This model has two family objects family(mixed_resp)
Extracts fitted values from a model including random effects.
## S3 method for class 'galamm' fitted(object, ...)## S3 method for class 'galamm' fitted(object, ...)
object |
An object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
A numerical vector with fit values for each row in the input data.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract fitted values and plot against x plot(hsced$x, fitted(mod))# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract fitted values and plot against x plot(hsced$x, fitted(mod))
Extract the fixed regression coefficients.
## S3 method for class 'galamm' fixef(object, ...)## S3 method for class 'galamm' fixef(object, ...)
object |
An object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
A named numeric vector containing the requested fixed effects.
ranef.galamm() for random effects, coef.galamm() for
coefficients more generally, and confint.galamm() for confidence intervals.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract fixed effects fixef(count_mod)# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract fixed effects fixef(count_mod)
Extract formula from fitted galamm object
## S3 method for class 'galamm' formula(x, ...)## S3 method for class 'galamm' formula(x, ...)
x |
Object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
The formula used to fit the model.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Mixed response model ------------------------------------------------------ # The mresp dataset contains a mix of binomial and Gaussian responses. # We need to estimate a factor loading which scales the two response types. loading_matrix <- matrix(c(1, NA), ncol = 1) # Define mapping to families. families <- c(gaussian, binomial) family_mapping <- ifelse(mresp$itemgroup == "a", 1, 2) # Fit the model mod <- galamm( formula = y ~ x + (0 + level | id), data = mresp, family = families, family_mapping = family_mapping, factor = "level", load.var = "itemgroup", lambda = loading_matrix ) # Formula formula(mod)# Mixed response model ------------------------------------------------------ # The mresp dataset contains a mix of binomial and Gaussian responses. # We need to estimate a factor loading which scales the two response types. loading_matrix <- matrix(c(1, NA), ncol = 1) # Define mapping to families. families <- c(gaussian, binomial) family_mapping <- ifelse(mresp$itemgroup == "a", 1, 2) # Fit the model mod <- galamm( formula = y ~ x + (0 + level | id), data = mresp, family = families, family_mapping = family_mapping, factor = "level", load.var = "itemgroup", lambda = loading_matrix ) # Formula formula(mod)
This function fits a generalized additive latent and mixed model
(GALAMMs), as described in
Sørensen et al. (2023).
The building blocks of these models are generalized additive mixed models
(GAMMs) (Wood 2017), of which
generalized linear mixed models
(Breslow and Clayton 1993; Harville 1977; Henderson 1975; Laird and Ware 1982)
are special cases. GALAMMs extend upon GAMMs by allowing factor structures,
as commonly used to model hypothesized latent traits underlying observed
measurements. In this sense, GALAMMs are an extension of generalized linear
latent and mixed models (GLLAMMs)
(Skrondal and Rabe-Hesketh 2004; Rabe-Hesketh et al. 2004)
which allows semiparametric estimation. The implemented algorithm used to
compute model estimates is described in
Sørensen et al. (2023),
and is an extension of the algorithm used for fitting generalized linear
mixed models by the lme4 package
(Bates et al. 2015). The syntax used to
define factor structures is based on that used by the PLmixed
package, which is detailed in
Rockwood and Jeon (2019).
galamm( formula, weights = NULL, data, family = gaussian, family_mapping = rep(1, nrow(data)), load.var = NULL, lambda = NULL, factor = NULL, factor_interactions = NULL, na.action = getOption("na.action"), start = NULL, control = galamm_control() )galamm( formula, weights = NULL, data, family = gaussian, family_mapping = rep(1, nrow(data)), load.var = NULL, lambda = NULL, factor = NULL, factor_interactions = NULL, na.action = getOption("na.action"), start = NULL, control = galamm_control() )
formula |
A formula specifying the model. Smooth terms are defined in
the style of the |
weights |
An optional formula object specifying an expression for the
residual variance. Defaults to |
data |
A data.frame containing all the variables specified by the model formula, with the exception of factor loadings. |
family |
A a list or character vector containing one or more model
families. For each element in |
family_mapping |
Optional vector mapping from the elements of
|
load.var |
Optional character specifying the name of the variable in
|
lambda |
Optional factor loading matrix. Numerical values indicate that
the given value is fixed, while
|
factor |
Optional character vector whose |
factor_interactions |
Optional list of length equal to the number of
columns in |
na.action |
Character of length one specifying a function which
indicates what should happen when the data contains |
start |
Optional named list of starting values for parameters. Possible
names of list elements are |
control |
Optional control object for the optimization procedure of
class |
A model object of class galamm, containing the following
elements:
call the matched call used when fitting the model.
random_effects a list containing the following two elements:
b random effects in original parametrization.
u random effects standardized to have identity covariance matrix.
model a list with various elements related to the model setup and
fit:
deviance deviance of final model.
deviance_residuals deviance residuals of the final model.
df degrees of freedom of model.
family a list of one or more family objects, as specified in the
family arguments to galamm.
factor_interactions List of formulas specifying interactions
between latent and observed variables, as provided to the argument
factor_interactions to galamm. If not provided, it is
NULL.
fit a numeric vector with fitted values.
fit_population a numeric vector with fitted values excluding
random effects.
hessian Hessian matrix of final model, i.e., the second
derivative of the log-likelihood with respect to all model parameters.
lmod Linear model object returned by lme4::lFormula, which
is used internally for setting up the models.
loglik Log-likelihood of final model.
n Number of observations.
pearson_residual Pearson residuals of final model.
reduced_hessian Logical specifying whether the full Hessian
matrix was computed, or a Hessian matrix with derivatives only with respect
to beta and lambda.
response A numeric vector containing the response values used when
fitting the model.
weights_object Object with weights used in model fitting. Is
NULL when no weights were used.
parameters A list object with model parameters and related
information:
beta_inds Integer vector specifying the indices of fixed
regression coefficients among the estimated model parameters.
dispersion_parameter One or more dispersion parameters of the
final model.
lambda_dummy Dummy matrix of factor loadings, which shows the
structure of the loading matrix that was supplied in the lambda
arguments.
lambda_inds Integer vector specifying the indices of factor
loadings among the estimated model parameters.
lambda_interaction_inds Integer vector specifying the indices
of regression coefficients for interactions between latent and observed
variables.
parameter_estimates Numeric vector of final parameter estimates.
parameter_names Names of all parameters estimates.
theta_inds Integer vector specifying the indices of variance
components among the estimated model parameters. Technically these are the
entries of the Cholesky decomposition of the covariance matrix.
weights_inds Integer vector specifying the indices of estimated
weights (used in heteroscedastic Gaussian models) among the estimated model
parameters.
gam List containing information about smooth terms in the model.
If no smooth terms are contained in the model, then it is a list of length
zero.
Bates DM, Mächler M, Bolker B, Walker S (2015).
“Fitting Linear Mixed-Effects Models Using Lme4.”
Journal of Statistical Software, 67(1), 1–48.
ISSN 1548-7660, doi:10.18637/jss.v067.i01.
Breslow NE, Clayton DG (1993).
“Approximate Inference in Generalized Linear Mixed Models.”
Journal of the American Statistical Association, 88(421), 9–25.
ISSN 0162-1459, doi:10.2307/2290687.
Harville DA (1977).
“Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems.”
Journal of the American Statistical Association, 72(358), 320–338.
ISSN 0162-1459, doi:10.2307/2286796.
Henderson CR (1975).
“Best Linear Unbiased Estimation and Prediction under a Selection Model.”
Biometrics, 31(2), 423–447.
ISSN 0006-341X, doi:10.2307/2529430.
Laird NM, Ware JH (1982).
“Random-Effects Models for Longitudinal Data.”
Biometrics, 38(4), 963–974.
ISSN 0006-341X, doi:10.2307/2529876.
Rabe-Hesketh S, Skrondal A, Pickles A (2004).
“Generalized Multilevel Structural Equation Modeling.”
Psychometrika, 69(2), 167–190.
ISSN 1860-0980, doi:10.1007/BF02295939.
Rockwood NJ, Jeon M (2019).
“Estimating Complex Measurement and Growth Models Using the R Package PLmixed.”
Multivariate Behavioral Research, 54(2), 288–306.
ISSN 0027-3171, doi:10.1080/00273171.2018.1516541.
Skrondal A, Rabe-Hesketh S (2004).
Generalized Latent Variable Modeling, Interdisciplinary Statistics Series.
Chapman and Hall/CRC, Boca Raton, Florida.
Sørensen Ø, Fjell AM, Walhovd KB (2023).
“Longitudinal Modeling of Age-Dependent Latent Traits with Generalized Additive Latent and Mixed Models.”
Psychometrika, 88(2), 456–486.
ISSN 1860-0980, doi:10.1007/s11336-023-09910-z.
Wood SN (2017).
Generalized Additive Models: An Introduction with R, 2 edition.
Chapman and Hall/CRC.
Other modeling functions:
s(),
t2()
# Mixed response model ------------------------------------------------------ # The mresp dataset contains a mix of binomial and Gaussian responses. # We need to estimate a factor loading which scales the two response types. loading_matrix <- matrix(c(1, NA), ncol = 1) # Define mapping to families. families <- c(gaussian, binomial) family_mapping <- ifelse(mresp$itemgroup == "a", 1, 2) # Fit the model mod <- galamm( formula = y ~ x + (0 + level | id), data = mresp, family = families, family_mapping = family_mapping, factor = "level", load.var = "itemgroup", lambda = loading_matrix ) # Summary information summary(mod) # Heteroscedastic model ----------------------------------------------------- # Residuals allowed to differ according to the item variable # We also set the initial value of the random intercept standard deviation # to 1 mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced, start = list(theta = 1) ) summary(mod) # Generalized additive mixed model with factor structures ------------------- # The cognition dataset contains simulated measurements of three latent # time-dependent processes, corresponding to individuals' abilities in # cognitive domains. We focus here on the first domain, and take a single # random timepoint per person: dat <- subset(cognition, domain == 1) dat <- split(dat, f = dat$id) dat <- lapply(dat, function(x) x[x$timepoint %in% sample(x$timepoint, 1), ]) dat <- do.call(rbind, dat) dat$item <- factor(dat$item) # At each timepoint there are three items measuring ability in the cognitive # domain. We fix the factor loading for the first measurement to one, and # estimate the remaining two. This is specified in the loading matrix. loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # We can now estimate the model. mod <- galamm( formula = y ~ 0 + item + sl(x, factor = "loading") + (0 + loading | id), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # We can plot the estimated smooth term plot_smooth(mod, shade = TRUE) # Interaction between observed and latent covariates ------------------------ # Define the loading matrix lambda <- matrix(c(1, NA, NA), ncol = 1) # Define the regression functions, one for each row in the loading matrix factor_interactions <- list(~1, ~1, ~x) # Fit the model mod <- galamm( formula = y ~ type + x:response + (0 + loading | id), data = latent_covariates, load.var = "type", lambda = lambda, factor = "loading", factor_interactions = factor_interactions ) # The summary output now include an interaction between the latent variable # and x, for predicting the third element in "type" summary(mod)# Mixed response model ------------------------------------------------------ # The mresp dataset contains a mix of binomial and Gaussian responses. # We need to estimate a factor loading which scales the two response types. loading_matrix <- matrix(c(1, NA), ncol = 1) # Define mapping to families. families <- c(gaussian, binomial) family_mapping <- ifelse(mresp$itemgroup == "a", 1, 2) # Fit the model mod <- galamm( formula = y ~ x + (0 + level | id), data = mresp, family = families, family_mapping = family_mapping, factor = "level", load.var = "itemgroup", lambda = loading_matrix ) # Summary information summary(mod) # Heteroscedastic model ----------------------------------------------------- # Residuals allowed to differ according to the item variable # We also set the initial value of the random intercept standard deviation # to 1 mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced, start = list(theta = 1) ) summary(mod) # Generalized additive mixed model with factor structures ------------------- # The cognition dataset contains simulated measurements of three latent # time-dependent processes, corresponding to individuals' abilities in # cognitive domains. We focus here on the first domain, and take a single # random timepoint per person: dat <- subset(cognition, domain == 1) dat <- split(dat, f = dat$id) dat <- lapply(dat, function(x) x[x$timepoint %in% sample(x$timepoint, 1), ]) dat <- do.call(rbind, dat) dat$item <- factor(dat$item) # At each timepoint there are three items measuring ability in the cognitive # domain. We fix the factor loading for the first measurement to one, and # estimate the remaining two. This is specified in the loading matrix. loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # We can now estimate the model. mod <- galamm( formula = y ~ 0 + item + sl(x, factor = "loading") + (0 + loading | id), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # We can plot the estimated smooth term plot_smooth(mod, shade = TRUE) # Interaction between observed and latent covariates ------------------------ # Define the loading matrix lambda <- matrix(c(1, NA, NA), ncol = 1) # Define the regression functions, one for each row in the loading matrix factor_interactions <- list(~1, ~1, ~x) # Fit the model mod <- galamm( formula = y ~ type + x:response + (0 + loading | id), data = latent_covariates, load.var = "type", lambda = lambda, factor = "loading", factor_interactions = factor_interactions ) # The summary output now include an interaction between the latent variable # and x, for predicting the third element in "type" summary(mod)
This function can be called for controling the optimization
procedure used when fitting GALAMMs using galamm.
galamm_control( optim_control = list(), method = c("L-BFGS-B", "Nelder-Mead"), maxit_conditional_modes = 10, pirls_tol_abs = 0.01, reduced_hessian = FALSE )galamm_control( optim_control = list(), method = c("L-BFGS-B", "Nelder-Mead"), maxit_conditional_modes = 10, pirls_tol_abs = 0.01, reduced_hessian = FALSE )
optim_control |
List containing optimization parameters. If |
method |
Character string defining the algorithm to be used for
maximizing the marginal log-likelihood. The default is |
maxit_conditional_modes |
Maximum number of iterations in penalized
iteratively reweighted least squares algorithm. Ignored if |
pirls_tol_abs |
Absolute convergence criterion for penalized iteratively reweighted least squares algorithm. Defaults to 0.01, which means that when the reduction in marginal likelihood between two iterations is below 0.01, the iterations stop. |
reduced_hessian |
Logical value. Defaults to |
Object of class galamm_control, which typically will be
provided as an argument to galamm.
Bates DM, Mächler M, Bolker B, Walker S (2015). “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software, 67(1), 1–48. ISSN 1548-7660, doi:10.18637/jss.v067.i01.
BROYDEN CG (1970). “The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations.” IMA Journal of Applied Mathematics, 6(1), 76–90. ISSN 0272-4960, doi:10.1093/imamat/6.1.76.
Byrd RH, Lu P, Nocedal J, Zhu C (1995). “A Limited Memory Algorithm for Bound Constrained Optimization.” SIAM Journal on Scientific Computing, 16(5), 1190–1208. ISSN 1064-8275, doi:10.1137/0916069.
Fletcher R (1970). “A New Approach to Variable Metric Algorithms.” The Computer Journal, 13(3), 317–322. ISSN 0010-4620, doi:10.1093/comjnl/13.3.317.
Goldfarb D (1970). “A Family of Variable-Metric Methods Derived by Variational Means.” Mathematics of Computation, 24(109), 23–26. ISSN 0025-5718, 1088-6842, doi:10.1090/S0025-5718-1970-0258249-6.
Nelder JA, Mead R (1965). “A Simplex Method for Function Minimization.” The Computer Journal, 7(4), 308–313. ISSN 0010-4620, doi:10.1093/comjnl/7.4.308.
Shanno DF (1970). “Conditioning of Quasi-Newton Methods for Function Minimization.” Mathematics of Computation, 24(111), 647–656. ISSN 0025-5718, 1088-6842, doi:10.1090/S0025-5718-1970-0274029-X.
Other optimization functions:
extract_optim_parameters.galamm()
# Define control object with quite a high degree of verbosity (trace = 6) # and using the last 20 BFGS updates to estimate the Hessian in L-BFGS-B. control <- galamm_control(optim_control = list(trace = 6, lmm = 20))# Define control object with quite a high degree of verbosity (trace = 6) # and using the last 20 BFGS updates to estimate the Hessian in L-BFGS-B. control <- galamm_control(optim_control = list(trace = 6, lmm = 20))
Simulated dataset with residual standard deviation that varies between items.
hscedhsced
hsced A data frame with 1200 rows and 5 columns:Subject ID.
Timepoint.
Item indicator.
Explanatory variable
Outcome.
There are no references for Rd macro \insertAllCites on this help page.
Other datasets:
cognition,
diet,
epilep,
latent_covariates,
latent_covariates_long,
lifespan,
mresp,
mresp_hsced
Simulated dataset for use in examples and testing with a latent covariate interacting with an observed covariate.
latent_covariateslatent_covariates
latent_covariates A data frame with 600 rows and 5 columns:Subject ID.
Type of observation in the y variable.
If it equals "measurement1" or "measurement2" then the
observation is a measurement of the latent variable. If it equals
"response", then the observation is the actual response.
Explanatory variable.
Observed response. Note, this includes both the actual response, and the measurements of the latent variable, since mathematically they are all treated as responses.
Dummy variable indicating whether the given row is a response or not.
Other datasets:
cognition,
diet,
epilep,
hsced,
latent_covariates_long,
lifespan,
mresp,
mresp_hsced
Simulated dataset for use in examples and testing with a latent covariate interacting with an observed covariate. In this data, each response has been measured six times for each subject.
latent_covariates_longlatent_covariates_long
latent_covariates_long A data frame with 800 rows and 5columns:
Subject ID.
Type of observation in the y variable.
If it equals "measurement1" or "measurement2" then the
observation is a measurement of the latent variable. If it equals
"response", then the observation is the actual response.
Explanatory variable.
Observed response. Note, this includes both the actual response, and the measurements of the latent variable, since mathematically they are all treated as responses.
Dummy variable indicating whether the given row is a response or not.
Other datasets:
cognition,
diet,
epilep,
hsced,
latent_covariates,
lifespan,
mresp,
mresp_hsced
This dataset is simulated based on the data used in Section 4.1 of Sørensen et al. (2023).
lifespanlifespan
lifespan A data frame with 54,457 rows and 7 columns:Subject ID.
Cognitive domain being measured. One of "epmem" for
episodic memory, "wmem" for working memory and "execfun" for
executive function/speed.
Integer indicating the timepoint number.
Age of participant at the timepoint.
The particular test at this observation.
Response. For "epmem" and "wmem" this is the number
of successes in 16 trials. For "execfun" it is the time in seconds
to complete the task.
Integer indicating whether the participant has taken the test at a previous timepoint.
Dummy variable for domain=="epmem".
Dummy variable for domain=="wmem".
Dummy variable for domain=="execfun".
Sørensen Ø, Fjell AM, Walhovd KB (2023). “Longitudinal Modeling of Age-Dependent Latent Traits with Generalized Additive Latent and Mixed Models.” Psychometrika, 88(2), 456–486. ISSN 1860-0980, doi:10.1007/s11336-023-09910-z.
Other datasets:
cognition,
diet,
epilep,
hsced,
latent_covariates,
latent_covariates_long,
mresp,
mresp_hsced
Extract Log-Likelihood of galamm Object
## S3 method for class 'galamm' logLik(object, ...)## S3 method for class 'galamm' logLik(object, ...)
object |
Object |
... |
Other arguments |
Object of class logLik
deviance.galamm() for a function returning deviance and
logLik() for the generic function.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract log likelihood logLik(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract log likelihood logLik(mod)
Very basic mixed response dataset with one set of normally distributed responses and one set of binomially distributed responses.
mrespmresp
mresp A data frame with 4000 rows and 5 columns:Subject ID.
Predictor variable.
Response.
Factor variable which equals "a" for the normally distributed responses and "b" for the binomially distributed response (with 1 trial).
Other datasets:
cognition,
diet,
epilep,
hsced,
latent_covariates,
latent_covariates_long,
lifespan,
mresp_hsced
Mixed response dataset with one set of normally distributed responses and one set of binomially distributed responses. The normally distributed response follow two different residual standard deviations.
mresp_hscedmresp_hsced
mresp A data frame with 4000 rows and 5 columns:Subject ID.
Predictor variable.
Response.
Factor variable which equals "a" for the normally distributed responses and "b" for the binomially distributed response (with 1 trial).
Grouping variable denoting which of the two residual standard deviations apply. Only relevant for the normally distributed responses.
Dummy variable indicating whether the observation on the given line is normally (Gaussian) distributed or not.
Other datasets:
cognition,
diet,
epilep,
hsced,
latent_covariates,
latent_covariates_long,
lifespan,
mresp
Extract the Number of Observations from a galamm Fit
## S3 method for class 'galamm' nobs(object, ...)## S3 method for class 'galamm' nobs(object, ...)
object |
An object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
A number
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Example model from lme4 data(sleepstudy, package = "lme4") fm1 <- galamm(Reaction ~ Days + (Days | Subject), data = sleepstudy) # There are 180 observations, which matches the number of rows in sleepstudy nobs(fm1) nrow(sleepstudy)# Example model from lme4 data(sleepstudy, package = "lme4") fm1 <- galamm(Reaction ~ Days + (Days | Subject), data = sleepstudy) # There are 180 observations, which matches the number of rows in sleepstudy nobs(fm1) nrow(sleepstudy)
Plots smooth terms of a fitted galamm object. This function is a thin
wrapper around mgcv::plot.gam
(Wood 2017).
## S3 method for class 'galamm' plot_smooth(object, ...)## S3 method for class 'galamm' plot_smooth(object, ...)
object |
Object of class |
... |
Other optional arguments, passed on to |
A plot is displayed on the screen.
Wood SN (2017). Generalized Additive Models: An Introduction with R, 2 edition. Chapman and Hall/CRC.
Other summary functions:
anova.galamm(),
plot.galamm(),
print.galamm(),
print.summary.galamm(),
summary.galamm()
# Generalized additive mixed model with factor structures ------------------- # The cognition dataset contains simulated measurements of three latent # time-dependent processes, corresponding to individuals' abilities in # cognitive domains. We focus here on the first domain, and take a single # random timepoint per person: dat <- subset(cognition, domain == 1) dat <- split(dat, f = dat$id) dat <- lapply(dat, function(x) x[x$timepoint %in% sample(x$timepoint, 1), ]) dat <- do.call(rbind, dat) dat$item <- factor(dat$item) # At each timepoint there are three items measuring ability in the cognitive # domain. We fix the factor loading for the first measurement to one, and # estimate the remaining two. This is specified in the loading matrix. loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # We can now estimate the model. mod <- galamm( formula = y ~ 0 + item + sl(x, factor = "loading") + (0 + loading | id), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # We can plot the estimated smooth term plot_smooth(mod, shade = TRUE) # We can turn off the rug at the bottom plot_smooth(mod, shade = TRUE, rug = FALSE)# Generalized additive mixed model with factor structures ------------------- # The cognition dataset contains simulated measurements of three latent # time-dependent processes, corresponding to individuals' abilities in # cognitive domains. We focus here on the first domain, and take a single # random timepoint per person: dat <- subset(cognition, domain == 1) dat <- split(dat, f = dat$id) dat <- lapply(dat, function(x) x[x$timepoint %in% sample(x$timepoint, 1), ]) dat <- do.call(rbind, dat) dat$item <- factor(dat$item) # At each timepoint there are three items measuring ability in the cognitive # domain. We fix the factor loading for the first measurement to one, and # estimate the remaining two. This is specified in the loading matrix. loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # We can now estimate the model. mod <- galamm( formula = y ~ 0 + item + sl(x, factor = "loading") + (0 + loading | id), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # We can plot the estimated smooth term plot_smooth(mod, shade = TRUE) # We can turn off the rug at the bottom plot_smooth(mod, shade = TRUE, rug = FALSE)
Diagnostic plots for galamm objects
## S3 method for class 'galamm' plot(x, ...)## S3 method for class 'galamm' plot(x, ...)
x |
An object of class |
... |
Optional arguments passed on to the |
A plot is displayed.
residuals.galamm() for extracting residuals and plot() for the
generic function.
Other summary functions:
anova.galamm(),
plot_smooth.galamm(),
print.galamm(),
print.summary.galamm(),
summary.galamm()
# Linear mixed model example from lme4 data("sleepstudy", package = "lme4") mod <- galamm(Reaction ~ Days + (Days | Subject), data = sleepstudy) # Diagnostic plot plot(mod)# Linear mixed model example from lme4 data("sleepstudy", package = "lme4") mod <- galamm(Reaction ~ Days + (Days | Subject), data = sleepstudy) # Diagnostic plot plot(mod)
Predictions are given at the population level, i.e., with random
effects set to zero. For fitted models including random effects, see
fitted.galamm. For mixed response models, only predictions on
the scale of the linear predictors is supported.
## S3 method for class 'galamm' predict(object, newdata = NULL, type = c("link", "response"), ...)## S3 method for class 'galamm' predict(object, newdata = NULL, type = c("link", "response"), ...)
object |
An object of class |
newdata |
Data from for which to evaluate predictions, in a
|
type |
Character argument specifying the type of prediction object to be returned. Case sensitive. |
... |
Optional arguments passed on to other methods. Currently used for
models with smooth terms, for which these arguments are forwarded to
|
A numeric vector of predicted values.
fitted.galamm() for model fits, residuals.galamm() for
residuals, and predict() for the generic function.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Plot response versus link: plot( predict(count_mod, type = "link"), predict(count_mod, type = "response") ) # Predict on a new dataset nd <- data.frame(lbas = c(.3, .2), treat = c(0, 1), lage = 0.2, v4 = -.2) predict(count_mod, newdata = nd) predict(count_mod, newdata = nd, type = "response")# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Plot response versus link: plot( predict(count_mod, type = "link"), predict(count_mod, type = "response") ) # Predict on a new dataset nd <- data.frame(lbas = c(.3, .2), treat = c(0, 1), lage = 0.2, v4 = -.2) predict(count_mod, newdata = nd) predict(count_mod, newdata = nd, type = "response")
Print method for GALAMM fits
## S3 method for class 'galamm' print(x, ...)## S3 method for class 'galamm' print(x, ...)
x |
An object of class |
... |
Further arguments passed on to other methods. Currently not used. |
Summary printed to screen. Invisibly returns the argument x.
summary.galamm() for the summary function and print() for the
generic.
Other summary functions:
anova.galamm(),
plot.galamm(),
plot_smooth.galamm(),
print.summary.galamm(),
summary.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) print(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) print(mod)
Print method for summary GALAMM fits
## S3 method for class 'summary.galamm' print(x, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'summary.galamm' print(x, digits = max(3, getOption("digits") - 3), ...)
x |
An object of class |
digits |
Number of digits to present in outputs. |
... |
Further arguments passed on to other methods. Currently used by
|
Summary printed to screen. Invisibly returns the argument x.
Some of the code for producing summary information has been derived
from the summary methods of mgcv (author: Simon Wood) and
lme4 (Bates et al. 2015)
(authors: Douglas M. Bates, Martin Maechler, Ben Bolker, and Steve Walker).
There are no references for Rd macro \insertAllCites on this help page.
summary.galamm() for the summary function and print() for the
generic function.
Other summary functions:
anova.galamm(),
plot.galamm(),
plot_smooth.galamm(),
print.galamm(),
summary.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) summary(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) summary(mod)
Print method for variance-covariance objects
## S3 method for class 'VarCorr.galamm' print( x, digits = max(3, getOption("digits") - 2), comp = c("Std.Dev.", "Variance"), corr = any(comp == "Std.Dev."), ... )## S3 method for class 'VarCorr.galamm' print( x, digits = max(3, getOption("digits") - 2), comp = c("Std.Dev.", "Variance"), corr = any(comp == "Std.Dev."), ... )
x |
An object of class |
digits |
Optional arguments specifying number of digits to use when printing. |
comp |
Character vector of length 1 or 2 specifying which variance components to print. Case sensitive. Can take one of the values "Std.Dev." and "Variance". |
corr |
Logical value indicating whether covariances or correlations should be printed. |
... |
Optional arguments passed on to other methods. Currently not used. |
The variance-covariance information is printed to the console and the
argument x is silently returned.
This function is derived from lme4:::print.VarCorr.merMod
written by Douglas M. Bates, Martin Maechler, Ben Bolker, and Steve Walker.
Bates DM, Mächler M, Bolker B, Walker S (2015). “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software, 67(1), 1–48. ISSN 1548-7660, doi:10.18637/jss.v067.i01.
VarCorr.galamm() for the function creating the variance-covariance
objects.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract information on variance and covariance VarCorr(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract information on variance and covariance VarCorr(mod)
Extract random effects from galamm object.
## S3 method for class 'galamm' ranef(object, ...)## S3 method for class 'galamm' ranef(object, ...)
object |
An object of class |
... |
Optional parameters passed on to other methods. Currently not used. |
An object of class ranef.galamm, containing the requested
random effects.
This function is derived from lme4::ranef.merMod, written by
Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker.
Bates DM, Mächler M, Bolker B, Walker S (2015). “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software, 67(1), 1–48. ISSN 1548-7660, doi:10.18637/jss.v067.i01.
fixef.galamm() for fixed effects and coef.galamm() for
coefficients more generally.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract random effects ranef(count_mod)# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract random effects ranef(count_mod)
Residuals of galamm objects
## S3 method for class 'galamm' residuals(object, type = c("pearson", "deviance"), ...)## S3 method for class 'galamm' residuals(object, type = c("pearson", "deviance"), ...)
object |
An object of class |
type |
Character of length one describing the type of residuals to be
returned. One of |
... |
Optional arguments passed on to other methods. Currently not used. |
Numeric vector of residual values.
fitted.galamm() for model fitted values, predict.galamm() for
model predictions, and plot.galamm() for diagnostic plots. The generic
function is residuals().
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract residuals residuals(count_mod)# Poisson GLMM count_mod <- galamm( formula = y ~ lbas * treat + lage + v4 + (1 | subj), data = epilep, family = poisson ) # Extract residuals residuals(count_mod)
Extracts response values from a model.
response(object, ...)response(object, ...)
object |
An object of class |
... |
Optional arguments passed on to other methods. Currently not used. |
A numerical vector with fit response values for each row in the input data.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
sigma.galamm(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Plot response versus fitted values plot(fitted(mod), response(mod))# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Plot response versus fitted values plot(fitted(mod), response(mod))
This is a very thin wrapper around mgcv::s. It enables
the specification of loading variables for smooth terms. The last letter "l",
which stands for "loading", has been added to avoid namespace conflicts with
mgcv and gamm4.
sl(..., factor = NULL)sl(..., factor = NULL)
... |
Arguments passed on to |
factor |
Optional character argument specifying the loading variable. Case sensitive. |
The documentation of the function mgcv::s should be consulted
for details on how to properly set up smooth terms. In particular, note
that these terms distinguish between ordered and unordered factor terms
in the by variable, which can be provided in ... and is
forwarded to mgcv::s.
An object of class xx.smooth.spec, where xx is a basis
identifying code given by the bs argument of s. It differs
from the smooth returned by mgcv::s in that it has an additional
attribute named "factor" which specifies any factor loading which
this smooth term should be multiplied with in order to produce the observed
outcome.
Wood SN (2003). “Thin Plate Regression Splines.” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 65(1), 95–114. ISSN 1369-7412.
Wood SN (2017). Generalized Additive Models: An Introduction with R, 2 edition. Chapman and Hall/CRC.
Other modeling functions:
galamm(),
t2()
# Linear mixed model with factor structures dat <- subset(cognition, domain == 1 & timepoint == 1) loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # Model with four thin-plate regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + sl(x, k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # Model with four cubic regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + sl(x, bs = "cr", k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" )# Linear mixed model with factor structures dat <- subset(cognition, domain == 1 & timepoint == 1) loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # Model with four thin-plate regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + sl(x, k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # Model with four cubic regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + sl(x, bs = "cr", k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" )
Extracts the square root of the dispersion parameter(s) from an object of
class galamm, returned from galamm. In the case of
conditionally Gaussian responses, this is the residual standard deviation.
When there are multiple dispersion parameters, e.g., with mixed response
type models, the square root of all of them are returned in a numeric vector.
## S3 method for class 'galamm' sigma(object, ...)## S3 method for class 'galamm' sigma(object, ...)
object |
An object of class |
... |
Optional parameters passed on to other methods. Currently not used. |
The square root of one or more dispersion parameters.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract residual standard deviation. sigma(mod) # The residual standard deviation applies to the base case. The variance # function shown in the model output shows the estimated multiplier for # various grouping levels: summary(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract residual standard deviation. sigma(mod) # The residual standard deviation applies to the base case. The variance # function shown in the model output shows the estimated multiplier for # various grouping levels: summary(mod)
Summary method for class "galamm".
## S3 method for class 'galamm' summary(object, ...)## S3 method for class 'galamm' summary(object, ...)
object |
An object of class |
... |
Further arguments passed on to other methods. Currently not used. |
A list of summary statistics of the fitted model of class
summary.galamm, containing the following elements:
AICtab a table of model fit measures, returned by
llikAIC.
call the matched call used when fitting the model.
fixef a matrix with fixed effect estimated, returned by
fixef.
gam List containing information about smooth terms in the model. If
no smooth terms are contained in the model, then it is a list of length
zero.
model a list with various elements related to the model setup and
fit. See ?galamm for details.
parameters A list object with model parameters and related
information. See ?galamm for details.
Lambda An object containing the estimated factor loadings. Returned
from factor_loadings.galamm. If there are no estimated factor
loadings, then this object is NULL.
random_effects a list containing the random effects.
See ?galamm for details.
VarCorr An object of class VarCorr.galamm, returned from
VarCorr.galamm.
weights An object containing information about estimated variance
functions, when there are heteroscedastic residuals. Otherwise the object
is NULL.
Some of the code for producing summary information has been derived
from the summary methods of mgcv (author: Simon Wood) and
lme4 (Bates et al. 2015)
(authors: Douglas M. Bates, Martin Maechler, Ben Bolker, and Steve Walker).
print.summary.galamm() for the print method and summary() for
the generic.
Other summary functions:
anova.galamm(),
plot.galamm(),
plot_smooth.galamm(),
print.galamm(),
print.summary.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) summary(mod)# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) summary(mod)
This is a very thin wrapper around mgcv::t2. It enables
the specification of loading variables for smooth terms. The last letter
"l", which stands for "loading", has been added to avoid namespace
conflicts with mgcv and gamm4.
t2l(..., factor = NULL)t2l(..., factor = NULL)
... |
Arguments passed on to |
factor |
Optional character of length one specifying the loading variable. Case sensitive. |
The documentation of the function mgcv::t2 should be consulted
for details on how to properly set up smooth terms. In particular, note
that these terms distinguish between ordered and unordered factor terms
in the by variable, which can be provided in ... and is
forwarded to mgcv::t2.
An object of class xx.smooth.spec, where xx is a basis
identifying code given by the bs argument of t2. It differs
from the smooth returned by mgcv::s in that it has an additional
attribute named "factor" which specifies any factor loading which
this smooth term should be multiplied with in order to produce the observed
outcome.
Wood SN (2003). “Thin Plate Regression Splines.” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 65(1), 95–114. ISSN 1369-7412.
Wood SN (2017). Generalized Additive Models: An Introduction with R, 2 edition. Chapman and Hall/CRC.
Other modeling functions:
galamm(),
s()
# Linear mixed model with factor structures dat <- subset(cognition, domain == 1 & timepoint == 1) loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # Model with four cubic regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + t2l(x, k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # Model with four thin-plate regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + t2l(x, bs = "tp", k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" )# Linear mixed model with factor structures dat <- subset(cognition, domain == 1 & timepoint == 1) loading_matrix <- matrix(c(1, NA, NA), ncol = 1) # Model with four cubic regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + t2l(x, k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" ) # Model with four thin-plate regression splines as basis functions mod <- galamm( formula = y ~ 0 + item + t2l(x, bs = "tp", k = 4, factor = "loading"), data = dat, load.var = "item", lambda = loading_matrix, factor = "loading" )
Extract variance and correlation components from model
## S3 method for class 'galamm' VarCorr(x, sigma = 1, ...)## S3 method for class 'galamm' VarCorr(x, sigma = 1, ...)
x |
An object of class |
sigma |
Numeric value used to multiply the standard deviations. Defaults to 1. |
... |
Other arguments passed onto other methods. Currently not used. |
An object of class c("VarCorr.galamm", "VarCorr.merMod").
print.VarCorr.galamm() for the print function.
Other details of model fit:
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm(),
vcov.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract information on variance and covariance VarCorr(mod) # Convert to data frame # (this invokes lme4's function as.data.frame.VarCorr.merMod) as.data.frame(VarCorr(mod))# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract information on variance and covariance VarCorr(mod) # Convert to data frame # (this invokes lme4's function as.data.frame.VarCorr.merMod) as.data.frame(VarCorr(mod))
Calculate variance-covariance matrix for GALAMM fit
## S3 method for class 'galamm' vcov(object, parm = "beta", ...)## S3 method for class 'galamm' vcov(object, parm = "beta", ...)
object |
Object of class |
parm |
The parameters for which the variance-covariance matrix should be calculated. Character vector with one or more of the elements "theta", "beta", "lambda", and "weights". Can also be an integer vector. When given as a character, it must be in only lowercase letters. |
... |
Further arguments passed on to other methods. Currently not used. |
A variance-covariance matrix.
confint.galamm() for the method computing confidence intervals.
See vcov() for the generic function.
Other details of model fit:
VarCorr(),
coef.galamm(),
confint.galamm(),
deviance.galamm(),
factor_loadings.galamm(),
family.galamm(),
fitted.galamm(),
fixef(),
formula.galamm(),
llikAIC(),
logLik.galamm(),
nobs.galamm(),
predict.galamm(),
print.VarCorr.galamm(),
ranef.galamm(),
residuals.galamm(),
response(),
sigma.galamm()
# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract covariance matrix for fixed regression coefficients vcov(mod, parm = "beta") # and then for weights, which gives us the variance. vcov(mod, parm = "weights")# Linear mixed model with heteroscedastic residuals mod <- galamm( formula = y ~ x + (1 | id), weights = ~ (1 | item), data = hsced ) # Extract covariance matrix for fixed regression coefficients vcov(mod, parm = "beta") # and then for weights, which gives us the variance. vcov(mod, parm = "weights")