# R/pass.lme.R In pass.lme: Power and Sample Size for Linear Mixed Effect Models

#### Documented in lme.Lb.dist.thetapass.lme.CLb.test

```.onAttach <- function(libname,pkgname) {
packageStartupMessage('Title: Power and Sample Size for Linear Mixed Effect Models')
packageStartupMessage('Author: Marco Chak Yan YU')
packageStartupMessage('Maintainer: Marco Chak Yan YU <[email protected]>')
}

##################################################
#' @title Calculate mean and variance  \cr
#' for linear combination of the  \cr
#' Best Linear Unbiased Estimator (BLUE) \cr
#' for Linear Mixed Effect (LME) Model
#' @author Marco Chak Yan YU \cr
#' Maintainer: Marco Chak Yan YU <[email protected]>
#' @description
#' Consider the following model: \cr \cr
#'       Y=XB+Zu+e, u~N(0,D), e~N(0,R) \cr
#'       Yi~N(XBi,Vi), Vi=Zi*D*Zi'+Ri, \cr \cr
#'           for each independent observation i \cr \cr \cr
#' estimate of fixed effect coefficient B, denoted by b: \cr \cr
#'       b=inv(sum(Xi'*inv(Vi)*Xi))*(sum(Xi'*inv(Vi)*Yi)) \cr \cr \cr
#' variance of b: \cr \cr
#'       var(b)=Vb/n=inv(sum(Xi'*inv(Vi)*Xi)) \cr \cr
#'           where Vb=inv(Xi'*inv(Vi)*Xi) \cr \cr \cr
#' distribution of any linear combinations L of b is given by: \cr \cr
#'       Lb~N(mu,Sigma/n) \cr \cr
#'             where mu = LB, Sigma = L*Vb*L' \cr \cr \cr
#' @param B      fixed effect beta in px1 matrix
#' @param D      list of qxq random effect variance matrix;
#'                  where the first element corresponding to the highest-level effect,
#'                  the last element corresponding to the level 1 effect \cr \cr
#' @param R      residual variance \cr \cr
#' @param L      lxp matrix, representing l-linear-combinations of beta interested, \cr
#'                  if L is not defined, it will be auto-created to select the last coefficient \cr \cr
#' @param X      nxp matrix representing the covariates for the fixed effects \cr \cr
#' @param Z      nxq matrix representing the covariates for each level of random effects \cr \cr
#' @param m      vector of repeated measures from the highest to lowest level
#'                  (level 1 effects are addressed by Z and X and no need to be specified) \cr \cr
#' @return
#' theta:  parameters (mu and Sigma) of the normal distribution for Lb
#' @examples
#' #Example 1
#' # calc BLUE for 1-level LME model,
#' # with covariates X, Z: (1,t), t=1,2,3
#' # for both fixed and random effects,
#' # with fixed effect coefficients B: (100,-0.5),
#' # random effect variance D: (2 1;1 2),
#' # residual variance R: 0.2
#' B <- matrix(c(100,-0.5),2,1)
#' D <- matrix(c(2,1,1,2),2,2)
#' R <- 0.2
#' X <- cbind(rep(1,3),1:3)
#' Z <- X
#' lme.Lb.dist.theta(B,D,R,X,Z)
#'
#' #Example 2
#' # calc BLUE for 3-levels LME model,
#' # with level 1 same as the above example
#' # with 3 repeated-measures in level 2
#' # and 5 repeated-measures in the highest level,
#' # assuming random effect variance for level 2 = (3 1;1 3),
#' # and random effect variance for highest level = (5 1;1 5)
#' D <- list(matrix(c(2,1,1,2),2,2),matrix(c(3,1,1,3),2,2), matrix(c(5,1,1,5),2,2))
#' lme.Lb.dist.theta(B,D,R,X,Z,m=c(5,3))
##
lme.Lb.dist.theta <- function(B,D,R,X,Z,m=NULL,L) {
##########
# format input
B <- matrix(c(B),ncol=1)
if (missing(L)) {
L <- matrix(c(rep(0,nrow(B)-1),1),nrow=1)
} else if (is.vector(L)) {
L <- matrix(c(L),nrow=1)
} else if (dim(L)!=dim(B)) {
L <- t(L)
}
if (!is.list(D)) D <- list(as.matrix(D))
##########
# input check
if (dim(B)!=dim(X)) {
stop('inconsistent matrix dimension between X and B, requires dim(X) == dim(B)')
} else if (dim(L)!=dim(B)) {
stop('inconsistent matrix dimension between L and B, requires dim(L) == dim(B)')
} else if (length(m)!=length(D)-1) {
stop('inconsistent input length between D and m, requires length(D)-1 == length(m)')
} else if (any(sapply(
1:length(D),function(i)
(dim(D[[i]])!=dim(D[[i]]))&
(dim(D[[i]])!=dim(Z))))) {
stop('inconsistent matrix dimension between D and Z, requires dim(Z) == dim(D[[i]]) == dim(D[[i]]) for all i')
} else if (dim(X)!=dim(Z)) {
stop('inconsistent matrix dimension between X and Z, requires dim(X) == dim(Z)')
} else if (length(c(R))!=1) {
stop('wrong input R, requires a single number')
} else {
##########
# perform the calculation
LB <- L%*%B
# calc Zi*D*Zi' in Vi=Zi*D*Zi'+Ri
Zi <- Z
Xi <- X
Vi <- Zi %*% matrix(D[[length(D)]],dim(D[[length(D)]])) %*% t(Zi)
if (length(m)>=1) {
for (i in seq(length(m),1,-1)) {
Zi <- kronecker(matrix(rep(1,m[i]),ncol=1),Zi)
Xi <- kronecker(matrix(rep(1,m[i]),ncol=1),Xi)
Vi <- kronecker(diag(rep(1,m[i])),Vi)
Vi <- Vi + Zi %*% matrix(D[[i]],dim(D[[i]])) %*% t(Zi)
}
}
# calc Vi=Zi*D*Zi'+Ri
Vi <- Vi + diag(R,dim(Vi))
# Vb/n=var(b)=inv(sum(Xi'*inv(Vi)*Xi))
Vb <- solve(t(Xi) %*% solve(Vi) %*% Xi)
# Lb~N(LB,L*Vb*L'/n)
VLb <- L %*% Vb %*% t(L)
# generate output
theta <- NULL
theta\$mu <- LB
theta\$Sigma <- VLb
return(theta)
}
}

##################################################
#' @title Calculate Power or Sample Size required \cr
#' for Contrasts of linear combinations \cr
#' of fixed effect parameters \cr
#' in Linear Mixed Effect (LME) Model
#' @author Marco Chak Yan YU \cr
#' Maintainer: Marco Chak Yan YU <[email protected]>
#' @description
#' Interested parameters/linear combinations LB from more than one
#' independent populations can be aggregrate togeter by
#' appending mu vertically and Sigma/n diagonally \cr \cr
#' Consider Lb~N(MU,SIGMA) as the aggregrated estimates\cr
#' Any comparison of interested parameters can be formulated by
#' multiplying a contrast matrix C on LB and set \cr \cr
#' H0: C*LB=d for any vector of value d to be tested \cr \cr
#' We then have \cr \cr
#'       C*Lb~N(C*MU,C*SIGMA*C') \cr \cr
#' and \cr \cr
#'       (C*Lb-d)'*inv(C*SIGMA*C')*(C*Lb-d)~chisq(q,lambda) \cr \cr
#' where degree of freedom q=rank(C*SIGMA*C'), \cr
#' non-centrality parameter lambda=(C*LB-d)'*inv(C*SIGMA*C')*(C*LB-d) \cr \cr \cr
#' Power of the test H0 is given by 1-beta=P(chisq(q,lambda)>qchisq(1-alpha,lambda)) \cr
#' Required sample size for desired power can be obtained by bisection method.
#' @param thetas  list of theta (LB and VLb), can be different for each group \cr \cr
#' @param C       Contrast of Matrix \cr \cr
#' @param d       Value vector to be tested for all contrast \cr \cr
#' @param alpha   significant level \cr \cr
#' @param power   desired power for sample size calculation \cr \cr
#' @param n       sample size for power calculation / \cr
#'                    or sample size ratio with power for sample size calculation
#'                       (NULL for balanced design) \cr \cr
#' @return
#'   solved.power  given sample size n, this gives the power for testing H0 \cr
#'   solved.n      given the desired power, this gives the sample size for H0 \cr
#' @examples
#' #Example 1 (test fixed effect coefficient 2=0) with power of 80%
#' # for 1-level LME model, with covariates X, Z: (1,t), t=1,2,3
#' # for both fixed and random effects, with fixed effect coefficients B: (100,-0.5),
#' # random effect variance D: (2 1;1 2), residual variance R: 0.2
#' B <- matrix(c(100,-0.5),2,1)
#' D <- matrix(c(2,1,1,2),2,2)
#' R <- 0.2
#' X <- cbind(rep(1,3),1:3)
#' Z <- X
#' theta <- lme.Lb.dist.theta(B,D,R,X,Z)
#' pass.lme.CLb.test(list(theta),alpha=0.05,power=0.8)
#' pass.lme.CLb.test(list(theta),alpha=0.05,n=66)
#'
#' #Example 2 (compare two fixed effect coefficient 2) with power of 80%
#' # Consider above model as a control group model,
#' # with an independent treatment group with model same as the control
#' # except a different fixed effect coefficient 2 for treatment
#' # = fixed effect coefficient 2 for control x 0.7
#' theta2 <- theta
#' theta2\$mu <- theta\$mu *0.7
#' C <- matrix(c(1,-1),1,2)
#' pass.lme.CLb.test(list(theta,theta2),C,alpha=0.05,power=0.8)
#' pass.lme.CLb.test(list(theta,theta2),C,alpha=0.05,n=1468)
#'
#' #Example 3 (compare two fixed effect coefficient 2) with power of 80%
#' # with sample size ratio, control:treatment = 1:2
#' pass.lme.CLb.test(list(theta,theta2),C,alpha=0.05,power=0.8,n=c(1,2))
#' pass.lme.CLb.test(list(theta,theta2),C,alpha=0.05,n=c(1101,2202))
#'
#' #Example 4 (repeated-measures ANOVA for comparing 3 group means) with power of 80%
#' # for 1-level LME model with mean for group 1, 2 and 3 are 100, 99, 102, respectively,
#' # each subject to be measured 2 times, with within-subject variance = 15, residual variance = 10
#' B <- 100
#' D <- 15
#' R <- 10
#' X <- matrix(1,2,1)
#' Z <- X
#' theta <- lme.Lb.dist.theta(B,D,R,X,Z)
#' theta2 <- theta
#' theta3 <- theta
#' theta2\$mu <- 99
#' theta3\$mu <- 102
#' C <- rbind(c(1,-1,0),c(1,0,-1))
#' pass.lme.CLb.test(list(theta,theta2,theta3),C,alpha=0.05,power=0.8)
#' pass.lme.CLb.test(list(theta,theta2,theta3),C,alpha=0.05,n=41)
##
pass.lme.CLb.test <- function(thetas,C=NULL,d=NULL,alpha=0.05,power=NULL,n=NULL) {
if (is.null(n)) {
n <- rep(1,length(thetas))
} else if ((length(n)==1)&(length(thetas)>1)) {
n <- rep(n,length(thetas))
}
LB <- NULL
for (i in 1:length(thetas)) {
LB <- rbind(LB,thetas[[i]]\$mu)
}
VLb <- matrix(0,dim(LB),dim(LB))
j <- 0
for (i in 1:length(thetas)) {
VLbi <- thetas[[i]]\$Sigma
VLb[(j+1):(j+dim(VLbi)),(j+1):(j+dim(VLbi))] <- VLbi/n[i]
j <- j+dim(VLbi)
}

if (is.null(C)) C <- diag(dim(LB))
if (is.null(d)) d <- matrix(0,nrow(C),1)
VCLb <- C%*%VLb%*%t(C)
df <- qr(VCLb)\$rank
chisq.ncp <- t(C%*%LB-d) %*% solve(VCLb) %*% (C%*%LB-d)

if (is.null(power)) {
solved.power <- 1-pchisq(qchisq(1-alpha,df,0),df,chisq.ncp)
return(solved.power)
} else {
solved.n <- 100
solved.n.lb <- 0
solved <- FALSE
while (!solved) {
if ((1-pchisq(qchisq(1-alpha,df,0),df,chisq.ncp*solved.n) < power)) {
solved.n.lb <- solved.n
solved.n <- solved.n*10
} else if ((1-pchisq(qchisq(1-alpha,df,0),df,chisq.ncp*solved.n) > power+0.001)) {
solved.n <- (solved.n + solved.n.lb)/2
} else {
solved <- TRUE
}
}
solved.n <- solved.n *n
return(solved.n)
}
}
```

## Try the pass.lme package in your browser

Any scripts or data that you put into this service are public.

pass.lme documentation built on Aug. 20, 2019, 5:13 p.m.