mar 7th stuff for jvw

# Stuff for Jenny from 3/7 meeting
# I lost my original code because I am the worst 

library (tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.4.0      ✔ purrr   1.0.1 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.3.0      ✔ stringr 1.5.0 
✔ readr   2.1.3      ✔ forcats 0.5.2 
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

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 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.

# Original code for feeding effects
behavior_23Feb23 <- read.csv("behavior_23Feb23.csv")

behavior <- behavior_23Feb23

m01 <- glm(formula = feeding~temp+day_of_experiment+infected+total_alive, family = "poisson", data = behavior)
summary(m01)

Call:
glm(formula = feeding ~ temp + day_of_experiment + infected + 
    total_alive, family = "poisson", data = behavior)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5306  -0.9771  -0.5120   0.5452   3.6374  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -6.863806   0.341015 -20.128  < 2e-16 ***
temp               0.184612   0.006789  27.193  < 2e-16 ***
day_of_experiment -0.008101   0.003073  -2.636 0.008388 ** 
infected           0.152720   0.052340   2.918 0.003525 ** 
total_alive        0.072359   0.018847   3.839 0.000123 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2710.2  on 1279  degrees of freedom
Residual deviance: 1402.7  on 1275  degrees of freedom
  (728 observations deleted due to missingness)
AIC: 3187.4

Number of Fisher Scoring iterations: 6

m01

Call:  glm(formula = feeding ~ temp + day_of_experiment + infected + 
    total_alive, family = "poisson", data = behavior)

Coefficients:
      (Intercept)               temp  day_of_experiment           infected  
        -6.863806           0.184612          -0.008101           0.152720  
      total_alive  
         0.072359  

Degrees of Freedom: 1279 Total (i.e. Null);  1275 Residual
  (728 observations deleted due to missingness)
Null Deviance:      2710 
Residual Deviance: 1403     AIC: 3187

plot(allEffects(m01)) 

# Change behavior from count to prop. with total_alive 
# Start with feeding and go from there lol 
behaviorprop<- behavior_23Feb23
xtabs(formula=~feeding, data=behaviorprop)

feeding
  0   1   2   3   4   5   6   7   8   9 
600 268 168  86  80  42  20  11   4   1 

# Try adding new column and making it feeding/total alive for each replicate
feedingprop = (behaviorprop$feeding/behaviorprop$total_alive)
behaviorprop<- cbind(behaviorprop,feedingprop)
# HELL YEAH I THINK THIS WORKS 
# Now redo that model, the 50+ warnings cannot hurt me if I can't see them 
m02 <- glm(formula = feedingprop~temp+day_of_experiment+infected+total_alive, family = "binomial", data = behaviorprop)

Warning in eval(family$initialize): non-integer #successes in a binomial glm!

m02

Call:  glm(formula = feedingprop ~ temp + day_of_experiment + infected + 
    total_alive, family = "binomial", data = behaviorprop)

Coefficients:
      (Intercept)               temp  day_of_experiment           infected  
         -8.49353            0.22110           -0.01178            0.20148  
      total_alive  
         -0.10265  

Degrees of Freedom: 1279 Total (i.e. Null);  1275 Residual
  (728 observations deleted due to missingness)
Null Deviance:      373.5 
Residual Deviance: 189.5    AIC: 504.8

plot(allEffects(m02))

# Interaction of number of bees feeding infected vs uninfected across thermal gradient 
# Feeding y, temp x, one line for infected and one for uninfected 
ggplot(data=behavior, aes(x=temp, y=feeding, color=as.character(infected)))+
  geom_point()+
  geom_jitter()+
  geom_smooth(method="lm")+
  theme_classic()+
  scale_color_discrete(na.translate=FALSE)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 728 rows containing non-finite values (`stat_smooth()`).

Warning: Removed 728 rows containing missing values (`geom_point()`).
Removed 728 rows containing missing values (`geom_point()`).

ggplot(data=behaviorprop, aes(x=temp, y=feedingprop, color=as.character(infected)))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()+
  scale_color_discrete(na.translate=FALSE)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 728 rows containing non-finite values (`stat_smooth()`).

Warning: Removed 728 rows containing missing values (`geom_point()`).

# Combine incubating and feeding by temperature, two lines for infected/uninfected 
feedingplot<- ggplot(behavior, aes(x=temp, y=feeding, color=as.character(infected)))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()+
  scale_color_discrete(na.translate=FALSE)

incubatingplot<- ggplot(behavior, aes(x=temp, y=incubating, color=as.character(infected)))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()+
  scale_color_discrete(na.translate=FALSE)

library(patchwork)
(feedingplot+incubatingplot)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 728 rows containing non-finite values (`stat_smooth()`).

Warning: Removed 728 rows containing missing values (`geom_point()`).

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 728 rows containing non-finite values (`stat_smooth()`).
Removed 728 rows containing missing values (`geom_point()`).