lmenssp.heavy: Function to obtain the maximum likelihood estimates of the...
In lmenssp: Linear Mixed Effects Models with Non-Stationary Stochastic Processes

Description Usage Arguments Details Value Author(s) References Examples

Obtains the maximum likelihood estimates of the parameters by expectatio-maximisation (E-M) algorithm for linear mixed effects models with random intercept and a stationary or non-stationary stochastic process component, under multivariate normal t distribution

lmenssp.heavy(formula, data, id, timeVar, init.em = NULL, maxiter.em = 100, 
tol.em = 0.001, process, silent = TRUE, dof.est = c(0.1, 10, 1e-04), 
tol.cd = 0.001, tol.lmenssp = 10^-5, init.lmenssp = NULL, maxiter.lmenssp = 100, 
silent.lmenssp = TRUE)

`formula`	a typical `R` formula for the fixed effects component of the model
`data`	a data frame from which the variables are to be extracted
`id`	a numerical vector for subject identification
`timeVar`	a numerical vector for the time variable
`init.em`	a vector of inital values for the E-M algorithm
`maxiter.em`	a numerical value for the maximum number of iterations for the E-M algorithm
`tol.em`	a numerical value for the maximum tolerance to assess the convergence of the E-M algorithm
`process`	a character string for the stochastic process: `"bm"` for Brownian motion, `"ibm"` for integrated Brownian motion, `"iou"` for integrated Ornstein-Uhlenbeck process, `"sgp-powered-power-method"` for stationary process with powered correlation function and `"sgp-matern-kappa"` for stationary process with Matern correlation function
`silent`	a character string, if set to `FALSE` the details of the E-M steps are printed when the algorithm is running
`dof.est`	a vector of three elements, to be passed to `optimize`
`tol.cd`	a numerical value for the tolerance of central-difference approximation
`tol.lmenssp`	a numerical value for the maximum tolerance to assess the convergence, to be passed to `lmenssp`
`init.lmenssp`	a vector of inital values, to be passed to `lmenssp`
`maxiter.lmenssp`	a numerical value of the number of iterations for the Fisher-Scoring or Nelder-Mead algorithms, to be passed to `lmenssp`
`silent.lmenssp`	a character string, if set to `FALSE` the details of the Fisher-Scoring steps are printed when the algorithm is running, to be passed to `lmenssp`

lmenssp.heavy calls lmenssp inside.

Whilst theoretical standard errors are calculated and reported only for the fixed effects, central-difference based standard errors are calculated and reported for all the parameter estimates.

There are more than one way of specifying init.em, it can be set to:

1) fixed effects, random effects parameters and degrees-of freedom,

2) only the degrees-of-freedom, and

3) NULL; for this specification, lmenssp.heavy finds the intials internally.

For the details of "process", see lmenssp.

In dof.est, first and second elements are the minimum and maximum values of the search and the third element is the tolerance. It is passed to optimize.

Returns a list of results.

Ozgur Asar, Peter J. Diggle

Asar O, Ritchie J, Kalra P, Diggle PJ (2015) Acute kidney injury amongst chronic kidney disease patients: a case-study in statistical modelling. To be submitted.

Pinheiro JC, Liu C, Wu YN. (2001) Efficient algorithms for robust estimation in linear mixed-effects models using the multivariate t distribution. Journal of Computational and Graphical Statistics 10, 249-276.

# loading the data set and subsetting it for the first 20 patients 
# for the sake illustration of the usage of the functions
data(data.sim.ibm.heavy)
data.sim.ibm.heavy.short <- data.sim.ibm.heavy[data.sim.ibm.heavy$id <= 20, ]

# estimating the parameters
# tol.em is set to 10^-1 and tol.lmenssp to 10^-2 only for illustration, 
# decrease these values in your applications
fit.heavy <- lmenssp.heavy(formula = log.egfr ~ sex + bage + fu + pwl, 
       data = data.sim.ibm.heavy.short, id = data.sim.ibm.heavy.short$id, 
       timeVar = data.sim.ibm.heavy.short$fu, init.em = 5, maxiter.em = 1000, 
       tol.em = 10^-1, process = "ibm", silent = FALSE, 
       dof.est = c(0.1, 10, 0.0001), tol.cd = 0.001, tol.lmenssp = 10^-2, 
       silent.lmenssp = FALSE)
fit.heavy