---
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"])
```
