smsn.lmm {skewlmm}

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())

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.

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).

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.

Author(s)

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

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., 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.

See Also

lmmControl, update, predict.SMSN, residuals.SMSN, plot.SMSN, smn.lmm, smn.clmm

Examples

#simple example
dat1 = as.data.frame(nlme::Orthodont)
fm1 = smsn.lmm(dat1, formFixed=distance ~ age, groupVar="Subject",
                control=lmmControl(max.iter=30))
fm1

#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)