Extract categorical coeffients and all p-values from a mixed model into a data table

user_18880691 0
Viewed 36

Here is a reproduceable code and sample data I want to achieve a final data table with 3 columns: 1. exposure quantile 2. OR/RR 3. PV

set.seed(42) 
n <- 100
dat = data.frame(ID = rep(c(1:25),times=4 ) ,
                 Score = rnorm(n, mean=0.3, sd=0.8))
                 
dat = dat %>%
  group_by(ID)%>%
  dplyr::mutate(exposure1 = rep(c(rnorm(1, mean=6, sd=1.8))),
                exposure2 = rep(c(rnorm(1, mean=3, sd=0.6))),
                age = rep(c(rnorm(1, mean=40, sd=15))))%>%
  ungroup()%>%
  dplyr::mutate(exposure1_quantile = cut(exposure1, breaks = 4, labels = c("Q1","Q2","Q3","Q4")),
                exposure2_quantile = cut(exposure2, breaks = 4, labels = c("Q1","Q2","Q3","Q4")))

exposures_var = c("exposure1_quantile","exposure2_quantile")
exposure_var_labels("exposure1 Q1","exposure1 Q2 ", "exposure 1 Q3", 
                    "exposure2 Q1","exposure2 Q2 ", "exposure2 Q3")
age="age"
outcome = "Score"
exposure_data_table = c()

for(i in 1:length(exposures_var)){
  exp = exposures_var[i]
  fixed_effects_formula = paste0(outcome, "~",exp,"+",age)
  fixed_effects_formula = as.formula(fixed_effects_formula)
  mixedmodel = lme(fixed =fixed_effects_formula, random = ~1|ID, data=dat, method = "ML")
  for(m in 2:4){
    v = mixedmodel$coefficients$fixed[m]
    
    vector = c(exp , v)
    #P=p value for every quantile (HOW TO ADD?)
    #exposure_name = exposure_var_labels[?] (HOW TO ADD LABEL)
    exposure_data_table = rbind(exposure_data_table, vector)
    
   
  }
  
}

exposure_data_table = as.data.table(exposure_data_table)
colnames(exposure_data_table)=c("Exposure","RR")#,"pv")
view(exposure_data_table)

I first used anova to try and get the pvalue but it didnt work.

mixed-models r for-loop
1 Answers

I think a tidymodels approach using lme would work well here:

library(nlme)
library(tidymodels)
library(multilevelmod)
library(data.table)

lme_spec <- 
  linear_reg() %>% 
  set_engine("lme", random = ~ 1 | ID)

Map(function(exp) {

  fixed_effects_formula <- as.formula(paste0("Score~",exp,"+ age +", 0))
 
  lme_spec %>% 
    fit(fixed_effects_formula, data = dat) %>%
    broom.mixed::tidy() %>%
    filter(effect == "fixed", grepl("exposure", term)) %>%
    select(term, estimate, std.error, p.value)
  }, exposures_var) %>%
  bind_rows() %>%
  as.data.table()
#>                    term    estimate std.error   p.value
#> 1: exposure1_quantileQ1 -0.16147364 0.3532834 0.6525497
#> 2: exposure1_quantileQ2  0.22318505 0.2719366 0.4214784
#> 3: exposure1_quantileQ3  0.24976757 0.3484126 0.4817411
#> 4: exposure1_quantileQ4  0.14177064 0.4020702 0.7280757
#> 5: exposure2_quantileQ1  0.28976458 0.4191198 0.4972840
#> 6: exposure2_quantileQ2  0.19907863 0.2699164 0.4693496
#> 7: exposure2_quantileQ3  0.35040767 0.2827229 0.2295436
#> 8: exposure2_quantileQ4 -0.09587234 0.3533819 0.7889412

Created on 2022-08-07 by the reprex package (v2.0.1)

user_125003150
Related