R: ML estimation of scale mixture of normal linear mixed models...

smn.clmm {skewlmm}

R Documentation

ML estimation of scale mixture of normal linear mixed models with censored responses

Description

It fits left, right, or interval censored scale mixture of normal linear mixed model with possible within-subject dependence structure, using the EM algorithm. It provides estimates and standard errors of parameters.

Usage

smn.clmm(data, formFixed, groupVar, formRandom = ~1, depStruct = "UNC",
         ci, lcl, ucl, timeVar = NULL, distr = "norm",
         nufix = FALSE, pAR = 1, control = lmmControl())

Arguments

`data`	A data frame containing the variables named in `formFixed`, `formRandom`, `groupVar`, `timeVar`, `ci`, `lcl`, and `ucl`.
`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, `"CAR1"` for continuous-time AR(1), and `"MA1"` for moving average of order 1.
`ci`	A character containing the name of the censoring indicator variable in `data`, which should be 1 if the respective observation is censored or missing, and 1 otherwise. If missing, it is assumed that none of the observations is censored.
`lcl`	A character containing the name of the lower censoring limit in `data`. If missing, it is assumed `lcl=-Inf`, i.e., no left limit.
`ucl`	A character containing the name of the upper censoring limit in `data`. If missing, it is assumed `ucl=Inf`, i.e., no right limit.
`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. `"norm"` for normal and `"t"` for Student-t.
`nufix`	`TRUE` or `FALSE` indicating if `nu` should be estimated for t distribution. If `nufix=TRUE`, nu must be especified through `lmmControl()`.
`pAR`	If `depStruct="ARp"`, `pAR` indicates the order of the autoregressive process that should be used (1 by default). 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 ~ SMN(0, D; H) and \epsilon_i ~ SMN(0, \Sigma_i; H), indexed by the same mixing distribution.

For details see Matos et al. (2013) and Lachos et al. (2019).

Value

An object of class "SMNclmm" representing the SMN-CLMM 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.
`yest`	Imputed values in the response variable.
`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.
`criteria`	A list with AIC, BIC, and SIC criterion.
`call`	The `smn.clmm` call that produced the object.
`data`	The data frame used on `smn.clmm` call.
`formula`	A list containing the formulas used on `smn.clmm` 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.
`timeVar`	A character indicating the name of the time variable, if any.
`fitted`	A vector of fitted values.

Author(s)

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

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.

Matos, L. A., Prates, M. O., Chen, M. H., and Lachos, V. H. (2013). Likelihood-based inference for mixed-effects models with censored response using the multivariate-t distribution. Statistica Sinica 23(3), 1323-1345.

Lachos, V. H., A. Matos, L., Castro, L. M., and Chen, M. H. (2019). Flexible longitudinal linear mixed models for multiple censored responses data. Statistics in medicine, 38(6), 1074-1102.

Examples

# Generating a sample for m=30 individuals with 5 times
# Considering 10% of left censoring responses
# AR(1) and normal distribution
nj1 = 5
m = 30
time = rep(1:nj1, times=m)
groups = as.factor(rep(1:m, each=nj1))
dat1 = rsmsn.clmm(time = time, ind = groups, x = cbind(1,time), z = rep(1,m*nj1),
                  sigma2=0.6, D=0.5*diag(1), beta=c(1,2), depStruct="ARp",
                  phi=0.4, pcens = .1, type = "left")

# Estimation using UNC
fm0 = smn.clmm(dat1, formFixed=y~x, groupVar="ind", ci="ci", ucl="ucl")
summary(fm0)

# Estimation using AR(1)
fm1 = update(fm0, depStruct="ARp")

# Estimation using AR(1) and t distribution
fm2 = update(fm1, distr="t")

# Comparing fitted models
criteria(list(fm0= fm0, fm1=fm1, fm2=fm2))

[Package skewlmm version 1.1.0 Index]