Toad diets, paralell regressions

running many single prediction regressions of oxidative stress

These Knitr Options are for appearance only

We always start with knitr options for building a good html report

Setup Code

Please modify the path below, to point at your working directory

#====================================================================#
# Setup Options
#====================================================================#

# remove all objects if restarting script
rm(list=ls())

# set tibble width for printing
options(tibble.width = Inf)

# remove scientific notation
options(scipen=999)

#==============================================================================#
# Install Packages / Load Packages
#==============================================================================#

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(lubridate)
library(readxl) #install if needed
library(DT)

Warning: package 'DT' was built under R version 4.5.3

#==============================================================================#
# Set paths
#==============================================================================#

# set all paths
path_main    <- "C:/Users/jsteibel/OneDrive/Documents/Cheryl"
path_data    <- str_c(path_main, "/", sep="")

# str_c() is from the stringr package, which is a part of 'tidyverse'
# this concatenates two pieces of the path. We will use it a lot.

#check:
path_data

[1] "C:/Users/jsteibel/OneDrive/Documents/Cheryl/"

Read data

First input the dataset and check

#==============================================================================#
# Set Inputs
#==============================================================================#

dt<-read_xlsx("Toad_OS_and_diet_comp_SAS.xlsx",na=".")

reformat to take means by treatment and species combination

Let’s read the results and describe them below. Also, rename column NA to “sodio” because NA is a forbidden name in R (used for missing data)

mn<-dt%>%select(Spp,Trt,GPx,PtC,DM:Tcar)%>%pivot_longer(
  cols=GPx:Tcar,names_to = 'variable',values_to = 'values'
)%>%group_by(Spp,Trt,variable)%>%summarize(mn=mean(values,na.rm=TRUE))

`summarise()` has grouped output by 'Spp', 'Trt'. You can override using the
`.groups` argument.

mn<-pivot_wider(mn,names_from = variable,values_from = mn)

mn<-rename(mn,sodio='NA')

Data seems OK at least, structured as expected

Regression to predict OS markers

Setup

Indicate here what are the OS markers and select predictors and responses

target_var<-c("GPx","PtC")

vrs<-colnames(mn)[-(1:2)]
predictor_vars<-vrs[!(vrs%in%target_var)]
target_var

[1] "GPx" "PtC"

predictor_vars

 [1] "CA"      "CAP"     "CP"      "CU"      "DM"      "FAT"     "FE"     
 [8] "K"       "MEmod"   "MEunmod" "MG"      "MN"      "sodio"   "OM"     
[15] "P"       "PF"      "S"       "SE"      "TDF"     "Tcar"    "VitE"   
[22] "ZN"      "bcar"    "lut"     "retinol" "zeax"   

Predict first OS marker

model_results <- tibble(predictor = predictor_vars) %>%
  mutate(
    # Create formulas dynamically
    formula_str = map_chr(predictor, ~ paste(target_var[1], "~", .x)),
    
    # Fit the linear models
    model = map(formula_str, ~ lm(as.formula(.x), data = mn)),
    
    # Extract model-level statistics (R-squared, p-value, etc.)
    glance = map(model, broom::glance),
    
    # Extract coefficient-level statistics (Estimates, std error, etc.)
    tidy = map(model, broom::tidy)
  )

slope_results <- model_results %>%
  unnest(tidy) %>%
  filter(term != "(Intercept)") %>%
  select(predictor, estimate, std.error, p.value)

# Combine metrics and include the slope's p-value
final_summary <- model_results %>%
  # 1. Unnest glance and drop conflicting columns immediately
  unnest(glance) %>% 
  select(predictor, r.squared, sigma, tidy) %>% 
  
  # 2. Unnest tidy safely
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  
  # 3. Select and rename columns, keeping the slope's p-value
  select(predictor, r.squared, sigma, slope = estimate, p.value)

# View the final combined table

datatable(final_summary)%>%formatRound(columns=c('r.squared','sigma','slope','p.value'),digits = 3)

Predict 2nd OS marker

model_results <- tibble(predictor = predictor_vars) %>%
  mutate(
    # Create formulas dynamically
    formula_str = map_chr(predictor, ~ paste(target_var[2], "~", .x)),
    
    # Fit the linear models
    model = map(formula_str, ~ lm(as.formula(.x), data = mn)),
    
    # Extract model-level statistics (R-squared, p-value, etc.)
    glance = map(model, broom::glance),
    
    # Extract coefficient-level statistics (Estimates, std error, etc.)
    tidy = map(model, broom::tidy)
  )

slope_results <- model_results %>%
  unnest(tidy) %>%
  filter(term != "(Intercept)") %>%
  select(predictor, estimate, std.error, p.value)

# Combine metrics and include the slope's p-value
final_summary <- model_results %>%
  # 1. Unnest glance and drop conflicting columns immediately
  unnest(glance) %>% 
  select(predictor, r.squared, sigma, tidy) %>% 
  
  # 2. Unnest tidy safely
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  
  # 3. Select and rename columns, keeping the slope's p-value
  select(predictor, r.squared, sigma, slope = estimate, p.value)

# View the final combined table

datatable(final_summary)%>%formatRound(columns=c('r.squared','sigma','slope','p.value'),digits = 3)