CLME (version 1.0-1)

clme.em: Constrained Expectation-Maximization Algorithm

Description

Implements an expectation-maximization (EM) algorithm under inequality constraints to estimate parameters and compute a test statistic.

Usage

clme.em.all(method, Y, X1, X2 = NULL, U = NULL, Nks = dim(X1)[1],
        Qs = dim(U)[2] , constraints, mq.phi = NULL,
        tsf=lrt.stat , tsf.ind=w.stat.ind , pav.alg, hp=FALSE , em.iter=500,
        em.eps = sqrt(.Machine$double.eps), verbose = FALSE)
clme.em.fixed(method, Y, X1, X2 = NULL, U = NULL, Nks = dim(X1)[1],
        Qs = dim(U)[2] , constraints, mq.phi = NULL,
        tsf=lrt.stat , tsf.ind=w.stat.ind , pav.alg, hp=FALSE , em.iter=500,
        em.eps = sqrt(.Machine$double.eps), verbose = FALSE)
clme.em.mixed(method, Y, X1, X2 = NULL, U = NULL, Nks = dim(X1)[1],
        Qs = dim(U)[2] , constraints, mq.phi = NULL,
        tsf=lrt.stat , tsf.ind=w.stat.ind , pav.alg, hp=FALSE , em.iter=500,
        em.eps = sqrt(.Machine$double.eps), verbose = FALSE)

Arguments

method

method for applying order constraints. Supported are QPE and PAVA

Y

$N \times 1$ vector of response data.

X1

$N \times p_1$ design matrix.

X2

optional. $N \times p_2$ matrix of covariates.

U

optional. $N \times c$ matrix of random effects.

Nks

optional. $K \times 1$ vector of group sizes.

Qs

optional. $Q \times 1$ vector of group sizes for random effects.

constraints

optional. List containing the constraints, defaults to empty list. See Details for further information.

mq.phi

optional. MINQUE estimates of variance parameters.

tsf

Function to calculate the test statistic. Defaults to likelihood ratio test statistic (lrt.stat). See Details for further information.

tsf.ind

Function to calculate the test statistic for individual constrats. Defaults to Williams type test statistic (w.stat.ind). See Details for further information.

pav.alg

algorithm to apply PAVA constraints. Only required when custom constraints are specified and method=PAVA. Ignored when method=QPE. See Details for further information.

hp

logical to indicate whether weights for PAVA should be use the full covariance matrix or the diagonal elements. Ignored when method=QPE.

em.eps

optional. Criterion for convergence for the EM algorithm.

em.iter

optional. Maximum number of iterations permitted for the EM algorithm.

verbose

If TRUE, function prints messages on progress of the EM algorithm.

Value

• The function returns a list with the elements:
• theta
  vector of theta.
• ssq
  vector of estimates of variance terms.
• tsq
  vector of estimates of variance terms for the random effects.
• cov.theta
  the covariance matrix of the unconstrained theta coefficients.
• ts.glb
  test statistic for the global hypothesis.
• ts.ind
  vector of test statistics for each of the constraints.

Details

Constraints can be implemented either through Restricted Maximum Likelihood Estimation (by setting method=QPE) or the Pool Adjacent Violators Algorithm (method=PAVA). For PAVA, the user must also specify an appropriate algorithm. Several default PAVA algorithms are included in the package CLME. See pava.functions for more details. Argument constraints is a list including at least the elements A and B. This argument can be generated by function create.constraints.

References

Farnan, L., Ivanova, A., and Peddada, S. D. (2014). Linear Mixed Efects Models under Inequality Constraints with Applications. PLOS ONE, 9(1). e84778. doi: 10.1371/journal.pone.0084778 http://www.plosone.org/article/info:doi/10.1371/journal.pone.0084778

See Also

CLME-package, constrained.lme, pava.functions, create.constraints, w.stat

Examples

set.seed( 42 )

n  <- 5
P1 <- 5

X1 <- diag(P1) %x% rep(1,n)
X2 <- as.matrix( rep(1,P1) %x% runif(n , 0,2) )
U  <- rep(1,P1) %x% diag(n)
X  <- as.matrix( cbind(X1,X2) )

tsq <- 1
ssq <- 0.7

Nks <- dim(X1)[1]
Qs  <- dim(U)[2]

xi <- rnorm( sum(Qs)  , 0 , rep(sqrt(tsq) , Qs)  )
ep <- rnorm( sum(Nks) , 0 , rep(sqrt(ssq) , Nks) )  

thetas <- c(2 , 3 , 3, 3 , 4 , 2 )
Y      <- X%*%thetas + U%*%xi + ep
const  <- create.constraints( X1=X1 , X2=X2 , 
          constraints=list(order='simple' , decreasing=FALSE) )

## Two runs of clme.em with different options

clme.out1 <- clme.em.all( method='QPE', Y=Y, X1=X1, X2=X2, U=U, 
                     constraints=const )

clme.out2 <- clme.em.mixed( method='QPE', Y=Y, X1=X1, X2=X2, U=U, 
                     constraints=const )
                    
                    
clme.out3 <- clme.em.all( method='PAVA', Y=Y, X1=X1, X2=X2, U=U,
                     constraints=const , tsf=w.stat ,
                    pav.alg=pava.simple.order )