---
title: "Figure supplement 2"
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 2: Analysis using BMI adjusted for height. 
Coefficients (and standard errors) are from multivariate linear regression models for each cohort including age, height, BMI adjusted for height and time of day (as available). NHANES I includes population weights. *, **, *** indicates significance at the 90%, 95%, and 99% level, respectively.


```{r message = FALSE, warning = FALSE}
lm_bmi_adj <- lm(BMI ~ heightCM, data= dataVeterans)
dataVeterans$bmi_adj <- residuals(lm_bmi_adj)
summary (lm_bmi_adj)

base_h <- 170
base_w <- 70
base_bmi <- base_w / (base_h/100)^2
base_bmi_adj <- base_bmi - (predict(lm_bmi_adj,data.frame(heightCM=base_h)))
 
 
```

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

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

```{r message = FALSE, warning = FALSE}
dataNHANES$bmi_adj <- residuals(lm(BMI ~ height_CM, data=dataNHANES))

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

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

```

```{r message = FALSE, warning = FALSE}
lmNHANES_BMI_adj <-lm(temp_C ~ age + height_CM  + bmi_adj + race +time_HR, data=dataNHANES[dataNHANES$sex == "male", ],weights=dataNHANES$sample_weights[dataNHANES$sex=="male" ])
summary(lmNHANES_BMI_adj)

NHANES_predict=data.frame (age=30, weight_KG=70, height_CM = 170, bmi_adj=base_bmi_adj, race="white", sex ="male", time_HR=12)
round(predict(lmNHANES_BMI_adj, NHANES_predict, interval="confidence"),2)

```


```{r message = FALSE, warning = FALSE}
ST_data$bmi_adj <- residuals(lm(BMI ~ height_CM, data=ST_data))
lmSTRIDE_BMI_adj_F <-lm(temp_C ~ age_years + height_CM  + bmi_adj + race +time_HR, data=ST_data[ST_data$GENDER  == "Female" , ])

summary (lmSTRIDE_BMI_adj_F)

STRIDE_predict=data.frame(age_years=30, weight_KG=70, height_CM = 170, time_HR=12, bmi_adj=base_bmi_adj, race="White", GENDER = "Female")

round(predict(lmSTRIDE_BMI_adj_F, STRIDE_predict, interval="confidence"),2)

```

```{r message = FALSE, warning = FALSE}
lmSTRIDE_BMI_adj <-lm(temp_C ~ age_years + height_CM  + bmi_adj + race + time_HR,
data=ST_data[ST_data$GENDER  == "Male" , ])

summary (lmSTRIDE_BMI_adj)


STRIDE_predict=data.frame(age_years=30, weight_KG=70, height_CM = 170, time_HR=12,
bmi_adj=base_bmi_adj, race="White", GENDER = "Male")
round(predict(lmSTRIDE_BMI_adj, STRIDE_predict, interval="confidence"),2)
```
