powerLME {powerEQTL}R Documentation

Power Calculation for Simple Linear Mixed Effects Model

Description

Power calculation for simple linear mixed effects model. This function can be used to calculate one of the 3 parameters (power, sample size, and minimum detectable slope) by setting the corresponding parameter as NULL and providing values for the other 2 parameters.

Usage

powerLME(
  slope, 
  n, 
  m, 
  sigma.y, 
  sigma.x,
  power = NULL,
  rho = 0.8, 
  FWER = 0.05,
  nTests = 1,
  n.lower = 2.01,
  n.upper = 1e+30)

Arguments

slope

numeric. Slope under alternative hypothesis.

n

integer. Total number of subjects.

m

integer. Number of observations per subject.

sigma.y

numeric. Standard deviation of the outcome y.

sigma.x

numeric. Standard deviation of the predictor x.

power

numeric. Desired power.

rho

numeric. Intra-class correlation (i.e., correlation between yijy_{ij} and yiky_{ik} for the jj-th and kk-th observations of the ii-th subject).

FWER

numeric. Family-wise Type I error rate.

nTests

integer. Number of tests (e.g., number of genes in differential expression analysis based on scRNAseq to compare gene expression between diseased subjects and healthy subjects).

n.lower

numeric. Lower bound of the total number of subjects. Only used when calculating toal number of subjects.

n.upper

numeric. Upper bound of the total number of subjects. Only used when calculating total number of subjects.

Details

We assume the following simple linear mixed effects model to characterize the association between the predictor x and the outcome y:

yij=β0i+β1xi+ϵij,y_{ij} = \beta_{0i} + \beta_1 * x_i + \epsilon_{ij},

where

β0iN(β0,σβ2),\beta_{0i} \sim N\left(\beta_0, \sigma^2_{\beta}\right),

and

ϵijN(0,σ2),\epsilon_{ij} \sim N\left(0, \sigma^2\right),

i=1,,ni=1,\ldots, n, j=1,,mj=1,\ldots, m, nn is the number of subjects, mm is the number of observations per subject, yijy_{ij} is the outcome value for the jj-th observation of the ii-th subject, xix_i is the predictor value for the ii-th subject. For example, xix_i is the binary variable indicating if the ii-th subject is a diseased subject or not.

We would like to test the following hypotheses:

H0:β1=0,H_0: \beta_1=0,

and

H1:β1=δ,H_1: \beta_1 = \delta,

where δ0\delta\neq 0.

We can derive the power calculation formula is

power=1Φ(zα/2a×b)+Φ(zα/2a×b),power=1- \Phi\left(z_{\alpha^{*}/2}-a\times b\right) +\Phi\left(-z_{\alpha^{*}/2} - a\times b\right),

where

a=σ^xσya= \frac{\hat{\sigma}_x }{\sigma_y}

and

b=δm(n1)1+(m1)ρ b=\frac{\delta\sqrt{m(n-1)}}{\sqrt{1+(m-1)\rho}}

and zα/2z_{\alpha^{*}/2} is the upper 100α/2100\alpha^{*}/2 percentile of the standard normal distribution, α=α/nTests\alpha^{*}=\alpha/nTests, nTests is the number of tests, σy=σβ2+σ2\sigma_y=\sqrt{\sigma^2_{\beta}+\sigma^2}, σ^x=i=1n(xixˉ)2/(n1)\hat{\sigma}_x=\sqrt{\sum_{i=1}^n\left(x_i-\bar{x}\right)^2/(n-1)}, and ρ=σβ2/(σβ2+σ2)\rho=\sigma^2_{\beta}/\left(\sigma^2_{\beta}+\sigma^2\right) is the intra-class correlation.

Value

power if the input parameter power = NULL.

sample size (total number of subjects) if the input parameter n = NULL;

minimum detectable slope if the input parameter slope = NULL.

Author(s)

Xianjun Dong <XDONG@rics.bwh.harvard.edu>, Xiaoqi Li<xli85@bwh.harvard.edu>, Tzuu-Wang Chang <Chang.Tzuu-Wang@mgh.harvard.edu>, Scott T. Weiss <restw@channing.harvard.edu>, Weiliang Qiu <weiliang.qiu@gmail.com>

References

Dong X, Li X, Chang T-W, Scherzer CR, Weiss ST, and Qiu W. powerEQTL: An R package and shiny application for sample size and power calculation of bulk tissue and single-cell eQTL analysis. Bioinformatics, 2021;, btab385

Examples

  n = 102
  m = 227868
  
  # calculate power
  power = powerLME(
    slope = 0.6, 
    n = n, 
    m = m,
    sigma.y = 0.29, 
    sigma.x = 0.308,
    power = NULL,
    rho = 0.8, 
    FWER = 0.05,
    nTests = 1e+6)

  print(power)
  
  # calculate sample size (total number of subjects)
  n = powerLME(
    slope = 0.6, 
    n = NULL, 
    m = m,
    sigma.y = 0.29, 
    sigma.x = 0.308,
    power = 0.9562555,
    rho = 0.8, 
    FWER = 0.05,
    nTests = 1e+6)

  print(n)
  
  # calculate slope
  slope = powerLME(
    slope = NULL, 
    n = n, 
    m = m,
    sigma.y = 0.29, 
    sigma.x = 0.308,
    power = 0.9562555,
    rho = 0.8, 
    FWER = 0.05,
    nTests = 1e+6)

  print(slope)
  



[Package powerEQTL version 0.3.4 Index]