---
title: "Figure 1 eLife submission"
output: pdf_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 1 A: Body temperature measurements by age as observed in three different time periods (cohorts)

```{r message = FALSE, warning = FALSE}
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
#Names("Black", "Orange", "Cyan", "Lt.Green", "Dk.Gray", "Red", "Blue", "Green", "Brown", "Purple", "Lt.Gray", "Yellow", "Pink”))
#c("#000000", "#FF9233", "#29D0D0", "#81C57A", "#575757", "#AD2323", "#2A4BD7", "#1D6914", "#814A19", "#8126C0", "#A0A0A0", "#FFEE33", "#FFCDF3")

plot_colors <- c("#FF9233", "#29D0D0", "#81C57A", "#575757", "#8126C0", "#29D0D0" ,"#575757", "#AD2323", "#2A4BD7", "#1D6914", "#814A19", "#8126C0", "#A0A0A0", "#FFEE33", "#FFCDF3") # gg_color_hue(3)

combinedData$plot_sex[combinedData$sex=="male"] <- "men"
 combinedData$plot_sex[combinedData$sex=="female"] <- "women"
 combinedData$plot_sex <- factor(combinedData$plot_sex,levels=c("men","women"))
 
 Combined <- ggplot(data = combinedData[ combinedData$age  >30  &  combinedData$age
 <80 & combinedData$race %in% c("white","black"),],aes(age, temp_C, group=period,colour=period)) + geom_smooth() + scale_colour_manual(values = plot_colors[c(2,3,1 )
 ]) + ylim(36,37.5) + facet_grid(plot_sex ~ race)
 
 Combined + theme(legend.position="bottom", axis.text=element_text(size=14),
         axis.title=element_text(size=14,face="bold"), legend.text=element_text(size=12), legend.title=element_text(size=14)) + labs(x= "Age, years", y = "Temperature,°C in white men(smoothed unadjusted data)", colour = "Time period of the measuremement")

```

Figure 1 Table B and C: Body temperature measurements by age as observed in three different time periods (cohorts)

B. Coefficients and standard errors from multivariate linear regression models for each cohort including age, weight, height, ethnicity group and time of day as available. Yellow cells are statistically significant at a p value of <0.01, orange cells are of borderline significance (p<0.1 but >0.05), and remaining uncolored cells are not statistically significant. C. Expected body temperature for 30-year old men and women with weight 70 kg and height 170 cm in each time period/cohort. 

We run separate regressions for each cohort and present it in a table in results : lmVets- regression from Union Army veterans cohort from 1860-1940, NHANES data from 1971-1975 and STRIDE data from 2000-2017. STRIDE and a part of NHANES also includes time of the day (Time_HR). The coeficients of regresstion models are then used in the  Figure 1 Table B. We also do predictions for men and women at age 30, body weight 70KG and height 170 CM, to compare the values across different cohorts. The coeficients are used in Figure 1 table C.


```{r message = FALSE, warning = FALSE}
lmVeterans <-lm(temp_C~ age + weightKG + heightCM + race, data= dataVeterans)
summary(lmVeterans)
lmVeterans_predict=data.frame(age=30, weightKG=70, heightCM = 170, race="white")

round(predict (lmVeterans, lmVeterans_predict, interval="confidence"),2)
```


NHANES analysis controlling for race and time of the day

In NHANES we included weights from the sample to try and get adjusted values, because we know that NHANES for example oversampled women in fertile age. But the coefficients in the weighted model were not very different from the initial model.

```{r message = FALSE, warning = FALSE}
lmNHANES_W_races <- lm(temp_C ~ age + weight_KG + height_CM + race+ time_HR, data=
 dataNHANES[dataNHANES$sex == "male" , ], weights=dataNHANES$sample_weights[dataNHANES$sex=="male"])
 summary(lmNHANES_W_races)
 
NHANES_predict=data.frame (age=30, weight_KG=70, height_CM = 170, race="white", sex ="male", time_HR=12)
round(predict(lmNHANES_W_races, NHANES_predict, interval="confidence"),2)

lmNHANES_W_races_F <- lm(temp_C ~ age + weight_KG + height_CM + race+ time_HR, data=dataNHANES[dataNHANES$sex == "female" , ], weights=dataNHANES$sample_weights[dataNHANES$sex=="female"])
summary(lmNHANES_W_races_F)

NHANES_predict=data.frame (age=30, weight_KG=70, height_CM = 170, race="white", sex ="female", time_HR=12)
 round(predict(lmNHANES_W_races_F, NHANES_predict, interval="confidence"),2)
```
 
