feeding models

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

behavior_23Feb23 <- read.csv("behavior_23Feb23.csv")

behavior <- behavior_23Feb23

Feeding Model 1

m1 <- glmer(formula = feeding~temp+day_of_experiment+infected+total_alive+(1|parent), family = "poisson", data = behavior)
plot(allEffects(m1)) 

vif(m1) #This checks for colinearity, we want the values below 2.5 
             temp day_of_experiment          infected       total_alive 
         1.038569          1.118786          1.029644          1.118236

Model 1.1: everything except random effect

m1.1 <- glm(formula = feeding~temp+day_of_experiment+infected+total_alive, family = "poisson", data = behavior)
plot(allEffects(m1.1))

vif(m1.1)
             temp day_of_experiment          infected       total_alive 
         1.093894          1.178129          1.090494          1.095426

Compare dAIC

ICtab(m1,m1.1)
     dAIC df
m1   0.0  6 
m1.1 4.9  5 
# keep random effect

Model 1.2: keep random effect, everything else except total_alive

m1.2 <- glmer(formula=feeding~temp+day_of_experiment+infected+(1|parent),family="poisson", data=behavior)
plot(allEffects(m1.2))

Model 1.3: keep everything except infected

m1.3 <- glmer(formula = feeding~temp+day_of_experiment+total_alive+(1|parent), family = "poisson", data = behavior)
plot(allEffects(m1.3))

Model 1.4: keep everything except day of experiment

m1.4 <- glmer(formula = feeding~temp+infected+total_alive+(1|parent), family = "poisson", data = behavior)
plot(allEffects(m1.4))

Check these with performance

#According to this m1.3 is the best overall, but m1.2 has the lowest AIC 
compare_performance(m1,m1.1,m1.2,m1.3,m1.4, rank=TRUE)
Warning: Following indices with missing values are not used for ranking:
  R2_conditional, R2_marginal, ICC, R2_Nagelkerke
# Comparison of Model Performance Indices

Name |    Model |  RMSE | Sigma | Score_log | Score_spherical | AIC weights | AICc weights | BIC weights | Performance-Score
----------------------------------------------------------------------------------------------------------------------------
m1.3 | glmerMod | 1.070 | 1.000 |    -1.233 |           0.023 |       0.385 |        0.387 |       0.772 |            86.76%
m1   | glmerMod | 1.066 | 1.000 |    -1.233 |           0.023 |       0.542 |        0.540 |       0.083 |            85.02%
m1.4 | glmerMod | 1.064 | 1.000 |    -1.236 |           0.023 |       0.024 |        0.024 |       0.048 |            49.58%
m1.2 | glmerMod | 1.083 | 1.000 |    -1.237 |           0.023 |       0.002 |        0.002 |       0.005 |            31.46%
m1.1 |      glm | 1.076 | 1.049 |    -1.241 |           0.023 |       0.046 |        0.046 |       0.092 |             9.14%
behaviorprop<- behavior_23Feb23

feedingprop = (behaviorprop$feeding/behaviorprop$total_alive)
behaviorprop<- cbind(behaviorprop,feedingprop)

m1.32 <- glmer(formula = feedingprop~temp+day_of_experiment+total_alive+(1|parent), family = "binomial", data = behaviorprop)
Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!
boundary (singular) fit: see help('isSingular')
m1.32
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: feedingprop ~ temp + day_of_experiment + total_alive + (1 | parent)
   Data: behaviorprop
      AIC       BIC    logLik  deviance  df.resid 
 307.7433  333.5163 -148.8716  297.7433      1275 
Random effects:
 Groups Name        Std.Dev.
 parent (Intercept) 0       
Number of obs: 1280, groups:  parent, 6
Fixed Effects:
      (Intercept)               temp  day_of_experiment        total_alive  
        -45.44305            1.15394            0.02816           -0.18938  
optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
plot(allEffects(m1.32))
Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Warning in eval(family$initialize): non-integer #successes in a binomial glm!

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

behaviorprop$notfeeding <- behaviorprop$total_alive-behaviorprop$feeding

m1.5 <- glmer(formula = cbind(feeding,notfeeding)~temp+day_of_experiment+infected+total_alive+(1|parent), family = "binomial", data = behaviorprop)
summary(m1.5)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(feeding, notfeeding) ~ temp + day_of_experiment + infected +  
    total_alive + (1 | parent)
   Data: behaviorprop

     AIC      BIC   logLik deviance df.resid 
  3246.2   3277.1  -1617.1   3234.2     1274 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2909 -0.6928 -0.3336  0.5829  8.7238 

Random effects:
 Groups Name        Variance Std.Dev.
 parent (Intercept) 0.04107  0.2027  
Number of obs: 1280, groups:  parent, 6

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -7.807860   0.421473 -18.525  < 2e-16 ***
temp               0.207411   0.008000  25.925  < 2e-16 ***
day_of_experiment -0.012638   0.003617  -3.494 0.000476 ***
infected           0.115561   0.066453   1.739 0.082037 .  
total_alive       -0.123430   0.025461  -4.848 1.25e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) temp   dy_f_x infctd
temp        -0.784                     
dy_f_xprmnt -0.302  0.082              
infected    -0.122  0.059  0.156       
total_alive -0.694  0.162  0.297 -0.008
plot(allEffects(m1.5))