TLDR: Given a dataset as shown in the last plot, we consider the following analyses:
1. For each subject extract the slope and perform a simple t-test on the slopes
2. Fit lmer(salery ~ slope + (1|subject)) and test if \(\beta_{slope} =0\)
3. Fit lmer(salery ~ slope + (slope|subject)) and test if \(\beta_{slope} =0\)
Then: 1 & 3 are valid (if data balanced, they are equally powerful) approaches and 2 is only valid if there is no individual slope.

Simulate Data

Salary of subject \(s\) at time \(t\):

\[ salary_{t}^{(s)} = t \cdot (slope + d_s) + intercept_s + \varepsilon_{t,s} \]

\(d_s \sim \mathcal N(0, subjSlopeSD^2)\)
\(intercept_s \sim \mathcal N(0, subjSD^2)\)
\(\varepsilon_{t,s} \sim \mathcal N(0, obsSD^2)\)

nsim <- 2000

library(mcreplicate) # for parallelization
library(ggplot2)
library(lmerTest)

Loading required package: lme4

Loading required package: Matrix


Attaching package: 'lmerTest'

The following object is masked from 'package:lme4':

    lmer

The following object is masked from 'package:stats':

    step

#' `nsub`: how many subjects (default: 6)
#' `nyears`: how many years(default: 10)
#' `obsSD`: standard deviation of noise (observation-level) (default: 15)
#' `subjSD`: standard deviation of individual effect (default: 4)
#' `slope`: shared increase of income per year (default: 5)
#' `subjSlopeSD`: subject specific standard deviation from slope (default: 2)
get_data <- function(nsub=6, nyears=10, 
                     obsSD=15, subjSD=4, 
                     slope=5, subjSlopeSD=2){
  subj_intercept <- rep(rnorm(nsub, 0, subjSD), each=nyears)
  subj_slope <- rep(rnorm(nsub, slope, subjSlopeSD), each=nyears)
  data.frame(
    subject = as.factor(rep(1:nsub, each = nyears)),
    year = rep(1:nyears, times = nsub),
    salary =  subj_intercept +  # subject effect
              subj_slope*(1:nyears) + # individual slope
              rnorm(nyears*nsub, 0, obsSD) # obsSD
  )
}

Plot Data

plot_data <- function(main, ...){
  ggplot(get_data(...), aes(x = year, y = salary, group = subject, color = subject)) +
    geom_line() +
    geom_point() +
    labs(x = "Year", y = "Salary") +
    theme_minimal() +
    ggtitle(main)
}

set.seed(123)
plot_data("Noisless, equal subjects", obsSD=0, subjSD=0, subjSlopeSD=0)

plot_data("Noisless, equal slopes", obsSD=0, subjSlopeSD=0)

plot_data("Noisless", obsSD=0)

plot_data("General case")

Analysis Methods

slopeTTest <- function(data){
  fits <- lmList(salary ~ year | subject, data)
  slopes <- coef(fits)[,"year"]
  t_test <- t.test(slopes, mu = 0)
  t_test$p.value
}

lmmRandItcpt <- function(data){
  lmm <- lmer(salary ~ year + (1|subject), data = data)
    summary(lmm)$coefficients[2, "Pr(>|t|)"]
}

lmmRandSlope <- function(data){
  lmm <- lmer(salary ~ year + (year|subject), data = data)
    summary(lmm)$coefficients[2, "Pr(>|t|)"]
}

p_values <- function(...){
  data <- get_data(...)
  c(
    slopeTTest = slopeTTest(data),
    lmmRandItcpt = lmmRandItcpt(data),
    lmmRandSlope = lmmRandSlope(data)
  ) # if you cange the amount of arguments, change also the object "Power"
}

get_power <- function(...){
  args <- list(...)
  PVALS <- as.data.frame(t(
  mc_replicate(nsim, do.call(p_values, args))
  ))
  # POWER:
  sapply(lapply(PVALS, function(x) x<0.05), mean)
}

Power Calculations

set.seed(123)
ARGS <- as.data.frame(rbind(
  # NULL:
  expand.grid( 
    nsub=6, 
    nyears=10, 
    obsSD=5, 
    subjSD=4, 
    slope=0, 
    subjSlopeSD=c(0,4)),
  # Change standard deviations and effects:
  expand.grid(
    nsub=6, 
    nyears=10, 
    obsSD=c(5,15), 
    subjSD=c(4, 12), 
    slope=c(2,4), 
    subjSlopeSD=c(1,4)),
  # Change samplesize (allocation):
  expand.grid(
    nsub=c(6,10,20), 
    nyears=c(10,20), 
    obsSD=5, 
    subjSD=4, 
    slope=2, 
    subjSlopeSD=1)[-1,]
))
rownames(ARGS) <- NULL

Power <- matrix(NA, nrow=nrow(ARGS), ncol=3)
colnames(Power) <- c("slopeTTest", 
  "lmmRandItcpt", "lmmRandSlope")

for (i in 1:nrow(ARGS)){
  args <- as.list(ARGS[i,])
  Power[i,] <- do.call(get_power, args)
}

Show Results:

as.data.frame(cbind(ARGS, Power)) |> knitr::kable()

nsub	nyears	obsSD	subjSD	slope	subjSlopeSD	slopeTTest	lmmRandItcpt	lmmRandSlope
6	10	5	4	0	0	0.0440	0.0435	0.0170
6	10	5	4	0	4	0.0520	0.5415	0.0515
6	10	5	4	2	1	0.9200	0.9985	0.9185
6	10	15	4	2	1	0.5385	0.7970	0.5120
6	10	5	12	2	1	0.9175	0.9995	0.9170
6	10	15	12	2	1	0.5500	0.7880	0.5335
6	10	5	4	4	1	1.0000	1.0000	1.0000
6	10	15	4	4	1	0.9775	0.9995	0.9780
6	10	5	12	4	1	1.0000	1.0000	1.0000
6	10	15	12	4	1	0.9805	0.9995	0.9805
6	10	5	4	2	4	0.1820	0.7715	0.1820
6	10	15	4	2	4	0.1515	0.5910	0.1475
6	10	5	12	2	4	0.1570	0.7380	0.1570
6	10	15	12	2	4	0.1510	0.5815	0.1480
6	10	5	4	4	4	0.5000	0.9630	0.5015
6	10	15	4	4	4	0.4555	0.9010	0.4470
6	10	5	12	4	4	0.4905	0.9680	0.4905
6	10	15	12	4	4	0.4405	0.9165	0.4350
10	10	5	4	2	1	0.9990	1.0000	0.9990
20	10	5	4	2	1	1.0000	1.0000	1.0000
6	20	5	4	2	1	0.9655	1.0000	0.9655
10	20	5	4	2	1	1.0000	1.0000	1.0000
20	20	5	4	2	1	1.0000	1.0000	1.0000