liu.liang.linear.power {longpower}

Linear mixed model sample size calculations from Liu & Liang (1997).

Description

This function performs the sample size calculation for a linear mixed model. See Liu and Liang (1997) for parameter definitions and other details.

Usage

liu.liang.linear.power(
  N = NULL,
  delta = NULL,
  u = NULL,
  v = NULL,
  sigma2 = 1,
  R = NULL,
  R.list = NULL,
  sig.level = 0.05,
  power = NULL,
  Pi = rep(1/length(u), length(u)),
  alternative = c("two.sided", "one.sided"),
  tol = .Machine$double.eps^2
)

Arguments

N

The total sample size. This formula can accommodate unbalanced group allocation via Pi. See Liu and Liang (1997) for more details

delta

group difference (possibly a vector of differences)

u

a list of covariate vectors or matrices associated with the parameter of interest

v

a respective list of covariate vectors or matrices associated with the nuisance parameter

sigma2

the error variance

R

the variance-covariance matrix for the repeated measures

R.list

a list of variance-covariance matrices for the repeated measures, if assumed different in two groups

sig.level

type one error

power

power

Pi

the proportion of covariates of each type

alternative

one- or two-sided test

tol

numerical tolerance used in root finding.

Details

The parameters u, v, and Pi are expected to be the same length and sorted with respect to each other. See Liu and Liang (1997) and package vignette for more details.

References

Liu, G. and Liang, K. Y. (1997) Sample size calculations for studies with correlated observations. Biometrics, 53(3), 937-47.

See Also

lmmpower

Examples

## Not run: 
browseVignettes(package = "longpower")

## End(Not run)

# Reproduces the table on page 29 of Diggle et al for
# difference in slopes between groups

n <- 3
t <- c(0,2,5)
u <- list(u1 = t, u2 = rep(0,n))
v <- list(v1 = cbind(1,1,t),
         v2 = cbind(1,0,t))         
rho <- c(0.2, 0.5, 0.8)
sigma2 <- c(100, 200, 300)
tab <- outer(rho, sigma2, 
      Vectorize(function(rho, sigma2){
        ceiling(liu.liang.linear.power(
          delta=0.5, u=u, v=v,
          sigma2=sigma2,
          R=rho, alternative="one.sided",
          power=0.80)$N/2)}))
colnames(tab) <- paste("sigma2 =", sigma2)
rownames(tab) <- paste("rho =", rho)
tab

# Reproduces the table on page 30 of Diggle et al for 
# difference in average response between groups.

n <- 3
u <- list(u1 = rep(1,n), u2 = rep(0,n))
v <- list(v1 = rep(1,n),
         v2 = rep(1,n))
rho <- c(0.2, 0.5, 0.8)
delta <- c(20, 30, 40, 50)/100
tab <- outer(rho, delta, 
     Vectorize(function(rho, delta){
       ceiling(liu.liang.linear.power(
         delta=delta, u=u, v=v,
         sigma2=1,
         R=rho, alternative="one.sided",
         power=0.80)$n[1])}))
colnames(tab) <- paste("delta =", delta)
rownames(tab) <- paste("rho =", rho)
tab

# An Alzheimer's Disease example using ADAS-cog pilot estimates
# var of random intercept
sig2.i <- 55
# var of random slope
sig2.s <- 24
# residual var
sig2.e <- 10
# covariance of slope and intercep
cov.s.i <- 0.8*sqrt(sig2.i)*sqrt(sig2.s)

cov.t <- function(t1, t2, sig2.i, sig2.s, cov.s.i){
        sig2.i + t1*t2*sig2.s + (t1+t2)*cov.s.i 
}

t <- seq(0,1.5,0.25)
n <- length(t)
R <- outer(t, t, function(x,y){cov.t(x,y, sig2.i, sig2.s, cov.s.i)})
R <- R + diag(sig2.e, n, n)
u <- list(u1 = t, u2 = rep(0,n))
v <- list(v1 = cbind(1,1,t),
         v2 = cbind(1,0,t))         

liu.liang.linear.power(delta=1.5, u=u, v=v, R=R, sig.level=0.05, power=0.80)
liu.liang.linear.power(N=416, u=u, v=v, R=R, sig.level=0.05, power=0.80)
liu.liang.linear.power(N=416, delta = 1.5, u=u, v=v, R=R, sig.level=0.05)
liu.liang.linear.power(N=416, delta = 1.5, u=u, v=v, R=R, power=0.80, sig.level = NULL)

# Reproduces total sample sizes, m, of Table 1 of Liu and Liang 1997
tab1 <- data.frame(cbind(
  n = c(rep(4, 4), rep(2, 4), 1),
  rho = c(0.0, 0.3, 0.5, 0.8)))
m <- c()
for(i in 1:nrow(tab1)){
  R <- matrix(tab1$rho[i], nrow = tab1$n[i], ncol = tab1$n[i])
  diag(R) <- 1
  m <- c(m, ceiling(liu.liang.linear.power(
    delta=0.5,
    u = list(u1 = rep(1, tab1$n[i]), # treatment
             u2 = rep(0, tab1$n[i])), # control       
    v = list(v1 = rep(1, tab1$n[i]), v2 = rep(1, tab1$n[i])), # intercept
    sigma2=1,
    R=R, alternative="two.sided",
    power=0.90)$N))
}
cbind(tab1, m)

# Reproduces total sample sizes, m, of Table 3.a. of Liu and Liang 1997
# with unbalanced design
tab3 <- data.frame(cbind(
  rho = rep(c(0.0, 0.3, 0.5, 0.8), 2),
  pi1 = c(rep(0.8, 4), rep(0.2, 4))))
m <- c()
for(i in 1:nrow(tab3)){
  R <- matrix(tab3$rho[i], nrow = 4, ncol = 4)
  diag(R) <- 1
  m <- c(m, ceiling(liu.liang.linear.power(
    delta=0.5,
    u = list(u1 = rep(1, 4), # treatment
             u2 = rep(0, 4)), # control       
    v = list(v1 = rep(1, 4), v2 = rep(1, 4)), # intercept
    sigma2=1,
    Pi = c(tab3$pi1[i], 1-tab3$pi1[i]),
    R=R, alternative="two.sided",
    power=0.90)$N))
}
cbind(tab3, m)