Repeated measures (mixed model) GWAS about bigstatsr HOT 5 CLOSED

library(ggplot2)
library(dplyr)
library(tidyr)
library(faux)  
library(GGally)

# simulating PCs from genotype data
pc <- 
  rnorm_multi(
    n    = 10000,
    mu   = 0,
    sd   = 1,
    r    = 0,
    varnames = c(paste0("PC",1:10)),
    empirical = F
  ) %>%
  mutate(IID  = row_number())

# simulating longitudinal outcome data with age covariate
df <- 
  rnorm_multi(
    n    = 10000,
    mu   = 50,
    sd   = 10,
    r    = 0.4,
    varnames = c(paste0("Y_Q0",1:7)),
    empirical = F
  ) %>%
  mutate(
    IID  = row_number(),
    Age_Q01 = rnorm(10000, mean=35, sd=5) 
  ) %>%
  mutate(
    Age_Q02 = Age_Q01 + 0.5,
    Age_Q03 = Age_Q01 + 1,
    Age_Q04 = Age_Q01 + 3,
    Age_Q05 = Age_Q01 + 5,
    Age_Q06 = Age_Q01 + 7,
    Age_Q07 = Age_Q01 + 8
  ) 

# adjusting phenotype for genotype PCs (this 2 step process is done in regenie)
yRes <-
  df %>% 
  full_join(pc, "IID") %>%
  mutate(
    across(matches("Y"), 
    ~ rstandard(lm(.x~PC1+PC2+PC3+PC4+PC5+PC6+PC7+PC8+PC9+PC10,df))
    )) %>%
  select(IID, matches("Y_|Age"))%>% 
  pivot_longer(
    cols = !IID, 
    names_to = c(".value", "Wave"), 
    names_sep = "_"
  ) 

# adjusting genotypes for PCs (again... this 2 step process is done in regenie)
gRes <- 
  rnorm_multi(
    n    = 10000,
    mu   = 0,
    sd   = 1,
    r    = 0,
    varnames = c(paste0("PC",1:10)),
    empirical = F
  ) %>%
  mutate(
    IID  = row_number(),
    SNP  = rbinom(10000, 2, 0.5)
  ) %>%
  mutate(
    SNPRes = rstandard(lm(as.formula(paste0("SNP ~ ",paste0("PC",1:10, collapse="+"))), data = .))
  ) %>%
  select(IID,SNPRes)

df <-  yRes %>% full_join(gRes)

# run mixed model... for each SNP I want the main effect of SNP and its interaction with age 
output <- 
  summary(
    lme4::lmer(Y ~ SNPRes + SNPRes*Age + (1|IID), data=df)
    )$coefficients