--- title: "Figure 3 source code eLife submission" output: html_document --- --- title: "Temperature paper analysis" output: html_document: default pdf_document: default --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ```{r message = FALSE, warning = FALSE} library(easyGgplot2) #ggplot2.multiplot(..., plotlist=NULL, cols=2) library(easyGgplot2) library(lme4) library(merTools) combinedData <- read.csv("~/Desktop/Temperature study/combined_data.csv") combinedData$time_HR_imputed <- combinedData$time_HR combinedData$time_HR_imputed[is.na(combinedData$time_HR_imputed)] <- 12 dataVeterans <- read.csv('~/Desktop/Temperature study/veterans_processed.csv') dataNHANES <- read.csv('~/Desktop/Temperature study/NHANES_processed.csv') ST_data <- read.csv('~/Desktop/Temperature study/STRIDE_processed.csv') dataVeterans$race <- factor(as.character(dataVeterans$race),levels = c("white", "black")) combinedData$race [combinedData$race=="asian"]<- "other" combinedData$race <- factor(combinedData$race,levels = c("white", "black","other" )) dataNHANES$race[dataNHANES$race == "unknown"] <- NA dataNHANES$race <- factor(dataNHANES$race,levels = c("white", "black","other" )) ST_data$race [ST_data$race=="Asian"]<- "Other" ST_data$race <- factor(ST_data$race,levels = c("White", "Black","Other")) ST_data$exam_month<-factor (ST_data$exam_month) dataVeterans$state_type <- factor(dataVeterans$state_type, levels = c("very cold", "cold","moderately_cold","moderate","warm", "hot")) dataNHANES$exammonth <-as.factor(dataNHANES$exammonth) dataVeterans$exammonth <-as.factor(dataVeterans$exammonth) state_temps <- read.csv("~/Desktop/Temperature study/state_temps.csv") state_codes <- read.csv("~/Desktop/Temperature study/state_codes.csv") state_codes$Code <- factor(as.character(state_codes$Code), levels=levels(dataVeterans$state)) state_temps$state_code <- state_codes$Code[match(as.character(state_temps$state),as.character(state_codes$State))] dataVeterans$ambient_temp <- (state_temps$value[match(paste(dataVeterans$state,(dataVeterans$examyear * 100 + as.numeric(dataVeterans$exammonth))), paste(state_temps$state_code,state_temps$date))] - 32) / 1.8 ST_data$state <- "CA" ST_data$ambient_temp <- (state_temps$value[match(paste(ST_data$state,(ST_data$exam_year * 100 + as.numeric(ST_data$exam_month))), paste(state_temps$state_code,state_temps$date))] - 32) / 1.8 ``` Figure 3: Modeled body temperature over time in three cohorts by birth year (black and white ethnicity groups). Body temperature decreases by birth year in white and black men and women. No data for women were available for the birth years from 1800-1890 ```{r message = FALSE, warning = FALSE} combined_subset <- combinedData$sex=="male" & combinedData$race =="white" & !is.na(combinedData$race) & !is.na(combinedData$height_CM) & !is.na(combinedData$weight_KG) & !is.na(combinedData$sex) & combinedData$sample_weights > 0 combined_subset_f <- combinedData$sex=="female" & combinedData$race =="white" & !is.na(combinedData$race) & !is.na(combinedData$height_CM) & !is.na(combinedData$weight_KG) & !is.na(combinedData$sex) & combinedData$sample_weights > 0 linear <-lm(temp_C ~ age + weight_KG + height_CM + I(birth_year-1800), data=combinedData[combined_subset, ] ) summary(linear) ``` ```{r message = FALSE, warning = FALSE} linear_f <-lm(temp_C ~ age + weight_KG + height_CM + I(birth_year-1800), data=combinedData[combined_subset_f, ] ) summary(linear_f) ``` ```{r message = FALSE, warning = FALSE} pred_df <- data.frame(age=rep(30,201), birth_year=seq(1800,2000), weight_KG = rep(70,201), height_CM = rep(170,201), race=rep("white",201), sex=rep("male",201)) pred_df$mean <- (predict(linear, pred_df,interval="confidence"))[,1] pred_df$ci025 <- (predict(linear, pred_df,interval="confidence"))[,2] pred_df$ci975 <- (predict(linear, pred_df,interval="confidence"))[,3] pred_df_f <- data.frame(age=rep(30,111), birth_year=seq(1890,2000), weight_KG = rep(70,111), height_CM = rep(170,111), race=rep("white",111), sex=rep("female",111)) pred_df_f$mean <- (predict(linear_f, pred_df_f,interval="confidence"))[,1] pred_df_f$ci025 <- (predict(linear_f, pred_df_f,interval="confidence"))[,2] pred_df_f$ci975 <- (predict(linear_f, pred_df_f,interval="confidence"))[,3] predictions <- rbind(pred_df, pred_df_f) ggplot() + geom_ribbon(data = pred_df, aes(birth_year, ymin=ci025,ymax=ci975),fill="blue",alpha=0.3) + geom_line(data = pred_df, aes(birth_year, mean),colour="blue") + geom_ribbon(data = pred_df_f, aes(birth_year, ymin=ci025,ymax=ci975),alpha=0.3,fill="red") + geom_line(data = pred_df_f, aes(birth_year, mean),colour="red") + labs(x= "Year of birth", y = "Temperature,°C") + scale_fill_identity(name = 'the fill', guide = 'legend',labels = c('m1')) + scale_colour_manual(name = 'the colour', values =c("blue"="blue","red"="red"), labels = c("Men","Women")) +ylim(36,37.5) ``` ```{r message = FALSE, warning = FALSE} predictions$mean[predictions$birth_year==1800 & predictions$sex == "male"] predictions$mean[predictions$birth_year==2000 & predictions$sex == "male"] (predictions$mean[predictions$birth_year==1800 & predictions$sex == "male"]- predictions$mean[predictions$birth_year==2000 & predictions$sex == "male"]) ``` ```{r message = FALSE, warning = FALSE} predictions$mean[predictions$birth_year==1890 & predictions$sex == "female"] predictions$mean[predictions$birth_year==2000 & predictions$sex == "female"] (predictions$mean[predictions$birth_year==1890 & predictions$sex == "female"]- predictions$mean[predictions$birth_year==2000 & predictions$sex == "female"]) ``` ```{r message = FALSE, warning = FALSE} combined_subset_b <- combinedData$sex=="male" & combinedData$race =="black" & !is.na(combinedData$race) & !is.na(combinedData$height_CM) & !is.na(combinedData$weight_KG) & !is.na(combinedData$sex) & combinedData$sample_weights > 0 linear_b <-lm(temp_C ~ age + weight_KG + height_CM + I(birth_year-1800), data=combinedData[combined_subset_b, ] ) summary(linear_b) ``` ```{r message = FALSE, warning = FALSE} combined_subset_b_f <- combinedData$sex=="female" & combinedData$race =="black" & !is.na(combinedData$race) & !is.na(combinedData$height_CM) & !is.na(combinedData$weight_KG) & !is.na(combinedData$sex) & combinedData$sample_weights > 0 linear_b_f <-lm(temp_C ~ age + weight_KG + height_CM + I(birth_year-1800), data=combinedData[combined_subset_b_f, ] ) summary(linear_b_f) ``` ```{r message = FALSE, warning = FALSE} pred_df_b <- data.frame(age=rep(30,201), birth_year=seq(1800,2000), weight_KG = rep(70,201), height_CM = rep(170,201), race=rep("black",201), sex=rep("male",201)) pred_df_b$mean <- (predict(linear_b, pred_df_b,interval="confidence"))[,1] pred_df_b$ci025 <- (predict(linear_b, pred_df_b,interval="confidence"))[,2] pred_df_b$ci975 <- (predict(linear_b, pred_df_b,interval="confidence"))[,3] pred_df_b_f <- data.frame(age=rep(30,111), birth_year=seq(1890,2000), weight_KG = rep(70,111), height_CM = rep(170,111), race=rep("black",111), sex=rep("female",111)) pred_df_b_f$mean <- (predict(linear_b_f, pred_df_b_f,interval="confidence"))[,1] pred_df_b_f$ci025 <- (predict(linear_b_f, pred_df_b_f,interval="confidence"))[,2] pred_df_b_f$ci975 <- (predict(linear_b_f, pred_df_b_f,interval="confidence"))[,3] predictions_b <- rbind(pred_df_b, pred_df_b_f) predictions_full <- rbind(predictions,predictions_b) ggplot(data = predictions_full) + geom_ribbon(aes(birth_year, ymin=ci025,ymax=ci975,fill=sex),alpha=0.3) + geom_line(aes(birth_year, mean,colour=sex),show.legend=FALSE) + scale_colour_manual(values = plot_colors[c(9,8)]) + scale_fill_manual(values = plot_colors [c(9,8 )], name="Sex",labels=c("Men","Women")) + facet_grid(. ~ race ) + labs(x= "Year of birth", y = "Temperature,°C") +ylim(36,37.5) + theme(legend.position="bottom",axis.title=element_text(size=14,face="bold"), legend.title=element_text(size=14), legend.text=element_text(size=14)) ``` Table B. Coefficients (and standard errors) used for the graph from multivariate linear regression including age, body weight, height and birth year. All cells are significant at greater than 99% significance level. ```{r message = FALSE, warning = FALSE} predictions_b$mean[predictions_b$birth_year==1800 & predictions_b$sex == "male"] predictions_b$mean[predictions_b$birth_year==2000 & predictions_b$sex == "male"] (predictions_b$mean[predictions_b$birth_year==1800 & predictions_b$sex == "male"]- predictions_b$mean[predictions_b$birth_year==2000 & predictions_b$sex == "male"]) ``` ```{r message = FALSE, warning = FALSE} predictions_b$mean[predictions_b$birth_year==1890 & predictions_b$sex == "female"] predictions_b$mean[predictions_b$birth_year==2000 & predictions_b$sex == "female"] (predictions_b$mean[predictions_b$birth_year==1890 & predictions_b$sex == "female"]- predictions_b$mean[predictions_b$birth_year==2000 & predictions_b$sex == "female"]) ```