stationary 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

behaviorprop<- behavior_23Feb23

Stationary Model

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

vif(m3)
             temp day_of_experiment          infected       total_alive 
         1.038561          1.158002          1.048192          1.140449

Model 3.1: everything except random effect

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

vif(m3.1)
             temp day_of_experiment          infected       total_alive 
         1.065553          1.166925          1.072500          1.111387

Assess whether or not to keep random effect with dAIC

ICtab(m3,m3.1)
     dAIC  df
m3     0.0 6 
m3.1 137.4 5 
# keep it

Model 3.2: everything except total_alive

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

Model 3.3: everything except infected

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

Model 3.4: everything except day_of_experiment

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

Performance check

compare_performance(m3,m3.1,m3.2,m3.3,m3.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
----------------------------------------------------------------------------------------------------------------------------
m3   | glmerMod | 2.136 | 1.000 |    -2.088 |           0.023 |       1.000 |        1.000 |       1.000 |            88.78%
m3.2 | glmerMod | 2.153 | 1.000 |    -2.112 |           0.023 |    8.34e-14 |     8.42e-14 |    1.10e-12 |            41.71%
m3.3 | glmerMod | 2.163 | 1.000 |    -2.099 |           0.023 |    2.12e-06 |     2.14e-06 |    2.79e-05 |            37.19%
m3.4 | glmerMod | 2.237 | 1.000 |    -2.174 |           0.023 |    2.59e-48 |     2.61e-48 |    3.40e-47 |            28.57%
m3.1 |      glm | 2.215 | 1.389 |    -2.153 |           0.023 |    1.43e-30 |     1.45e-30 |    1.89e-29 |            12.78%
# Model 3 (original model with everything) works best
stationaryprop = (behaviorprop$stationary/behaviorprop$total_alive)
behaviorprop<- cbind(behaviorprop,stationaryprop)

m3.0 <- glmer(formula = stationaryprop~temp+day_of_experiment+infected+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')
plot(allEffects(m3.0)) 
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!

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

behaviorprop$notstationary <- behaviorprop$total_alive-behaviorprop$stationary

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

     AIC      BIC   logLik deviance df.resid 
  5891.7   5922.6  -2939.8   5879.7     1274 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9717 -1.0662 -0.1844  0.9518  5.2412 

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

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        0.855742   0.332095   2.577  0.00997 ** 
temp              -0.051084   0.003628 -14.080  < 2e-16 ***
day_of_experiment -0.048572   0.003111 -15.612  < 2e-16 ***
infected           0.320460   0.046328   6.917 4.61e-12 ***
total_alive        0.020088   0.024605   0.816  0.41425    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) temp   dy_f_x infctd
temp        -0.487                     
dy_f_xprmnt -0.354  0.130              
infected    -0.137  0.050  0.196       
total_alive -0.796  0.180  0.320  0.046
plot(allEffects(m3.5))