---
title: "Figure supplement 4 data source eLife submission"
output: html_document
---

```{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 1-figure supplement 4. Model of mean body temperature over time in three cohorts by birth year controlling for the time of the day of the temperature measurement (white men and women).
A. Graph shows a decrease of predicted body temperature by birth year for both white men (blue line) and white women (red line) based on the linear regression coefficients displayed in B. No data for women were available for the period 1800-1890. For all measurements from UACWV, where time of the day was not available, the values of 12:00 PM (noon) were imputed.
B. Coefficients (and standard errors) from linear regression model including age, weight in kg, height in cm and birth year.

```{r message = FALSE, warning = FALSE}
linear_HR <-lm(temp_C ~ age + weight_KG + height_CM + time_HR_imputed+ I(birth_year-1800),  data=combinedData[combined_subset, ] )
summary(linear_HR)
```
```{r message = FALSE, warning = FALSE}
linear_f_HR <-lm(temp_C ~ age + weight_KG + height_CM +time_HR_imputed + I(birth_year-1800),  data=combinedData[combined_subset_f, ] )
summary(linear_f_HR)
```
```{r message = FALSE, warning = FALSE}
pred_df_HR <- 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), time_HR_imputed=12)
pred_df_HR$mean <- (predict(linear_HR, pred_df_HR,interval="confidence"))[,1]
pred_df_HR$ci025 <- (predict(linear_HR, pred_df_HR,interval="confidence"))[,2]
pred_df_HR$ci975 <- (predict(linear_HR, pred_df_HR,interval="confidence"))[,3]

pred_df_f_HR <- 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),time_HR_imputed=12)
pred_df_f_HR$mean <- (predict(linear_f_HR, pred_df_f_HR,interval="confidence"))[,1]
pred_df_f_HR$ci025 <- (predict(linear_f_HR, pred_df_f_HR,interval="confidence"))[,2]
pred_df_f_HR$ci975 <- (predict(linear_f_HR, pred_df_f_HR,interval="confidence"))[,3]


predictions_HR <- rbind(pred_df_HR, pred_df_f_HR)


ggplot() + geom_ribbon(data = pred_df_HR, aes(birth_year, ymin=ci025,ymax=ci975),fill="blue",alpha=0.3) + geom_line(data = pred_df_HR, aes(birth_year, mean),colour="blue") + geom_ribbon(data = pred_df_f_HR, aes(birth_year, ymin=ci025,ymax=ci975),alpha=0.3,fill="red") + geom_line(data = pred_df_f_HR, 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_HR$mean[predictions_HR$birth_year==1800 & predictions_HR$sex == "male"]
 predictions_HR$mean[predictions_HR$birth_year==2000 & predictions_HR$sex == "male"]
(predictions_HR$mean[predictions_HR$birth_year==1800 & predictions_HR$sex == "male"]- predictions_HR$mean[predictions_HR$birth_year==2000 & predictions_HR$sex == "male"])
```
 
```{r message = FALSE, warning = FALSE}
 predictions_HR$mean[predictions_HR$birth_year==1890 & predictions_HR$sex == "female"]
 predictions_HR$mean[predictions_HR$birth_year==2000 & predictions_HR$sex == "female"]
(predictions_HR$mean[predictions_HR$birth_year==1890 & predictions_HR$sex == "female"]- predictions_HR$mean[predictions_HR$birth_year==2000 & predictions_HR$sex == "female"])
```
