CLME (version 1.0-1)

constrained.lme: Constrained Inference for Linear Mixed Effects Models.

Description

Implements constrained inference for linear fixed or mixed effects models using distribution-free bootstrap methodology.

Usage

constrained.lme( method="QPE" , Y , X1 , X2=NULL , U=NULL ,       
          Nks=dim(X1)[1] , Qs=dim(U)[2] , constraints=list() ,
          nsim=1000 , em.eps=sqrt(.Machine$double.eps) , em.iter=500, 
          mq.eps=sqrt(.Machine$double.eps) , mq.iter=500 , 
          verbose=c(FALSE,FALSE,FALSE) , tsf=lrt.stat , 
          tsf.ind=w.stat.ind , pav.alg=NULL, hp=FALSE , seed=NULL )

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.

nsim

optional. Number of bootstrap samples to generate.

em.eps

optional. Criterion for convergence for the EM algorithm.

em.iter

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

mq.eps

optional. Criterion for convergence for the MINQUE algorithm.

mq.iter

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

verbose

optional. Vector of 3 logicals. The first causes printing of iteration step, the second two are passed as the verbose argument to the functions minque and

tsf

Function to calculate the global 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.

seed

set the seed for the RNG.

Value

• The output of constrained.lme is an object of the class clme, which is list with elements:
• method
  the method that was used to apply order constraints (QPE or PAVA)
• theta
  estimates of $\theta$ coefficients
• ssq
  vector of estimate(s) of residual variance, $\sigma^{2}_{i}$.
• tsq
  vector of estimate(s) of random effects variance, $\tau^{2}_{i}$.
• cov.theta
  the covariance matrix of $\theta$
• ts.glb
  test statistic for the global hypothesis.
• ts.ind
  vector of test statistics for each of the constraints.
• p.value
  p-value for the global hypothesis
• p.value.ind
  vector of p-values for each of the constraints
• constraints
  list containing the constraints (A) and the contrast for the global test (B).
• est.order
  sentence describing the estimated order (or whether custom constraints were specified).

Details

If any random effects are included, the function then computes MINQUE estimates of variance parameters. After, clme.em is run to obtain the observed values, followed by resid.boot. Finally, clme.em is run on each of the bootstrap samples and a $p$-value is computed.

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

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  <- list(order='simple' , decreasing=FALSE)

args(constrained.lme)

# Note that 'nsim' has been reduced for illustration

clme.out1 <- constrained.lme( method='QPE', Y=Y, X1=X1, X2=X2, U=U,
                     constraints=const, nsim=10 , tsf=lrt.stat )

clme.out2 <- constrained.lme( method='PAVA', Y=Y, X1=X1, X2=X2, U=U,
                     constraints=const , nsim=10 )


# Let program estimate the order ('constraints' not specified)
# Takes longer to run
clme.out3 <- constrained.lme( method='PAVA', Y=Y, X1=X1, X2=X2, U=U,
                     nsim=10 , em.eps=0.0001)