Load packages

library (tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.4     ✔ tibble    3.1.8
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(ggplot2)
library(lme4)

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

library(glmmTMB)

Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
TMB was built with Matrix version 1.5.3
Current Matrix version is 1.5.1
Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package

Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
glmmTMB was built with TMB version 1.9.1
Current TMB version is 1.9.2
Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)

library(effects)

Loading required package: carData
lattice theme set by effectsTheme()
See ?effectsTheme for details.

library(car)

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

library(bbmle)

Loading required package: stats4

Attaching package: 'bbmle'

The following object is masked from 'package:dplyr':

    slice

library(performance)
library(tidyr)
library(readxl)
library(rcompanion)
library(viridis)

Loading required package: viridisLite

library(cowplot)

Attaching package: 'cowplot'

The following object is masked from 'package:lubridate':

    stamp

library(ggeffects)

Attaching package: 'ggeffects'

The following object is masked from 'package:cowplot':

    get_title

library(emmeans)
library(dplyr)
library(viridis)

Load data

ramping <-read.csv("/Users/madelinelohr/Desktop/thesis/thesis.bodytemp/body_temp_092123.csv")
behavior <- read.csv("/Users/madelinelohr/Desktop/thesis/thesis.bodytemp/body_temp_092123.csv")
data.rate<-read.csv("/Users/madelinelohr/Desktop/thesis/thesis.resp/respiration_rate2023.csv")

Thermoregulation Analysis

Data Cleaning

ramping1<- ramping[is.na(ramping$drop),]

ramping1$max_body_temp = as.numeric(ramping1$max_body_temp)
ramping1$infected = as.factor(ramping1$infected)
ramping1$set_temp_C = as.numeric(ramping1$set_temp_C)
ramping1$Ambient = as.character(ramping1$set_temp_C)
ramping1$replicate = as.numeric(ramping1$replicate )

thermoreg = (ramping1$max_body_temp - ramping1$set_temp_C)  
ramping1$size<-cut(ramping1$body_size, breaks=3)
ramping1$size = as.character(ramping1$size)

ramping2 <- cbind(ramping1,thermoreg)

ramping3=subset(ramping2, Ambient!="27")

ramping3$size<-cut(ramping3$body_size, breaks=3)
ramping3$size = as.character(ramping3$size)
#use "size" for figures, use "body_size" for stats



ramping4=subset(ramping3, Ambient!="23")

Linear Mixed Effects Model for Thermoregulation

m4.1 <- lmer(formula = thermoreg~set_temp_C+infected*body_size+(1|parent_colony)+(1|date), data = ramping4)

plot(allEffects(m4.1))

vif(m4.1)

        set_temp_C           infected          body_size infected:body_size 
          1.000094         104.880938           2.129224         111.983491 

check_model(m4.1)

Anova(m4.1)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: thermoreg
                      Chisq Df Pr(>Chisq)    
set_temp_C         986.9217  1  < 2.2e-16 ***
infected             3.8882  1  0.0486252 *  
body_size            9.8621  1  0.0016872 ** 
infected:body_size  13.3467  1  0.0002589 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Homeostatic Baseline

homebase = (ramping4$max_body_temp - 32)  
ramping5 <- cbind(ramping4,homebase)
m9 <- lmer(formula = homebase~set_temp_C+infected*body_size+(1|parent_colony)+(1|date), data = ramping5)
vif(m9)

        set_temp_C           infected          body_size infected:body_size 
          1.000094         104.880938           2.129224         111.983491 

Anova(m9)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: homebase
                      Chisq Df Pr(>Chisq)    
set_temp_C         835.6008  1  < 2.2e-16 ***
infected             3.8882  1  0.0486252 *  
body_size            9.8621  1  0.0016872 ** 
infected:body_size  13.3467  1  0.0002589 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Behavior Analysis

Data Cleaning

behavior1<- behavior[is.na(behavior$drop),]

behavior2 <- behavior1 %>% 
  filter(behavior != "N/A" | fanning == "N" | fanning == "Y")

behavior3 <- behavior2[-264,]

behavior4=subset(behavior3, set_temp_C!=23)

behavior5=subset(behavior4, set_temp_C!=27)
behavior5$max_body_temp = as.numeric(behavior5$max_body_temp)
behavior5$infected = as.factor(behavior5$infected)
behavior5$set_temp_C = as.numeric(behavior5$set_temp_C)
behavior5$Ambient = as.character(behavior5$set_temp_C)
behavior5$replicate = as.numeric(behavior5$replicate)

behavior6 <- behavior5
behavior6$behavior = as.factor(behavior6$behavior)
behavior6Ambient = as.factor(behavior6$Ambient)
moving = ifelse(behavior6$behavior=="M", "1", "0")

behavior7 <- cbind(behavior6, moving)
behavior7$moving = as.factor(behavior7$moving)
tempreg = (behavior7$max_body_temp - behavior7$set_temp_C)
behavior7 <- cbind(behavior7,tempreg)
behavior7$fanning = as.factor(behavior7$fanning)

behavior8 <- behavior7 
behavior8$fanning = ifelse(behavior8$fanning=="Y", "1", "0")
behavior8$fanning=as.numeric(behavior8$fanning)
behavior8 <- behavior8[,-20]
behavior8$moving = as.numeric(behavior8$moving)

Model Selection: Moving Behavior

moving_mod4<-lmer(moving ~ infected + set_temp_C + (1|parent_colony), data=behavior8)

move_lm <- lm(moving~infected+body_size+set_temp_C, data=behavior8)

move_lm2 <- lm(moving~infected+set_temp_C, data=behavior8)

ICtab(moving_mod4, move_lm, move_lm2)

            dAIC df
move_lm      0.0 5 
move_lm2     5.3 4 
moving_mod4 26.7 5 

vif(move_lm2)

  infected set_temp_C 
   1.00005    1.00005 

plot(allEffects(move_lm2))

check_model(move_lm2)

ANOVA on Selected Model for Moving

Anova(move_lm2)

Anova Table (Type II tests)

Response: moving
            Sum Sq  Df F value Pr(>F)
infected     0.001   1  0.0029 0.9573
set_temp_C   0.029   1  0.1471 0.7014
Residuals  139.081 713               

Model Selection: Fanning

behavior8$fanning_binary <- ifelse(behavior7$fanning== "N", 0, 1)

ramping4$fan_binary <- ifelse(ramping4$fanning== "N", 0, 1) 

sick_fan=subset(ramping4, infected==1)

fan_lm <- lm(fanning_binary ~ infected+set_temp_C+body_size, data=behavior8)

fan_lm2 <- lm(fanning_binary ~ infected+set_temp_C, data=behavior8)

fan_lm3 <- lm(fanning_binary ~ infected, data=behavior8)

fan_lm4 <- lm(fanning_binary ~ set_temp_C, data=behavior8)

fan_lm5 <- lm(fanning_binary~infected*set_temp_C, data=behavior8)

ICtab(fan_lm,fan_lm2,fan_lm3,fan_lm4,fan_lm5)

        dAIC df
fan_lm4  0.0 3 
fan_lm2  1.8 4 
fan_lm5  3.2 5 
fan_lm   3.6 5 
fan_lm3 45.1 3 

vif(fan_lm2)

  infected set_temp_C 
   1.00005    1.00005 

plot(allEffects(fan_lm2))

check_model(fan_lm2)

ANOVA on Selected Model for Fanning

Anova(fan_lm2)

Anova Table (Type II tests)

Response: fanning_binary
            Sum Sq  Df F value   Pr(>F)    
infected    0.0099   1  0.2473   0.6191    
set_temp_C  1.8757   1 46.6263 1.84e-11 ***
Residuals  28.6822 713                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Respiration Analysis

Data Cleaning

data.rate$rate = (data.rate$rate)*(-1)
data.rate$rate2 = data.rate$rate/data.rate$size_g
data.rate$rate2 = (data.rate$rate2)*(-1)
rate_umass <- subset(data.rate, X>=82)
rate_mhc <- subset(data.rate, X<=81)

Linear Models

Body Size

size_resp <- lm(rate~size_g, data=data.rate)
plot(allEffects(size_resp))

Anova(size_resp)

Anova Table (Type II tests)

Response: rate
           Sum Sq  Df F value  Pr(>F)    
size_g    0.03516   1  16.221 8.5e-05 ***
Residuals 0.36631 169                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Stable Temperature

Infection Treatment

um1 <-lm(rate2~infected, data=rate_umass)
plot(allEffects(um1))

Anova(um1)

Anova Table (Type II tests)

Response: rate2
             Sum Sq Df F value Pr(>F)
infected  0.0000017  1  0.0201 0.8875
Residuals 0.0073870 88               

Pathogen Load

um2 <-lm(rate2~Crithidia, data=subset(rate_umass, infected==1))
plot(allEffects(um2))

Anova(um2)

Anova Table (Type II tests)

Response: rate2
              Sum Sq Df F value   Pr(>F)   
Crithidia 0.00025138  1  7.5639 0.008677 **
Residuals 0.00142906 43                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Thermal Ramping

Infection Treatment

moho<-lm(rate2~infected*set_temp_C, data=rate_mhc)
plot(allEffects(moho))

Anova(moho)

Anova Table (Type II tests)

Response: rate2
                       Sum Sq Df F value Pr(>F)
infected            0.0000710  1  1.6878 0.1978
set_temp_C          0.0000498  1  1.1829 0.2802
infected:set_temp_C 0.0000487  1  1.1563 0.2856
Residuals           0.0032405 77               

Pathogen Load

moho2 <- lm(rate2~Crithidia+set_temp_C, data=subset(rate_mhc, infected==1))
plot(allEffects(moho2))

Anova(moho2)

Anova Table (Type II tests)

Response: rate2
               Sum Sq Df F value Pr(>F)
Crithidia  0.00000950  1  0.2543 0.6169
set_temp_C 0.00000089  1  0.0239 0.8778
Residuals  0.00145642 39