smsn.lmm: ML estimation of scale mixture of skew-normal linear mixed models

Description

It fits a scale mixture of skew-normal linear mixed model with possible within-subject dependence structure, using an EM-type algorithm (via the DAAREM method, by default). It provides estimates and standard errors of parameters.

Usage

smsn.lmm(data, formFixed, groupVar, formRandom = ~1, depStruct = "UNC",
         timeVar = NULL, distr = "sn", covRandom = "pdSymm",
         skewind, pAR = 1, control = lmmControl())

Value

An object of class "SMSN" representing the SMSN-LMM fit. Generic functions such as print and summary have methods to show the results of the fit. The functions fitted and ranef can be used to extract some of its components.

Specifically, the following components are returned:

theta: Named vector with parameter estimates.
iter: Number of iterations runned.
estimates: A named list containing parameter estimates.
uhat: Estimated weights.
loglik.track: Vector containing the log-likelihood value from each iteration of the estimation procedure.
random.effects: Estimated random effects.
std.error: A vector with standard errors.
loglik: Value of the log-likelihood at last iteration.
elapsedTime: Time elapsed in processing, in seconds.
error: Convergence criterion at last iteration.
call: The smsn.lmm call that produced the object.
criteria: A list with AIC and BIC criterion.
data: The data frame used on smsn.lmm call.
formula: A list containing the formulas used on smsn.lmm call.
depStruct: A character indicating which dependence structure was used.
covRandom: A character indicating which structure was used for the random effects scale matrix.
distr: A character indicating which distribution was used.
N: The number of observations used.
n: The number of individuals/groups used.
groupVar: A character indicating the name of the grouping variable.
control: The control list used in the function's call.
timeVar: A character indicating the name of the time variable, if any.
fitted: A vector of fitted values, if calc.bi=TRUE.

Arguments

data: A data frame containing the variables named in formFixed, formRandom, groupVar, and timeVar.
formFixed: A two-sided linear formula object describing the fixed effects part of the model, with the response on the left of a ~ operator and the covariates, separated by + operators, on the right.
groupVar: A character containing the name of the variable which represents the subjects or groups in data.
formRandom: A one-sided linear formula object describing the random effects part of the model, with the covariates, separated by + operators, on the right of a ~ operator. By default, a model with random intercept is considered.
depStruct: A character indicating which dependence structure should be used. "UNC" for conditionally uncorrelated ("CI" is also accepted), "ARp" for AR(p) -- p is length(phi)--, "CS" for compound symmetry, "DEC" for DEC, and "CAR1" for continuous-time AR(1).
timeVar: A character containing the name of the variable which represents the time in data. Meaningless if depStruct="UNC" or depStruct="CS". For other structures, if is.null(timeVar) the observations are considered equally spaced and ordered. If depStruct="ARp", timeVar must be an index, preferably starting at 1.
distr: A character indicating which distribution should be used. "sn" for skew-normal, "st" for skew-t, "ssl" for skew-slash, and "scn" for skew-contaminated normal.
covRandom: A character indicating which structure should be used for the random effects scale matrix (either "pdSymm" (default), for a general positive-definite matrix, or "pdDiag", for a diagonal matrix).
skewind: A vector of length equal to the number of random effects, containing 0's and 1's, indicating which elements of the skewness parameter vector should be estimated (optional, default is an all-ones vector).
pAR: If depStruct="ARp", pAR indicates the order of the autoregressive process that should be used. Otherwise, it is meaningless.
control: An object resulting from the function lmmControl(), containing additional options for the estimation algorithm.

Author

Fernanda L. Schumacher, Larissa A. Matos and Victor H. Lachos

Details

It fits the model \(Y_i = X_i \beta + Z_i b_i + \epsilon_i\), for \(i=1,\dots,n\), where \(Y_i\) is a vector with \(n_i\) observed continuous responses, \(b_i ~ SMSN(c\Delta, D,\lambda;H)\) and \(\epsilon_i ~ SMN(0, \Sigma_i; H)\), indexed by the same mixing distribution.

For efficiency, the DAAREM method is used for parameter estimation. In case of numerical errors, please try passing to lmmControl(algorithm = "DAAREM") to the control argument.

For details see Schumacher, Lachos & Matos (2021).

References

Henderson, N. C. and R. Varadhan (2019). Damped anderson acceleration with restarts and monotonicity control for accelerating EM and EM-like algorithms. Journal of Computational and Graphical Statistics 28(4), 834-846.

Lachos, V. H., P. Ghosh, and R. B. Arellano-Valle (2010). Likelihood based inference for skew-normal independent linear mixed models. Statistica Sinica 20, 303-322.

Schumacher, F. L., Matos, L. A., and Lachos, V. H. (2025). "skewlmm: An R Package for Fitting Skewed and Heavy-Tailed Linear Mixed Models." Journal of Statistical Software, 115(7), 1--32.

Schumacher, F. L., Lachos, V. H., and Matos, L. A. (2021). Scale mixture of skew-normal linear mixed models with within-subject serial dependence. Statistics in Medicine 40(7), 1790-1810.

Examples

Run this code

#simple example
dat1 = as.data.frame(nlme::Orthodont)
fm1 = smsn.lmm(dat1, formFixed=distance ~ age, groupVar="Subject",
                control=lmmControl(max.iter=30))
fm1
# \donttest{
#fitting for several distributions / dependence structures
fm1 = smsn.lmm(dat1, formFixed=distance ~ age+Sex, groupVar="Subject")
fm2 = smsn.lmm(dat1, formFixed=distance ~ age+Sex, groupVar="Subject",
               distr="st")
fm3 = smsn.lmm(dat1, formFixed=distance ~ age+Sex, groupVar="Subject",
               distr="ssl")
fm4 = smsn.lmm(dat1, formFixed=distance ~ age+Sex, groupVar="Subject",
               depStruct="ARp", pAR=1)
criteria(list(fm1=fm1, fm2=fm2, fm3=fm3, fm4=fm4))
summary(fm3)

#some diagnostic tools
plot(fm3)
acf3 = acfresid(fm3, calcCI=TRUE, MCiter=100)
plot(acf3)
plot(mahalDist(fm3), nlabels=2)
healy.plot(fm3)
# }

Run the code above in your browser using DataLab