Mixed models to assess the relationship between human disturbance and risk-taking behaviors
R codes of the mixed models used to assess the relationship between human disturbance and three risk-taking behaviors: sum of active defensive behaviors used, escape latency and re-emergence of the turtle after escaping.
Packages used in this file
library(Hmisc)
library(writexl)
library(ggplot2)
library(dplyr)
library(lmerTest)
library(optimx)
library(lme4)
library(PerformanceAnalytics)
library(effects)
library(MuMIn)
library(ggeffects)
library(splines)
library(glmtoolbox)
library(afex)
library(nloptr)
library(dfoptim)
library(psych)
library(ordinal)Setting of the directory
You need to set the directory by selecting the file where the datasets are located.
Load dataset in R
df <- read.table("Behaviors.txt", header = T)The dataset has 1135 entries corresponding to the number of trial made on the 730 turtle sampled. However, all three behaviors were not always recorded at each trial. Thus, the number of observations for each behavior is less than the number of entries.
Preparation of the data for the mixed models
Verification of the format of each data
We want to make sure that the data are in the good format for the analyses.
# Risk-taking behaviors that will be used as response variable in their respective mixed model.
df$active.defenses <- as.numeric(df$active.defenses)
df$escape <- as.numeric(df$escape)
df$log.escape <- log(df$escape + 1) # log(x+1) transformation of the escape latency to be able to use a Gaussian distribution in the mixed model.
df$emergence <- as.numeric(df$emergence)
# Predictor variables
df$Year <- as.factor(df$Year)
df$JulianDay <- as.numeric(df$JulianDay)
df$Hour.active.defenses <- as.numeric(df$Hour.active.defenses)
df$Hour.Platform <- as.numeric(df$Hour.Platform)
df$Order.trial <- as.numeric(df$Order.trial)
df$CL <- as.numeric(df$CL)
df$Wind <- as.numeric(df$Wind)
df$ID.Temperature <- as.numeric(df$ID.Temperature)
df$Distance.channel <- as.numeric(df$Distance.channel)
df$Vessels.activities.2019 <- as.numeric(df$Vessels.activities.2019)
df$House.200 <- as.numeric(df$House.200)
df$House.300 <- as.numeric(df$House.300)
df$House.400 <- as.numeric(df$House.400)
df$Urban.200 <- as.numeric(df$Urban.200)
df$Urban.600 <- as.numeric(df$Urban.600)
df$Urban.900 <- as.numeric(df$Urban.900)Standardization of the predictor variables
All continuous predictor variables were standardized (mean zero, unit variance) before model selection.
df$JulianDay.scaled <- scale(df$JulianDay, center=TRUE, scale=TRUE)
df$Hour.active.defenses.scaled <- scale(df$Hour.active.defenses, center=TRUE, scale=TRUE)
df$Hour.Platform.scaled <- scale(df$Hour.Platform, center=TRUE, scale=TRUE)
df$Order.trial.scaled <- scale(df$Order.trial, center=TRUE, scale=TRUE)
df$CL.scaled <- scale(df$CL, center=TRUE, scale=TRUE)
df$Wind.scaled <- scale(df$Wind, center=TRUE, scale=TRUE)
df$ID.Temperature.scaled <- scale(df$ID.Temperature, center=TRUE, scale=TRUE)
df$Distance.channel.scaled <- scale(df$Distance.channel, center=TRUE, scale=TRUE)
df$Vessels.activities.2019.scaled <- scale(df$Vessels.activities.2019, center=TRUE, scale=TRUE)
df$House.200.scaled <- scale(df$House.200, center=TRUE, scale=TRUE)
df$House.300.scaled <- scale(df$House.300, center=TRUE, scale=TRUE)
df$House.400.scaled <- scale(df$House.400, center=TRUE, scale=TRUE)
df$Urban.200.scaled <- scale(df$Urban.200, center=TRUE, scale=TRUE)
df$Urban.600.scaled <- scale(df$Urban.600, center=TRUE, scale=TRUE)
df$Urban.900.scaled <- scale(df$Urban.900, center=TRUE, scale=TRUE)Descriptive statistics of the preditor variables
Only applicable for the continuous preditor variables
df.description <- psych::describe(df)
df.description## vars n mean sd median trimmed mad
## Year* 1 1135 1.48 0.50 1.00 1.48 0.00
## Site* 2 1135 11.37 6.90 12.00 11.41 10.38
## JulianDay 3 1135 185.89 29.71 189.00 186.39 38.55
## Hour.active.defenses 4 962 1244.37 249.94 1235.00 1221.68 244.63
## ID 5 1135 379.73 207.78 387.00 381.34 257.97
## Order.trial 6 1135 1.60 0.97 1.00 1.40 0.00
## Sex* 7 1109 1.63 0.48 2.00 1.67 0.00
## CL 8 1135 138.24 18.69 140.00 139.12 16.31
## active.defenses 9 1117 1.08 0.92 1.00 1.01 1.48
## Hour.Platform 10 1103 1295.05 249.25 1251.00 1278.94 240.18
## Wind 11 948 1.91 1.81 2.00 1.68 1.48
## ID.Temperature 12 1097 22.94 4.69 23.30 23.11 4.89
## escape 13 1115 65.74 117.54 26.00 35.63 23.72
## emergence 14 1071 0.30 0.46 0.00 0.25 0.00
## Tank.Site* 15 1135 1.22 0.42 1.00 1.15 0.00
## Distance.channel 16 1135 744.17 1104.15 438.04 528.47 479.74
## Vessels.activities.2019 17 1135 14.99 7.46 10.79 13.85 1.30
## House.200 18 1135 0.97 1.28 1.00 0.70 1.48
## House.300 19 1135 2.91 4.33 1.00 1.89 1.48
## House.400 20 1135 6.41 8.41 2.00 4.82 2.97
## Urban.200 21 1135 0.05 0.14 0.00 0.01 0.00
## Urban.600 22 1135 0.14 0.23 0.04 0.09 0.06
## Urban.900 23 1135 0.16 0.23 0.05 0.11 0.04
## dummy 24 1135 1.00 0.00 1.00 1.00 0.00
## spatial.cor 25 1135 10.78 5.96 12.00 10.83 7.41
## log.escape 26 1115 3.40 1.16 3.30 3.32 1.08
## JulianDay.scaled 27 1135 0.00 1.00 0.10 0.02 1.30
## Hour.active.defenses.scaled 28 962 0.00 1.00 -0.04 -0.09 0.98
## Hour.Platform.scaled 29 1103 0.00 1.00 -0.18 -0.06 0.96
## Order.trial.scaled 30 1135 0.00 1.00 -0.62 -0.21 0.00
## CL.scaled 31 1135 0.00 1.00 0.09 0.05 0.87
## Wind.scaled 32 948 0.00 1.00 0.05 -0.13 0.82
## ID.Temperature.scaled 33 1097 0.00 1.00 0.08 0.04 1.04
## Distance.channel.scaled 34 1135 0.00 1.00 -0.28 -0.20 0.43
## Vessels.activities.2019.scaled 35 1135 0.00 1.00 -0.56 -0.15 0.17
## House.200.scaled 36 1135 0.00 1.00 0.02 -0.21 1.16
## House.300.scaled 37 1135 0.00 1.00 -0.44 -0.23 0.34
## House.400.scaled 38 1135 0.00 1.00 -0.52 -0.19 0.35
## Urban.200.scaled 39 1135 0.00 1.00 -0.34 -0.31 0.00
## Urban.600.scaled 40 1135 0.00 1.00 -0.46 -0.23 0.27
## Urban.900.scaled 41 1135 0.00 1.00 -0.47 -0.22 0.19
## min max range skew kurtosis se
## Year* 1.00 2.00 1.00 0.08 -2.00 0.01
## Site* 1.00 22.00 21.00 -0.05 -1.39 0.20
## JulianDay 128.00 234.00 106.00 -0.11 -1.13 0.88
## Hour.active.defenses 730.00 2020.00 1290.00 0.80 0.59 8.06
## ID 2.00 742.00 740.00 -0.11 -1.06 6.17
## Order.trial 1.00 7.00 6.00 1.81 3.42 0.03
## Sex* 1.00 2.00 1.00 -0.56 -1.69 0.01
## CL 30.00 188.00 158.00 -0.63 1.24 0.55
## active.defenses 0.00 4.00 4.00 0.48 -0.55 0.03
## Hour.Platform 730.00 2020.00 1290.00 0.54 -0.11 7.50
## Wind 0.00 6.00 6.00 0.79 -0.26 0.06
## ID.Temperature 9.50 34.00 24.50 -0.30 -0.30 0.14
## escape 1.00 600.00 599.00 3.34 11.15 3.52
## emergence 0.00 1.00 1.00 0.87 -1.24 0.01
## Tank.Site* 1.00 2.00 1.00 1.34 -0.20 0.01
## Distance.channel 15.13 6140.65 6125.52 3.39 13.33 32.77
## Vessels.activities.2019 8.74 34.07 25.33 1.17 -0.19 0.22
## House.200 0.00 5.00 5.00 1.63 2.06 0.04
## House.300 0.00 17.00 17.00 1.91 2.84 0.13
## House.400 0.00 28.00 28.00 1.39 0.57 0.25
## Urban.200 0.00 0.57 0.57 2.92 6.78 0.00
## Urban.600 0.00 0.78 0.78 1.86 1.96 0.01
## Urban.900 0.00 0.81 0.81 1.84 2.04 0.01
## dummy 1.00 1.00 0.00 NaN NaN 0.00
## spatial.cor 1.00 20.00 19.00 -0.07 -1.25 0.18
## log.escape 0.69 6.40 5.71 0.59 0.13 0.03
## JulianDay.scaled -1.95 1.62 3.57 -0.11 -1.13 0.03
## Hour.active.defenses.scaled -2.06 3.10 5.16 0.80 0.59 0.03
## Hour.Platform.scaled -2.27 2.91 5.18 0.54 -0.11 0.03
## Order.trial.scaled -0.62 5.55 6.17 1.81 3.42 0.03
## CL.scaled -5.79 2.66 8.46 -0.63 1.24 0.03
## Wind.scaled -1.06 2.26 3.32 0.79 -0.26 0.03
## ID.Temperature.scaled -2.87 2.36 5.22 -0.30 -0.30 0.03
## Distance.channel.scaled -0.66 4.89 5.55 3.39 13.33 0.03
## Vessels.activities.2019.scaled -0.84 2.56 3.39 1.17 -0.19 0.03
## House.200.scaled -0.76 3.15 3.91 1.63 2.06 0.03
## House.300.scaled -0.67 3.26 3.93 1.91 2.84 0.03
## House.400.scaled -0.76 2.57 3.33 1.39 0.57 0.03
## Urban.200.scaled -0.34 3.65 3.99 2.92 6.78 0.03
## Urban.600.scaled -0.64 2.80 3.44 1.86 1.96 0.03
## Urban.900.scaled -0.70 2.84 3.54 1.84 2.04 0.03Visualization of data distribution
This step is important to make sure that we will use the most adequate data distribution in our models and also to find possible errors or aberrant values in the predictor variables.
Response variables: Risk-taking behaviors
Sum of active defensive behaviors
ggplot(df, aes(x=active.defenses)) + geom_histogram() # histogram
sum(is.na(df$active.defenses)) # number of NAs## [1] 18dotchart(df$active.defenses, main = "Sum of active defensive behaviors") #distribution of the data
table(df$active.defenses) # Number of observations per value##
## 0 1 2 3 4
## 340 431 261 82 3Escape latency
without log(x+1) transformation
ggplot(df, aes(x=escape)) + geom_histogram() # histogram
sum(is.na(df$escape)) # number of NAs## [1] 20dotchart(df$escape, main = "Escape latency") #distribution of the data
table(df$escape) # Number of observations per value##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 12 16 24 31 32 35 24 21 30 32 34 20 19 29 14 19 27 19 22
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 21 21 15 20 13 26 13 15 19 15 12 18 6 13 12 8 12 10 7 3
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 6 11 4 6 4 13 8 9 3 8 5 2 9 4 9 4 3 7 2 5
## 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
## 3 2 4 3 3 2 1 3 3 4 6 3 4 1 4 1 3 2 4 1
## 82 83 84 85 86 88 89 90 91 92 93 95 96 97 98 99 100 101 107 108
## 4 1 1 3 3 2 3 1 3 3 2 2 2 1 3 1 4 3 1 2
## 110 111 114 115 116 119 120 121 122 124 125 126 127 128 129 131 135 138 142 143
## 5 2 2 1 1 2 1 1 1 2 1 1 2 2 3 3 3 3 1 1
## 145 148 149 150 152 159 162 164 170 171 172 175 177 180 181 182 188 189 192 194
## 1 1 2 1 2 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1
## 204 206 207 211 214 216 220 222 225 230 231 232 233 238 250 251 252 266 269 272
## 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1
## 282 283 286 289 291 295 305 308 309 311 321 324 328 330 331 341 348 350 367 369
## 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1
## 378 410 411 437 438 459 471 472 479 483 510 512 542 545 556 584 600
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 29with log(x+1) transformation
ggplot(df, aes(x=log.escape)) + geom_histogram() # histogram
sum(is.na(df$log.escape)) # number of NAs## [1] 20dotchart(df$log.escape, main = "Escape latency") #distribution of the data
table(df$log.escape) # Number of observations per value##
## 0.693147180559945 1.09861228866811 1.38629436111989 1.6094379124341
## 1 12 16 24
## 1.79175946922805 1.94591014905531 2.07944154167984 2.19722457733622
## 31 32 35 24
## 2.30258509299405 2.39789527279837 2.484906649788 2.56494935746154
## 21 30 32 34
## 2.63905732961526 2.70805020110221 2.77258872223978 2.83321334405622
## 20 19 29 14
## 2.89037175789616 2.94443897916644 2.99573227355399 3.04452243772342
## 19 27 19 22
## 3.09104245335832 3.13549421592915 3.17805383034795 3.2188758248682
## 21 21 15 20
## 3.25809653802148 3.29583686600433 3.3322045101752 3.36729582998647
## 13 26 13 15
## 3.40119738166216 3.43398720448515 3.46573590279973 3.49650756146648
## 19 15 12 18
## 3.52636052461616 3.55534806148941 3.58351893845611 3.61091791264422
## 6 13 12 8
## 3.63758615972639 3.66356164612965 3.68887945411394 3.71357206670431
## 12 10 7 3
## 3.73766961828337 3.76120011569356 3.78418963391826 3.80666248977032
## 6 11 4 6
## 3.8286413964891 3.85014760171006 3.87120101090789 3.89182029811063
## 4 13 8 9
## 3.91202300542815 3.93182563272433 3.95124371858143 3.97029191355212
## 3 8 5 2
## 3.98898404656427 4.00733318523247 4.02535169073515 4.04305126783455
## 9 4 9 4
## 4.06044301054642 4.07753744390572 4.0943445622221 4.11087386417331
## 3 7 2 5
## 4.14313472639153 4.15888308335967 4.17438726989564 4.18965474202643
## 3 2 4 3
## 4.20469261939097 4.21950770517611 4.23410650459726 4.24849524204936
## 3 2 1 3
## 4.26267987704132 4.27666611901606 4.29045944114839 4.30406509320417
## 3 4 6 3
## 4.31748811353631 4.33073334028633 4.34380542185368 4.35670882668959
## 4 1 4 1
## 4.36944785246702 4.38202663467388 4.39444915467244 4.40671924726425
## 3 2 4 1
## 4.4188406077966 4.43081679884331 4.44265125649032 4.45434729625351
## 4 1 1 3
## 4.46590811865458 4.48863636973214 4.49980967033027 4.51085950651685
## 3 2 3 1
## 4.52178857704904 4.53259949315326 4.54329478227 4.56434819146784
## 3 3 2 2
## 4.57471097850338 4.58496747867057 4.59511985013459 4.60517018598809
## 2 1 3 1
## 4.61512051684126 4.62497281328427 4.68213122712422 4.69134788222914
## 4 3 1 2
## 4.70953020131233 4.71849887129509 4.74493212836325 4.75359019110636
## 5 2 2 1
## 4.76217393479776 4.78749174278205 4.79579054559674 4.80402104473326
## 1 2 1 1
## 4.81218435537242 4.8283137373023 4.83628190695148 4.84418708645859
## 1 2 1 1
## 4.85203026391962 4.85981240436167 4.86753445045558 4.88280192258637
## 2 2 3 3
## 4.91265488573605 4.93447393313069 4.96284463025991 4.969813299576
## 3 3 1 1
## 4.98360662170834 5.00394630594546 5.01063529409626 5.01727983681492
## 1 1 2 1
## 5.03043792139244 5.07517381523383 5.09375020080676 5.10594547390058
## 2 1 1 1
## 5.14166355650266 5.14749447681345 5.15329159449778 5.17048399503815
## 1 1 2 2
## 5.18178355029209 5.19849703126583 5.2040066870768 5.20948615284142
## 1 1 1 1
## 5.24174701505964 5.24702407216049 5.26269018890489 5.27299955856375
## 1 1 1 1
## 5.32300997913841 5.33271879326537 5.33753807970132 5.35658627467201
## 1 1 1 1
## 5.37063802812766 5.37989735354046 5.39816270151775 5.40717177146012
## 1 1 1 1
## 5.42053499927229 5.44241771052179 5.44673737166631 5.4510384535657
## 1 1 1 2
## 5.4553211153577 5.47646355193151 5.52545293913178 5.52942908751142
## 1 1 1 1
## 5.53338948872752 5.58724865840025 5.59842195899838 5.60947179518496
## 1 2 1 1
## 5.64544689764324 5.64897423816121 5.65948221575962 5.66988092298052
## 1 1 1 1
## 5.67675380226828 5.69035945432406 5.72358510195238 5.73334127689775
## 1 1 1 1
## 5.73657229747919 5.74300318780948 5.77455154554441 5.78382518232974
## 1 2 1 2
## 5.79605775076537 5.80211837537706 5.80513496891649 5.8348107370626
## 1 1 1 1
## 5.85507192220243 5.86078622346587 5.90808293816893 5.91350300563827
## 1 1 1 1
## 5.93753620508243 6.01859321449623 6.02102334934953 6.08221891037645
## 1 1 1 1
## 6.08449941307517 6.13122648948314 6.15697898558556 6.15909538849193
## 1 1 1 1
## 6.17378610390194 6.18208490671663 6.2363695902037 6.24027584517077
## 1 1 1 1
## 6.29710931993394 6.30261897574491 6.32256523992728 6.37161184723186
## 1 1 1 1
## 6.39859493453521
## 29Emergence of the turtle after escaping
ggplot(df, aes(x=emergence)) + geom_histogram() # histogram
sum(is.na(df$emergence)) # number of NAs## [1] 64dotchart(df$emergence, main = "Emergence of the turtle after escaping") #distribution of the data
table(df$emergence) # Number of observations per value##
## 0 1
## 750 321Predictor variables
Here, we will give you one example, but the visualization was made on all the predictor variables
Shortest aquatic distance to the navigation channel
ggplot(df, aes(x=Distance.channel)) + geom_histogram() # histogram## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
sum(is.na(df$Distance.channel)) # number of NAs## [1] 0dotchart(df$Distance.channel, main = "Shortest distance to the navigation channel") #distribution of the data
table(df$Distance.channel) # Number of observation per value##
## 15.13313295 59.54172954 96.26934503 98.2983472 99.38102454 100.1995941
## 63 15 36 25 50 78
## 114.46075 175.433099 175.9749596 192.048037 202.7871069 316.8369367
## 30 40 51 32 37 46
## 388.8580992 438.0430606 479.398553 516.2279418 639.9552277 656.2130671
## 62 74 59 40 42 74
## 783.367761 1557.737778 2006.62065 6140.654572
## 54 77 118 32Verification of model assumptions
Sum of active defensive behaviors
Multicollinearity: Generalized variance inflation factor (GVIF^(1/(2*df)))
GVIF^(1/(2*df)) over two indicate the presence of multicollinearity
mod.active.defenses.vif <- lm(active.defenses ~ House.200.scaled + Urban.900.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled, data = df, na.action=na.exclude)
gvif(mod.active.defenses.vif)## GVIF df GVIF^(1/(2*df))
## House.200.scaled 3.4036 1 1.8449
## Urban.900.scaled 4.3233 1 2.0793
## Distance.channel.scaled 1.5719 1 1.2538
## Vessels.activities.2019.scaled 1.6432 1 1.2819
## Order.trial.scaled 1.2413 1 1.1141
## Sex 1.2647 1 1.1246
## Tank.Site 2.2465 1 1.4988
## Year 1.4893 1 1.2204
## JulianDay.scaled 2.8594 1 1.6910
## Hour.active.defenses.scaled 1.6055 1 1.2671
## ID.Temperature.scaled 2.1006 1 1.4493
## CL.scaled 1.2635 1 1.1241GVIF^(1/(2*df)) of Urban.900.scaled is over two. We need to delete it and recalculate the factors.
mod.active.defenses.vif.2 <- lm(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled, data = df, na.action=na.exclude)
gvif(mod.active.defenses.vif.2)## GVIF df GVIF^(1/(2*df))
## House.200.scaled 1.5730 1 1.2542
## Distance.channel.scaled 1.5632 1 1.2503
## Vessels.activities.2019.scaled 1.5994 1 1.2647
## Order.trial.scaled 1.1896 1 1.0907
## Sex 1.2139 1 1.1018
## Tank.Site 2.2362 1 1.4954
## Year 1.4890 1 1.2202
## JulianDay.scaled 2.6970 1 1.6423
## Hour.active.defenses.scaled 1.6031 1 1.2661
## ID.Temperature.scaled 2.0675 1 1.4379
## CL.scaled 1.2606 1 1.1228Calculation of Pearson and Spearman correlation coefficients
Pearson correlation coefficients
# Visualization of the correlations
cor.pearson.active.defenses <- df[, c(18,23,6,3,4,12,8,16,17)]
chart.Correlation(cor.pearson.active.defenses, histogram=TRUE, pch=19)
# Creation of the correlation table
table.corr.pearson.active.defenses <- rcorr(as.matrix(df[,c(18,23,6,3,4,12,8,16,17)]), type="pearson")
table.cor.pearson.active.defenses.r <- table.corr.pearson.active.defenses$r # Pearson correlation coefficients
table.cor.pearson.active.defenses.p <- table.corr.pearson.active.defenses$P # P values of the correlations
table.cor.pearson.active.defenses.r## House.200 Urban.900 Order.trial JulianDay
## House.200 1.0000000 0.79775182 0.14441075 -0.4024391
## Urban.900 0.7977518 1.00000000 0.09768635 -0.5580418
## Order.trial 0.1444107 0.09768635 1.00000000 0.1637192
## JulianDay -0.4024391 -0.55804179 0.16371919 1.0000000
## Hour.active.defenses -0.2829645 -0.28477450 -0.09146472 0.3068386
## ID.Temperature -0.1886545 -0.34461657 0.04136506 0.4406745
## CL 0.1519129 0.22519047 0.03401162 -0.1720253
## Distance.channel -0.1016703 -0.22397571 0.09912351 0.1010733
## Vessels.activities.2019 -0.1076827 -0.34645587 0.01856539 0.3303980
## Hour.active.defenses ID.Temperature CL
## House.200 -0.282964533 -0.188654524 0.15191293
## Urban.900 -0.284774503 -0.344616571 0.22519047
## Order.trial -0.091464721 0.041365059 0.03401162
## JulianDay 0.306838617 0.440674474 -0.17202529
## Hour.active.defenses 1.000000000 -0.007616587 -0.13533410
## ID.Temperature -0.007616587 1.000000000 -0.06344290
## CL -0.135334096 -0.063442901 1.00000000
## Distance.channel 0.012363808 0.253073722 -0.01469142
## Vessels.activities.2019 0.127469693 0.318224913 -0.08012667
## Distance.channel Vessels.activities.2019
## House.200 -0.10167029 -0.10768266
## Urban.900 -0.22397571 -0.34645587
## Order.trial 0.09912351 0.01856539
## JulianDay 0.10107332 0.33039803
## Hour.active.defenses 0.01236381 0.12746969
## ID.Temperature 0.25307372 0.31822491
## CL -0.01469142 -0.08012667
## Distance.channel 1.00000000 0.48906628
## Vessels.activities.2019 0.48906628 1.00000000table.cor.pearson.active.defenses.p## House.200 Urban.900 Order.trial JulianDay
## House.200 NA 0.000000e+00 1.032023e-06 0.000000e+00
## Urban.900 0.000000e+00 NA 9.832099e-04 0.000000e+00
## Order.trial 1.032023e-06 9.832099e-04 NA 2.905276e-08
## JulianDay 0.000000e+00 0.000000e+00 2.905276e-08 NA
## Hour.active.defenses 0.000000e+00 0.000000e+00 4.523448e-03 0.000000e+00
## ID.Temperature 3.016367e-10 0.000000e+00 1.709756e-01 0.000000e+00
## CL 2.716231e-07 1.643130e-14 2.522456e-01 5.455805e-09
## Distance.channel 6.027679e-04 2.264855e-14 8.257794e-04 6.493434e-04
## Vessels.activities.2019 2.786523e-04 0.000000e+00 5.320837e-01 0.000000e+00
## Hour.active.defenses ID.Temperature CL
## House.200 0.000000e+00 3.016367e-10 2.716231e-07
## Urban.900 0.000000e+00 0.000000e+00 1.643130e-14
## Order.trial 4.523448e-03 1.709756e-01 2.522456e-01
## JulianDay 0.000000e+00 0.000000e+00 5.455805e-09
## Hour.active.defenses NA 8.171484e-01 2.537130e-05
## ID.Temperature 8.171484e-01 NA 3.563983e-02
## CL 2.537130e-05 3.563983e-02 NA
## Distance.channel 7.017246e-01 0.000000e+00 6.210013e-01
## Vessels.activities.2019 7.350096e-05 0.000000e+00 6.916929e-03
## Distance.channel Vessels.activities.2019
## House.200 6.027679e-04 2.786523e-04
## Urban.900 2.264855e-14 0.000000e+00
## Order.trial 8.257794e-04 5.320837e-01
## JulianDay 6.493434e-04 0.000000e+00
## Hour.active.defenses 7.017246e-01 7.350096e-05
## ID.Temperature 0.000000e+00 0.000000e+00
## CL 6.210013e-01 6.916929e-03
## Distance.channel NA 0.000000e+00
## Vessels.activities.2019 0.000000e+00 NASpearman correlation coefficients
# Creation of the correlation table
table.corr.spearman.active.defenses <- rcorr(as.matrix(df[,c(18,23,6,3,4,12,8,16,17)]), type="spearman")
table.cor.spearman.active.defenses.r <- table.corr.spearman.active.defenses$r # Pearson correlation coefficients
table.cor.spearman.active.defenses.p <- table.corr.spearman.active.defenses$P # P values of the correlations
table.cor.spearman.active.defenses.r## House.200 Urban.900 Order.trial JulianDay
## House.200 1.00000000 0.4535131 0.06902811 -0.3156439
## Urban.900 0.45351306 1.0000000 -0.12653876 -0.7159501
## Order.trial 0.06902811 -0.1265388 1.00000000 0.2317717
## JulianDay -0.31564388 -0.7159501 0.23177171 1.0000000
## Hour.active.defenses -0.25020535 -0.2947812 -0.10794554 0.2546564
## ID.Temperature -0.09263193 -0.4414567 0.01087592 0.3952979
## CL 0.08753173 0.1558850 0.02353021 -0.1687139
## Distance.channel -0.20074713 -0.5183493 0.24200697 0.5566049
## Vessels.activities.2019 0.03760662 -0.3465167 -0.02161513 0.1889507
## Hour.active.defenses ID.Temperature CL
## House.200 -0.25020535 -0.09263193 0.08753173
## Urban.900 -0.29478119 -0.44145674 0.15588500
## Order.trial -0.10794554 0.01087592 0.02353021
## JulianDay 0.25465640 0.39529790 -0.16871392
## Hour.active.defenses 1.00000000 0.06180481 -0.14279438
## ID.Temperature 0.06180481 1.00000000 -0.06566684
## CL -0.14279438 -0.06566684 1.00000000
## Distance.channel 0.28241437 0.38823857 -0.11317974
## Vessels.activities.2019 0.04809711 0.19316565 -0.12635136
## Distance.channel Vessels.activities.2019
## House.200 -0.2007471 0.03760662
## Urban.900 -0.5183493 -0.34651666
## Order.trial 0.2420070 -0.02161513
## JulianDay 0.5566049 0.18895071
## Hour.active.defenses 0.2824144 0.04809711
## ID.Temperature 0.3882386 0.19316565
## CL -0.1131797 -0.12635136
## Distance.channel 1.0000000 0.36891367
## Vessels.activities.2019 0.3689137 1.00000000table.cor.spearman.active.defenses.p## House.200 Urban.900 Order.trial JulianDay
## House.200 NA 0.000000e+00 2.003175e-02 0.000000e+00
## Urban.900 0.000000e+00 NA 1.906574e-05 0.000000e+00
## Order.trial 2.003175e-02 1.906574e-05 NA 2.664535e-15
## JulianDay 0.000000e+00 0.000000e+00 2.664535e-15 NA
## Hour.active.defenses 3.330669e-15 0.000000e+00 7.979565e-04 8.881784e-16
## ID.Temperature 2.132226e-03 0.000000e+00 7.189809e-01 0.000000e+00
## CL 3.163985e-03 1.304308e-07 4.283820e-01 1.073361e-08
## Distance.channel 8.776313e-12 0.000000e+00 0.000000e+00 0.000000e+00
## Vessels.activities.2019 2.055095e-01 0.000000e+00 4.669252e-01 1.393752e-10
## Hour.active.defenses ID.Temperature CL
## House.200 3.330669e-15 2.132226e-03 3.163985e-03
## Urban.900 0.000000e+00 0.000000e+00 1.304308e-07
## Order.trial 7.979565e-04 7.189809e-01 4.283820e-01
## JulianDay 8.881784e-16 0.000000e+00 1.073361e-08
## Hour.active.defenses NA 6.038717e-02 8.749538e-06
## ID.Temperature 6.038717e-02 NA 2.964369e-02
## CL 8.749538e-06 2.964369e-02 NA
## Distance.channel 0.000000e+00 0.000000e+00 1.328682e-04
## Vessels.activities.2019 1.360378e-01 1.108240e-10 1.961956e-05
## Distance.channel Vessels.activities.2019
## House.200 8.776313e-12 2.055095e-01
## Urban.900 0.000000e+00 0.000000e+00
## Order.trial 0.000000e+00 4.669252e-01
## JulianDay 0.000000e+00 1.393752e-10
## Hour.active.defenses 0.000000e+00 1.360378e-01
## ID.Temperature 0.000000e+00 1.108240e-10
## CL 1.328682e-04 1.961956e-05
## Distance.channel NA 0.000000e+00
## Vessels.activities.2019 0.000000e+00 NAThe Pearson correlation coefficient of 0.8 between House.200 and Urban.900 confirm the deletion of Urban.900 with the calculation of the GVIF^(1/(2*df))
Verification of the assumptions with the initial model
mod.ad.full <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|Site) + (1|ID), data = df, REML = TRUE, na.action=na.exclude)
summary(mod.ad.full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled +
## ID.Temperature.scaled + CL.scaled + (1 | Site) + (1 | ID)
## Data: df
##
## REML criterion at convergence: 2263.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.47741 -0.53022 -0.06241 0.47665 2.84307
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.44852 0.6697
## Site (Intercept) 0.04453 0.2110
## Residual 0.32432 0.5695
## Number of obs: 889, groups: ID, 701; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.82717 0.08462 52.60258 9.775 2.01e-13
## House.200.scaled -0.13595 0.06086 15.86951 -2.234 0.04026
## Distance.channel.scaled -0.15182 0.05886 15.72871 -2.579 0.02037
## Vessels.activities.2019.scaled 0.20682 0.06839 12.89670 3.024 0.00985
## Order.trial.scaled 0.01857 0.03329 356.16926 0.558 0.57720
## SexM 0.29936 0.07440 673.21237 4.024 6.37e-05
## Tank.SiteT -0.19781 0.12592 477.39445 -1.571 0.11686
## Year2020 0.05110 0.08049 382.27678 0.635 0.52585
## JulianDay.scaled -0.01245 0.06886 45.61071 -0.181 0.85735
## Hour.active.defenses.scaled 0.01094 0.03737 563.75470 0.293 0.76971
## ID.Temperature.scaled 0.00352 0.04598 317.05133 0.077 0.93902
## CL.scaled 0.11759 0.03982 641.96377 2.953 0.00326
##
## (Intercept) ***
## House.200.scaled *
## Distance.channel.scaled *
## Vessels.activities.2019.scaled **
## Order.trial.scaled
## SexM ***
## Tank.SiteT
## Year2020
## JulianDay.scaled
## Hour.active.defenses.scaled
## ID.Temperature.scaled
## CL.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020 JlnDy.
## Hs.200.scld 0.080
## Dstnc.chnn. 0.142 0.115
## Vssl..2019. -0.056 -0.060 -0.483
## Ordr.trl.sc 0.179 -0.034 -0.023 0.068
## SexM -0.570 0.037 -0.016 0.005 -0.042
## Tank.SiteT -0.145 0.006 -0.143 0.042 0.080 -0.040
## Year2020 -0.385 -0.193 -0.149 0.090 -0.131 0.041 -0.249
## JulnDy.scld 0.122 0.248 0.253 -0.300 -0.085 0.065 -0.401 0.062
## Hr.ctv.dfn. 0.096 0.061 0.075 -0.077 0.138 0.014 -0.542 0.178 0.149
## ID.Tmprtr.s 0.018 0.091 -0.107 -0.049 -0.002 0.009 0.367 -0.295 -0.427
## CL.scaled -0.261 -0.034 -0.035 0.009 -0.043 0.378 -0.061 0.070 0.090
## Hr.c.. ID.Tm.
## Hs.200.scld
## Dstnc.chnn.
## Vssl..2019.
## Ordr.trl.sc
## SexM
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.ctv.dfn.
## ID.Tmprtr.s -0.029
## CL.scaled 0.061 -0.036plot(mod.ad.full) 
hist(residuals(mod.ad.full)) 
qqnorm(resid(mod.ad.full))
Escape latency
Multicollinearity: Generalized variance inflation factor (GVIF^(1/(2*df)))
GVIF^(1/(2*df)) over two indicate the presence of multicollinearity
mod.escape.vif <- lm(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled, data = df, na.action=na.exclude)
gvif(mod.escape.vif)## GVIF df GVIF^(1/(2*df))
## House.400.scaled 1.8537 1 1.3615
## Urban.200.scaled 1.6299 1 1.2767
## Distance.channel.scaled 1.9686 1 1.4031
## Vessels.activities.2019.scaled 1.9477 1 1.3956
## Order.trial.scaled 1.2602 1 1.1226
## Sex 1.1486 1 1.0717
## Tank.Site 3.3981 1 1.8434
## Year 1.2799 1 1.1313
## JulianDay.scaled 2.2605 1 1.5035
## Hour.Platform.scaled 1.7462 1 1.3214
## ID.Temperature.scaled 1.3960 1 1.1815
## CL.scaled 1.2457 1 1.1161
## Wind.scaled 1.9035 1 1.3797All the GVIF^(1/(2*df)) are under two.
Calculation of Pearson and Spearman correlation coefficients
Pearson correlation coefficients
# Visualization of the correlations
cor.pearson.escape <- df[, c(20,21,6,3,10,12,8,16,17,11)]
chart.Correlation(cor.pearson.escape, histogram=TRUE, pch=19)
# Creation of the correlation table
table.corr.pearson.escape <- rcorr(as.matrix(df[,c(20,21,6,3,10,12,8,16,17,11)]), type="pearson")
table.cor.pearson.escape.r <- table.corr.pearson.escape$r # Pearson correlation coefficients
table.cor.pearson.escape.p <- table.corr.pearson.escape$P # P values of the correlations
table.cor.pearson.escape.r## House.400 Urban.200 Order.trial JulianDay
## House.400 1.000000000 0.4326890387 -0.009439911 -0.3542009
## Urban.200 0.432689039 1.0000000000 0.197608295 -0.3859154
## Order.trial -0.009439911 0.1976082954 1.000000000 0.1637192
## JulianDay -0.354200905 -0.3859154480 0.163719193 1.0000000
## Hour.Platform -0.124917325 -0.1546110815 -0.105384108 0.2253572
## ID.Temperature -0.155963844 -0.2236247483 0.041365059 0.4406745
## CL 0.167195575 0.1848282077 0.034011619 -0.1720253
## Distance.channel -0.253479294 -0.0588967867 0.099123509 0.1010733
## Vessels.activities.2019 0.022988129 -0.2109359759 0.018565394 0.3303980
## Wind 0.012066735 -0.0004708797 0.026545568 -0.2562203
## Hour.Platform ID.Temperature CL
## House.400 -0.12491732 -0.15596384 0.16719557
## Urban.200 -0.15461108 -0.22362475 0.18482821
## Order.trial -0.10538411 0.04136506 0.03401162
## JulianDay 0.22535717 0.44067447 -0.17202529
## Hour.Platform 1.00000000 0.03449218 -0.11957931
## ID.Temperature 0.03449218 1.00000000 -0.06344290
## CL -0.11957931 -0.06344290 1.00000000
## Distance.channel 0.10732729 0.25307372 -0.01469142
## Vessels.activities.2019 0.09445235 0.31822491 -0.08012667
## Wind -0.34038488 0.02051645 0.07352587
## Distance.channel Vessels.activities.2019 Wind
## House.400 -0.25347929 0.02298813 0.0120667353
## Urban.200 -0.05889679 -0.21093598 -0.0004708797
## Order.trial 0.09912351 0.01856539 0.0265455684
## JulianDay 0.10107332 0.33039803 -0.2562202731
## Hour.Platform 0.10732729 0.09445235 -0.3403848794
## ID.Temperature 0.25307372 0.31822491 0.0205164543
## CL -0.01469142 -0.08012667 0.0735258655
## Distance.channel 1.00000000 0.48906628 -0.0292223219
## Vessels.activities.2019 0.48906628 1.00000000 -0.1317897600
## Wind -0.02922232 -0.13178976 1.0000000000table.cor.pearson.escape.p## House.400 Urban.200 Order.trial JulianDay
## House.400 NA 0.000000e+00 7.507239e-01 0.000000e+00
## Urban.200 0.000000e+00 NA 1.862999e-11 0.000000e+00
## Order.trial 7.507239e-01 1.862999e-11 NA 2.905276e-08
## JulianDay 0.000000e+00 0.000000e+00 2.905276e-08 NA
## Hour.Platform 3.178493e-05 2.468277e-07 4.553125e-04 3.641532e-14
## ID.Temperature 2.085079e-07 6.705747e-14 1.709756e-01 0.000000e+00
## CL 1.457405e-08 3.516900e-10 2.522456e-01 5.455805e-09
## Distance.channel 0.000000e+00 4.728322e-02 8.257794e-04 6.493434e-04
## Vessels.activities.2019 4.391003e-01 6.998846e-13 5.320837e-01 0.000000e+00
## Wind 7.105978e-01 9.884478e-01 4.142753e-01 1.110223e-15
## Hour.Platform ID.Temperature CL
## House.400 3.178493e-05 2.085079e-07 1.457405e-08
## Urban.200 2.468277e-07 6.705747e-14 3.516900e-10
## Order.trial 4.553125e-04 1.709756e-01 2.522456e-01
## JulianDay 3.641532e-14 0.000000e+00 5.455805e-09
## Hour.Platform NA 2.602924e-01 6.857072e-05
## ID.Temperature 2.602924e-01 NA 3.563983e-02
## CL 6.857072e-05 3.563983e-02 NA
## Distance.channel 3.559317e-04 0.000000e+00 6.210013e-01
## Vessels.activities.2019 1.687374e-03 0.000000e+00 6.916929e-03
## Wind 0.000000e+00 5.333715e-01 2.357961e-02
## Distance.channel Vessels.activities.2019 Wind
## House.400 0.0000000000 4.391003e-01 7.105978e-01
## Urban.200 0.0472832150 6.998846e-13 9.884478e-01
## Order.trial 0.0008257794 5.320837e-01 4.142753e-01
## JulianDay 0.0006493434 0.000000e+00 1.110223e-15
## Hour.Platform 0.0003559317 1.687374e-03 0.000000e+00
## ID.Temperature 0.0000000000 0.000000e+00 5.333715e-01
## CL 0.6210013281 6.916929e-03 2.357961e-02
## Distance.channel NA 0.000000e+00 3.687865e-01
## Vessels.activities.2019 0.0000000000 NA 4.697953e-05
## Wind 0.3687865313 4.697953e-05 NASpearman correlation coefficients
# Creation of the correlation table
table.corr.spearman.escape <- rcorr(as.matrix(df[,c(20,21,6,3,10,12,8,16,17,11)]), type="spearman")
table.cor.spearman.escape.r <- table.corr.spearman.escape$r # Pearson correlation coefficients
table.cor.spearman.escape.p <- table.corr.spearman.escape$P # P values of the correlations
table.cor.spearman.escape.r## House.400 Urban.200 Order.trial JulianDay
## House.400 1.000000000 0.45335084 0.005982156 -0.2436402
## Urban.200 0.453350837 1.00000000 0.017583899 -0.4426339
## Order.trial 0.005982156 0.01758390 1.000000000 0.2317717
## JulianDay -0.243640190 -0.44263394 0.231771710 1.0000000
## Hour.Platform -0.145961882 -0.05526323 -0.090744731 0.1850147
## ID.Temperature -0.058602842 -0.19567629 0.010875919 0.3952979
## CL 0.130176306 0.13710741 0.023530214 -0.1687139
## Distance.channel -0.277569869 -0.10029872 0.242006973 0.5566049
## Vessels.activities.2019 0.131077533 -0.30491389 -0.021615128 0.1889507
## Wind -0.016490187 0.06133348 -0.033315669 -0.2682615
## Hour.Platform ID.Temperature CL
## House.400 -0.145961882 -0.05860284 0.13017631
## Urban.200 -0.055263225 -0.19567629 0.13710741
## Order.trial -0.090744731 0.01087592 0.02353021
## JulianDay 0.185014667 0.39529790 -0.16871392
## Hour.Platform 1.000000000 0.08971374 -0.11614837
## ID.Temperature 0.089713738 1.00000000 -0.06566684
## CL -0.116148369 -0.06566684 1.00000000
## Distance.channel 0.277712076 0.38823857 -0.11317974
## Vessels.activities.2019 -0.007978281 0.19316565 -0.12635136
## Wind -0.412698498 0.10203512 0.08813909
## Distance.channel Vessels.activities.2019 Wind
## House.400 -0.2775699 0.131077533 -0.01649019
## Urban.200 -0.1002987 -0.304913889 0.06133348
## Order.trial 0.2420070 -0.021615128 -0.03331567
## JulianDay 0.5566049 0.188950707 -0.26826152
## Hour.Platform 0.2777121 -0.007978281 -0.41269850
## ID.Temperature 0.3882386 0.193165645 0.10203512
## CL -0.1131797 -0.126351359 0.08813909
## Distance.channel 1.0000000 0.368913666 -0.27032411
## Vessels.activities.2019 0.3689137 1.000000000 -0.08976149
## Wind -0.2703241 -0.089761492 1.00000000table.cor.spearman.escape.p## House.400 Urban.200 Order.trial JulianDay
## House.400 NA 0.000000e+00 8.404506e-01 0.000000e+00
## Urban.200 0.000000e+00 NA 5.539910e-01 0.000000e+00
## Order.trial 8.404506e-01 5.539910e-01 NA 2.664535e-15
## JulianDay 0.000000e+00 0.000000e+00 2.664535e-15 NA
## Hour.Platform 1.126180e-06 6.655125e-02 2.556506e-03 5.972247e-10
## ID.Temperature 5.232635e-02 6.280887e-11 7.189809e-01 0.000000e+00
## CL 1.084961e-05 3.553267e-06 4.283820e-01 1.073361e-08
## Distance.channel 0.000000e+00 7.147639e-04 0.000000e+00 0.000000e+00
## Vessels.activities.2019 9.413183e-06 0.000000e+00 4.669252e-01 1.393752e-10
## Wind 6.120910e-01 5.906359e-02 3.055012e-01 0.000000e+00
## Hour.Platform ID.Temperature CL
## House.400 1.126180e-06 5.232635e-02 1.084961e-05
## Urban.200 6.655125e-02 6.280887e-11 3.553267e-06
## Order.trial 2.556506e-03 7.189809e-01 4.283820e-01
## JulianDay 5.972247e-10 0.000000e+00 1.073361e-08
## Hour.Platform NA 3.357286e-03 1.105565e-04
## ID.Temperature 3.357286e-03 NA 2.964369e-02
## CL 1.105565e-04 2.964369e-02 NA
## Distance.channel 0.000000e+00 0.000000e+00 1.328682e-04
## Vessels.activities.2019 7.912606e-01 1.108240e-10 1.961956e-05
## Wind 0.000000e+00 1.899518e-03 6.618099e-03
## Distance.channel Vessels.activities.2019 Wind
## House.400 0.0000000000 9.413183e-06 0.612090962
## Urban.200 0.0007147639 0.000000e+00 0.059063589
## Order.trial 0.0000000000 4.669252e-01 0.305501234
## JulianDay 0.0000000000 1.393752e-10 0.000000000
## Hour.Platform 0.0000000000 7.912606e-01 0.000000000
## ID.Temperature 0.0000000000 1.108240e-10 0.001899518
## CL 0.0001328682 1.961956e-05 0.006618099
## Distance.channel NA 0.000000e+00 0.000000000
## Vessels.activities.2019 0.0000000000 NA 0.005680785
## Wind 0.0000000000 5.680785e-03 NANo correlation over 0.8
Verification of the assumptions with the initial model
mod.escape.full <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|Site) + (1|ID), data = df, REML = TRUE, na.action=na.exclude)
summary(mod.escape.full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled +
## ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | Site) +
## (1 | ID)
## Data: df
##
## REML criterion at convergence: 2599.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3139 -0.5407 -0.0599 0.4794 3.2714
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.48062 0.6933
## Site (Intercept) 0.02718 0.1648
## Residual 0.72167 0.8495
## Number of obs: 864, groups: ID, 579; Site, 19
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.334345 0.106375 69.268107 31.345
## House.400.scaled 0.027680 0.083392 12.715007 0.332
## Urban.200.scaled 0.109313 0.080455 19.686313 1.359
## Distance.channel.scaled 0.015454 0.066316 17.121120 0.233
## Vessels.activities.2019.scaled -0.034816 0.075154 12.110833 -0.463
## Order.trial.scaled -0.256260 0.042111 627.012641 -6.085
## SexM 0.004507 0.093523 522.440378 0.048
## Tank.SiteT 0.378908 0.161077 412.099292 2.352
## Year2020 -0.061584 0.102827 101.467957 -0.599
## JulianDay.scaled 0.170890 0.093519 21.993588 1.827
## Hour.Platform.scaled -0.091577 0.049252 401.699831 -1.859
## ID.Temperature.scaled -0.117107 0.062842 139.656326 -1.864
## CL.scaled 0.167121 0.050554 501.046102 3.306
## Wind.scaled -0.114240 0.051968 548.743968 -2.198
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## House.400.scaled 0.74535
## Urban.200.scaled 0.18961
## Distance.channel.scaled 0.81849
## Vessels.activities.2019.scaled 0.65139
## Order.trial.scaled 2.02e-09 ***
## SexM 0.96158
## Tank.SiteT 0.01913 *
## Year2020 0.55057
## JulianDay.scaled 0.08125 .
## Hour.Platform.scaled 0.06371 .
## ID.Temperature.scaled 0.06449 .
## CL.scaled 0.00101 **
## Wind.scaled 0.02835 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itplot(mod.escape.full) 
hist(residuals(mod.escape.full)) 
qqnorm(resid(mod.escape.full))
Emergence of the turtle after escaping
Multicollinearity: Generalized variance inflation factor (GVIF^(1/(2*df)))
GVIF^(1/(2*df)) over two indicate the presence of multicollinearity
mod.emergence.vif <- glm(emergence ~ House.300.scaled + Urban.600.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled, data = df, family = binomial, na.action=na.exclude)
gvif(mod.emergence.vif)## GVIF df GVIF^(1/(2*df))
## House.300.scaled 3.4918 1 1.8686
## Urban.600.scaled 4.5572 1 2.1348
## Distance.channel.scaled 1.5181 1 1.2321
## Vessels.activities.2019.scaled 1.5182 1 1.2322
## Order.trial.scaled 1.2502 1 1.1181
## Sex 1.1649 1 1.0793
## Tank.Site 1.9810 1 1.4075
## Year 1.2676 1 1.1259
## JulianDay.scaled 2.1176 1 1.4552
## Hour.Platform.scaled 1.4752 1 1.2146
## ID.Temperature.scaled 1.6635 1 1.2898
## CL.scaled 1.1858 1 1.0889GVIF^(1/(2*df)) of Urban.600.scaled is over two. We need to delete it and recalculate the factors.
mod.emergence.vif <- glm(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled, data = df, family = binomial, na.action=na.exclude)
gvif(mod.emergence.vif)## GVIF df GVIF^(1/(2*df))
## House.300.scaled 1.3616 1 1.1669
## Distance.channel.scaled 1.5179 1 1.2320
## Vessels.activities.2019.scaled 1.5162 1 1.2313
## Order.trial.scaled 1.2071 1 1.0987
## Sex 1.1408 1 1.0681
## Tank.Site 1.9730 1 1.4046
## Year 1.2305 1 1.1093
## JulianDay.scaled 1.9477 1 1.3956
## Hour.Platform.scaled 1.4804 1 1.2167
## ID.Temperature.scaled 1.6126 1 1.2699
## CL.scaled 1.1705 1 1.0819All the GVIF^(1/(2*df)) are now under two.
Calculation of Pearson and Spearman correlation coefficients
Pearson correlation coefficients
# Visualization of the correlations
cor.pearson.emergence <- df[, c(19,22,6,3,10,12,8,16,17)]
chart.Correlation(cor.pearson.emergence, histogram=TRUE, pch=19)
# Creation of the correlation table
table.corr.pearson.emergence <- rcorr(as.matrix(df[,c(19,22,6,3,10,12,8,16,17)]), type="pearson")
table.cor.pearson.emergence.r <- table.corr.pearson.emergence$r # Pearson correlation coefficients
table.cor.pearson.emergence.p <- table.corr.pearson.emergence$P # P values of the correlations
table.cor.pearson.emergence.r## House.300 Urban.600 Order.trial JulianDay
## House.300 1.00000000 0.8172697 0.09704351 -0.4205959
## Urban.600 0.81726971 1.0000000 0.11849156 -0.5377757
## Order.trial 0.09704351 0.1184916 1.00000000 0.1637192
## JulianDay -0.42059588 -0.5377757 0.16371919 1.0000000
## Hour.Platform -0.12340041 -0.1913218 -0.10538411 0.2253572
## ID.Temperature -0.32028077 -0.3236473 0.04136506 0.4406745
## CL 0.12697087 0.2248484 0.03401162 -0.1720253
## Distance.channel -0.17931867 -0.1776160 0.09912351 0.1010733
## Vessels.activities.2019 -0.17880670 -0.2672546 0.01856539 0.3303980
## Hour.Platform ID.Temperature CL
## House.300 -0.12340041 -0.32028077 0.12697087
## Urban.600 -0.19132180 -0.32364730 0.22484835
## Order.trial -0.10538411 0.04136506 0.03401162
## JulianDay 0.22535717 0.44067447 -0.17202529
## Hour.Platform 1.00000000 0.03449218 -0.11957931
## ID.Temperature 0.03449218 1.00000000 -0.06344290
## CL -0.11957931 -0.06344290 1.00000000
## Distance.channel 0.10732729 0.25307372 -0.01469142
## Vessels.activities.2019 0.09445235 0.31822491 -0.08012667
## Distance.channel Vessels.activities.2019
## House.300 -0.17931867 -0.17880670
## Urban.600 -0.17761598 -0.26725456
## Order.trial 0.09912351 0.01856539
## JulianDay 0.10107332 0.33039803
## Hour.Platform 0.10732729 0.09445235
## ID.Temperature 0.25307372 0.31822491
## CL -0.01469142 -0.08012667
## Distance.channel 1.00000000 0.48906628
## Vessels.activities.2019 0.48906628 1.00000000table.cor.pearson.emergence.p## House.300 Urban.600 Order.trial JulianDay
## House.300 NA 0.000000e+00 1.062247e-03 0.000000e+00
## Urban.600 0.000000e+00 NA 6.289997e-05 0.000000e+00
## Order.trial 1.062247e-03 6.289997e-05 NA 2.905276e-08
## JulianDay 0.000000e+00 0.000000e+00 2.905276e-08 NA
## Hour.Platform 3.967329e-05 1.492557e-10 4.553125e-04 3.641532e-14
## ID.Temperature 0.000000e+00 0.000000e+00 1.709756e-01 0.000000e+00
## CL 1.784482e-05 1.776357e-14 2.522456e-01 5.455805e-09
## Distance.channel 1.172760e-09 1.688853e-09 8.257794e-04 6.493434e-04
## Vessels.activities.2019 1.309164e-09 0.000000e+00 5.320837e-01 0.000000e+00
## Hour.Platform ID.Temperature CL
## House.300 3.967329e-05 0.00000000 1.784482e-05
## Urban.600 1.492557e-10 0.00000000 1.776357e-14
## Order.trial 4.553125e-04 0.17097560 2.522456e-01
## JulianDay 3.641532e-14 0.00000000 5.455805e-09
## Hour.Platform NA 0.26029241 6.857072e-05
## ID.Temperature 2.602924e-01 NA 3.563983e-02
## CL 6.857072e-05 0.03563983 NA
## Distance.channel 3.559317e-04 0.00000000 6.210013e-01
## Vessels.activities.2019 1.687374e-03 0.00000000 6.916929e-03
## Distance.channel Vessels.activities.2019
## House.300 1.172760e-09 1.309164e-09
## Urban.600 1.688853e-09 0.000000e+00
## Order.trial 8.257794e-04 5.320837e-01
## JulianDay 6.493434e-04 0.000000e+00
## Hour.Platform 3.559317e-04 1.687374e-03
## ID.Temperature 0.000000e+00 0.000000e+00
## CL 6.210013e-01 6.916929e-03
## Distance.channel NA 0.000000e+00
## Vessels.activities.2019 0.000000e+00 NASpearman correlation coefficients
# Creation of the correlation table
table.corr.spearman.emergence <- rcorr(as.matrix(df[,c(19,22,6,3,10,12,8,16,17)]), type="spearman")
table.cor.spearman.emergence.r <- table.corr.spearman.emergence$r # Spearman correlation coefficients
table.cor.spearman.emergence.p <- table.corr.spearman.emergence$P # P values of the correlations
table.cor.spearman.emergence.r## House.300 Urban.600 Order.trial JulianDay
## House.300 1.00000000 0.6433532 0.11880047 -0.2160482
## Urban.600 0.64335317 1.0000000 -0.11561227 -0.7348046
## Order.trial 0.11880047 -0.1156123 1.00000000 0.2317717
## JulianDay -0.21604821 -0.7348046 0.23177171 1.0000000
## Hour.Platform -0.11310710 -0.2285147 -0.09074473 0.1850147
## ID.Temperature -0.15297542 -0.3819516 0.01087592 0.3952979
## CL 0.07524865 0.1925009 0.02353021 -0.1687139
## Distance.channel -0.11598368 -0.4994053 0.24200697 0.5566049
## Vessels.activities.2019 0.05023256 -0.1063352 -0.02161513 0.1889507
## Hour.Platform ID.Temperature CL
## House.300 -0.113107098 -0.15297542 0.07524865
## Urban.600 -0.228514667 -0.38195164 0.19250093
## Order.trial -0.090744731 0.01087592 0.02353021
## JulianDay 0.185014667 0.39529790 -0.16871392
## Hour.Platform 1.000000000 0.08971374 -0.11614837
## ID.Temperature 0.089713738 1.00000000 -0.06566684
## CL -0.116148369 -0.06566684 1.00000000
## Distance.channel 0.277712076 0.38823857 -0.11317974
## Vessels.activities.2019 -0.007978281 0.19316565 -0.12635136
## Distance.channel Vessels.activities.2019
## House.300 -0.1159837 0.050232564
## Urban.600 -0.4994053 -0.106335150
## Order.trial 0.2420070 -0.021615128
## JulianDay 0.5566049 0.188950707
## Hour.Platform 0.2777121 -0.007978281
## ID.Temperature 0.3882386 0.193165645
## CL -0.1131797 -0.126351359
## Distance.channel 1.0000000 0.368913666
## Vessels.activities.2019 0.3689137 1.000000000table.cor.spearman.emergence.p## House.300 Urban.600 Order.trial JulianDay
## House.300 NA 0.000000e+00 6.016428e-05 1.871836e-13
## Urban.600 0.000000e+00 NA 9.471172e-05 0.000000e+00
## Order.trial 6.016428e-05 9.471172e-05 NA 2.664535e-15
## JulianDay 1.871836e-13 0.000000e+00 2.664535e-15 NA
## Hour.Platform 1.670322e-04 1.554312e-14 2.556506e-03 5.972247e-10
## ID.Temperature 3.563436e-07 0.000000e+00 7.189809e-01 0.000000e+00
## CL 1.121550e-02 6.175815e-11 4.283820e-01 1.073361e-08
## Distance.channel 8.988816e-05 0.000000e+00 0.000000e+00 0.000000e+00
## Vessels.activities.2019 9.073633e-02 3.324111e-04 4.669252e-01 1.393752e-10
## Hour.Platform ID.Temperature CL
## House.300 1.670322e-04 3.563436e-07 1.121550e-02
## Urban.600 1.554312e-14 0.000000e+00 6.175815e-11
## Order.trial 2.556506e-03 7.189809e-01 4.283820e-01
## JulianDay 5.972247e-10 0.000000e+00 1.073361e-08
## Hour.Platform NA 3.357286e-03 1.105565e-04
## ID.Temperature 3.357286e-03 NA 2.964369e-02
## CL 1.105565e-04 2.964369e-02 NA
## Distance.channel 0.000000e+00 0.000000e+00 1.328682e-04
## Vessels.activities.2019 7.912606e-01 1.108240e-10 1.961956e-05
## Distance.channel Vessels.activities.2019
## House.300 8.988816e-05 9.073633e-02
## Urban.600 0.000000e+00 3.324111e-04
## Order.trial 0.000000e+00 4.669252e-01
## JulianDay 0.000000e+00 1.393752e-10
## Hour.Platform 0.000000e+00 7.912606e-01
## ID.Temperature 0.000000e+00 1.108240e-10
## CL 1.328682e-04 1.961956e-05
## Distance.channel NA 0.000000e+00
## Vessels.activities.2019 0.000000e+00 NAThe Pearson correlation coefficients over 0.8 between House.200 and Urban.900 confirm the deletion of Urban.900 with the calculation of the GVIF^(1/(2*df))
Verification of the assumptions with the initial model
Testing for the best optimizer
In preliminary tests, we observed that we had convergence problems for the model with emergence of the turtle as the response variable. Thus, here, we are testing different model optimizer to include in our initial model before the verification of the model assumptions.
#model without optimizer
mod.emergence.full <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID) + (1|Site), data = df, family = binomial, na.action=na.exclude)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.127006 (tol = 0.002, component 1)gm_all <- allFit(mod.emergence.full)## bobyqa : [OK]
## Nelder_Mead :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0991471 (tol = 0.002, component 1)## [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00217746 (tol = 0.002, component 1)## [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]t(sapply(gm_all,fixef)) ## (Intercept) House.300.scaled
## bobyqa -1.039658 -0.2392512
## Nelder_Mead -1.042755 -0.2425455
## nlminbwrap -1.039652 -0.2392475
## nmkbw -1.039433 -0.2392864
## optimx.L-BFGS-B -1.039518 -0.2392432
## nloptwrap.NLOPT_LN_NELDERMEAD -1.039492 -0.2392655
## nloptwrap.NLOPT_LN_BOBYQA -1.039796 -0.2392890
## Distance.channel.scaled
## bobyqa -0.2677396
## Nelder_Mead -0.2617173
## nlminbwrap -0.2677339
## nmkbw -0.2677091
## optimx.L-BFGS-B -0.2676984
## nloptwrap.NLOPT_LN_NELDERMEAD -0.2676344
## nloptwrap.NLOPT_LN_BOBYQA -0.2677534
## Vessels.activities.2019.scaled Order.trial.scaled
## bobyqa 0.2919486 0.2971546
## Nelder_Mead 0.2830340 0.2961752
## nlminbwrap 0.2919385 0.2971546
## nmkbw 0.2919320 0.2971694
## optimx.L-BFGS-B 0.2919241 0.2971528
## nloptwrap.NLOPT_LN_NELDERMEAD 0.2918503 0.2971417
## nloptwrap.NLOPT_LN_BOBYQA 0.2919436 0.2971338
## SexM Tank.SiteT Year2020 JulianDay.scaled
## bobyqa -0.06305510 -1.891496 0.2065169 0.5615806
## Nelder_Mead -0.07079884 -1.883243 0.2158798 0.5611798
## nlminbwrap -0.06304787 -1.891486 0.2065190 0.5615721
## nmkbw -0.06323512 -1.891309 0.2064152 0.5616165
## optimx.L-BFGS-B -0.06305881 -1.891512 0.2064315 0.5616041
## nloptwrap.NLOPT_LN_NELDERMEAD -0.06307785 -1.891350 0.2064084 0.5614811
## nloptwrap.NLOPT_LN_BOBYQA -0.06312965 -1.891319 0.2066913 0.5615796
## Hour.Platform.scaled ID.Temperature.scaled
## bobyqa 0.04908472 0.05329818
## Nelder_Mead 0.04477836 0.05480058
## nlminbwrap 0.04908572 0.05331641
## nmkbw 0.04905547 0.05326379
## optimx.L-BFGS-B 0.04912007 0.05326761
## nloptwrap.NLOPT_LN_NELDERMEAD 0.04907895 0.05346595
## nloptwrap.NLOPT_LN_BOBYQA 0.04897886 0.05330129
## CL.scaled
## bobyqa -0.3144895
## Nelder_Mead -0.3198453
## nlminbwrap -0.3144818
## nmkbw -0.3145370
## optimx.L-BFGS-B -0.3144582
## nloptwrap.NLOPT_LN_NELDERMEAD -0.3144930
## nloptwrap.NLOPT_LN_BOBYQA -0.3145129sapply(gm_all,logLik) ## bobyqa Nelder_Mead
## -542.5624 -542.5706
## nlminbwrap nmkbw
## -542.5624 -542.5624
## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
## -542.5624 -542.5624
## nloptwrap.NLOPT_LN_BOBYQA
## -542.5624We decided to use optimx.L-BFGS-B optimizer to include in our initial model.
mod.emergence.full <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID) + (1|Site), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
summary(mod.emergence.full)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled +
## Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled +
## (1 | ID) + (1 | Site)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1113.1 1181.9 -542.6 1085.1 990
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9600 -0.4582 -0.3078 0.6286 3.9749
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.85404 1.3616
## Site (Intercept) 0.05486 0.2342
## Number of obs: 1004, groups: ID, 684; Site, 22
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03952 0.26366 -3.943 8.06e-05 ***
## House.300.scaled -0.23924 0.14124 -1.694 0.09029 .
## Distance.channel.scaled -0.26770 0.14832 -1.805 0.07109 .
## Vessels.activities.2019.scaled 0.29192 0.13949 2.093 0.03637 *
## Order.trial.scaled 0.29715 0.11282 2.634 0.00844 **
## SexM -0.06306 0.23593 -0.267 0.78925
## Tank.SiteT -1.89151 0.41471 -4.561 5.09e-06 ***
## Year2020 0.20643 0.26491 0.779 0.43583
## JulianDay.scaled 0.56160 0.17940 3.131 0.00175 **
## Hour.Platform.scaled 0.04912 0.12524 0.392 0.69489
## ID.Temperature.scaled 0.05327 0.13938 0.382 0.70234
## CL.scaled -0.31446 0.14102 -2.230 0.02576 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020 JlnDy.
## Hs.300.scld 0.046
## Dstnc.chnn. 0.191 0.169
## Vssl..2019. -0.229 -0.144 -0.449
## Ordr.trl.sc 0.049 -0.146 -0.069 0.115
## SexM -0.540 0.075 -0.016 -0.010 -0.069
## Tank.SiteT 0.054 0.047 -0.135 -0.066 -0.095 0.007
## Year2020 -0.564 -0.067 -0.147 0.227 -0.060 0.014 -0.195
## JulnDy.scld -0.077 0.252 0.232 -0.222 -0.102 0.037 -0.484 0.085
## Hr.Pltfrm.s 0.111 -0.035 0.044 -0.012 0.243 -0.011 -0.515 0.051 0.103
## ID.Tmprtr.s -0.015 0.139 -0.118 -0.079 0.049 0.017 0.306 -0.214 -0.401
## CL.scaled 0.040 -0.003 -0.028 -0.106 -0.161 0.344 0.179 -0.166 -0.042
## Hr.Pl. ID.Tm.
## Hs.300.scld
## Dstnc.chnn.
## Vssl..2019.
## Ordr.trl.sc
## SexM
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.Pltfrm.s
## ID.Tmprtr.s -0.060
## CL.scaled -0.020 -0.057plot(mod.emergence.full) 
hist(residuals(mod.emergence.full)) 
qqnorm(resid(mod.emergence.full))
Model selection
Sum of active defensive behaviors
Random variable
First, we are testing if the identity of the turtle and the sampling site need to be added in the initial model. We are using likelihood ratio tests to test if the addition of these random variables make a significant effect on the initial model. We are using a dummy variable (same value for all the observations) to create a null mixed model to compared with the different combinations of mixed models.
REML is set at FALSE for the selection of random and fixed effects
see this link for more informations: https://stats.stackexchange.com/questions/272633/how-to-decide-whether-to-set-reml-to-true-or-false
Creation of the different mixed models
## Null mixed model
mod.ad.null <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|dummy), data = df, na.action=na.exclude, control=lmerControl(check.nlev.gtr.1="ignore"), REML = FALSE)
## only site identity as random variable
mod.ad.dummy.site <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|dummy) + (1|Site), data = df, control=lmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude, REML = FALSE)
## only turtle identity as random variable
mod.ad.dummy.ID <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|dummy) + (1|ID), data = df, control=lmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude, REML = FALSE)
## Turtle identity without the dummy variable
mod.ad.ID <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
## Turtle and site identity without the dummy variable
mod.ad.ID.site <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)Likelihood ratio tests between the mixed models
The less complex model is always first.
anova(mod.ad.null, mod.ad.dummy.site)## Data: df
## Models:
## mod.ad.null: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy)
## mod.ad.dummy.site: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy) + (1 | Site)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.null 14 2312.9 2380.0 -1142.5 2284.9
## mod.ad.dummy.site 15 2308.7 2380.5 -1139.3 2278.7 6.2863 1 0.01217 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod.ad.null, mod.ad.dummy.ID)## Data: df
## Models:
## mod.ad.null: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy)
## mod.ad.dummy.ID: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.null 14 2312.9 2380.0 -1142.5 2284.9
## mod.ad.dummy.ID 15 2251.5 2323.4 -1110.8 2221.5 63.423 1 1.667e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod.ad.ID, mod.ad.ID.site)## Data: df
## Models:
## mod.ad.ID: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## mod.ad.ID.site: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.ID 14 2249.5 2316.6 -1110.8 2221.5
## mod.ad.ID.site 15 2244.8 2316.6 -1107.4 2214.8 6.7234 1 0.009515 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The test is considered significant when chisq is over 3.84
For this model, we will need to keep turtle and site identity as random variables.
Predictor variables
We are selecting the final model with the backward selection procedure. At each step, we deleted the variable with the higher p value. The deletion of each predictor variable is confirmed with a likelihood ratio test. We need to create a new dataset at each step to have only the rows with complete observations for all the preditors. It allow us to perform the likelihood ratio test that does not run between two models with not the same number of observations.
mod.ad.full <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
summary(mod.ad.full) ## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled +
## ID.Temperature.scaled + CL.scaled + (1 | Site) + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2244.8 2316.6 -1107.4 2214.8 874
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.50194 -0.54956 -0.06127 0.49751 2.86035
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.44568 0.6676
## Site (Intercept) 0.02425 0.1557
## Residual 0.32234 0.5678
## Number of obs: 889, groups: ID, 701; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.824885 0.077979 99.511870 10.578
## House.200.scaled -0.131909 0.052073 24.150494 -2.533
## Distance.channel.scaled -0.148630 0.050410 25.171019 -2.948
## Vessels.activities.2019.scaled 0.196832 0.057428 18.770258 3.427
## Order.trial.scaled 0.021650 0.032774 364.841612 0.661
## SexM 0.298798 0.073782 678.403121 4.050
## Tank.SiteT -0.180610 0.122661 450.731154 -1.472
## Year2020 0.064083 0.077865 360.271219 0.823
## JulianDay.scaled 0.007927 0.062518 59.663331 0.127
## Hour.active.defenses.scaled 0.015288 0.036601 564.780898 0.418
## ID.Temperature.scaled 0.005124 0.044321 298.022016 0.116
## CL.scaled 0.114473 0.039361 639.693627 2.908
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## House.200.scaled 0.01821 *
## Distance.channel.scaled 0.00680 **
## Vessels.activities.2019.scaled 0.00286 **
## Order.trial.scaled 0.50929
## SexM 5.72e-05 ***
## Tank.SiteT 0.14160
## Year2020 0.41105
## JulianDay.scaled 0.89953
## Hour.active.defenses.scaled 0.67632
## ID.Temperature.scaled 0.90804
## CL.scaled 0.00376 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020 JlnDy.
## Hs.200.scld 0.099
## Dstnc.chnn. 0.171 0.124
## Vssl..2019. -0.081 -0.080 -0.484
## Ordr.trl.sc 0.188 -0.046 -0.031 0.082
## SexM -0.613 0.045 -0.018 0.008 -0.039
## Tank.SiteT -0.153 0.001 -0.168 0.050 0.069 -0.041
## Year2020 -0.412 -0.220 -0.177 0.124 -0.115 0.043 -0.244
## JulnDy.scld 0.127 0.244 0.271 -0.295 -0.103 0.073 -0.435 0.052
## Hr.ctv.dfn. 0.095 0.065 0.082 -0.081 0.140 0.015 -0.537 0.177 0.133
## ID.Tmprtr.s 0.014 0.097 -0.119 -0.066 0.005 0.008 0.381 -0.299 -0.463
## CL.scaled -0.278 -0.037 -0.036 0.011 -0.043 0.377 -0.062 0.071 0.105
## Hr.c.. ID.Tm.
## Hs.200.scld
## Dstnc.chnn.
## Vssl..2019.
## Ordr.trl.sc
## SexM
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.ctv.dfn.
## ID.Tmprtr.s -0.036
## CL.scaled 0.063 -0.047We will delete the turtle temperature (ID.Temperature.scaled)
mod.ad.1 <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$active.defenses, df$House.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Sex, df$Tank.Site, df$Year, df$JulianDay.scaled, df$ID.Temperature.scaled ,df$Hour.active.defenses.scaled, df$CL.scaled),]
mod.full.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.ad.1.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.ad.1.adjust, mod.full.adjust)## Data: df.adjust
## Models:
## mod.ad.1.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + CL.scaled + (1 | Site) + (1 | ID)
## mod.full.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + ID.Temperature.scaled + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.1.adjust 14 2242.8 2309.9 -1107.4 2214.8
## mod.full.adjust 15 2244.8 2316.6 -1107.4 2214.8 0.0133 1 0.9081summary(mod.ad.1)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled +
## CL.scaled + (1 | Site) + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2335.5 2403.1 -1153.7 2307.5 912
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.33687 -0.57126 -0.07135 0.50007 2.75759
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4007 0.6330
## Site (Intercept) 0.0248 0.1575
## Residual 0.3592 0.5994
## Number of obs: 926, groups: ID, 708; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.823364 0.077573 100.658603 10.614
## House.200.scaled -0.132758 0.051597 24.410584 -2.573
## Distance.channel.scaled -0.149222 0.050024 25.295331 -2.983
## Vessels.activities.2019.scaled 0.196718 0.057258 18.998104 3.436
## Order.trial.scaled 0.022680 0.032252 404.619126 0.703
## SexM 0.306544 0.072637 685.451716 4.220
## Tank.SiteT -0.172916 0.112791 323.450859 -1.533
## Year2020 0.058627 0.074170 348.137361 0.790
## JulianDay.scaled 0.007504 0.054813 42.013756 0.137
## Hour.active.defenses.scaled 0.018249 0.035376 544.948570 0.516
## CL.scaled 0.112766 0.038546 640.985682 2.925
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## House.200.scaled 0.01657 *
## Distance.channel.scaled 0.00624 **
## Vessels.activities.2019.scaled 0.00277 **
## Order.trial.scaled 0.48232
## SexM 2.77e-05 ***
## Tank.SiteT 0.12624
## Year2020 0.42981
## JulianDay.scaled 0.89176
## Hour.active.defenses.scaled 0.60617
## CL.scaled 0.00356 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020 JlnDy.
## Hs.200.scld 0.096
## Dstnc.chnn. 0.176 0.139
## Vssl..2019. -0.078 -0.072 -0.496
## Ordr.trl.sc 0.191 -0.058 -0.019 0.073
## SexM -0.609 0.046 -0.018 0.008 -0.041
## Tank.SiteT -0.181 -0.038 -0.149 0.093 0.110 -0.048
## Year2020 -0.428 -0.208 -0.223 0.106 -0.124 0.051 -0.134
## JulnDy.scld 0.159 0.326 0.245 -0.365 -0.105 0.088 -0.341 -0.108
## Hr.ctv.dfn. 0.097 0.067 0.093 -0.095 0.103 0.011 -0.552 0.176 0.157
## CL.scaled -0.272 -0.033 -0.040 0.006 -0.050 0.378 -0.048 0.055 0.094
## Hr.c..
## Hs.200.scld
## Dstnc.chnn.
## Vssl..2019.
## Ordr.trl.sc
## SexM
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.ctv.dfn.
## CL.scaled 0.065We will delete the Julian Day (JulianDay.scaled)
mod.ad.2 <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + Hour.active.defenses.scaled + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$active.defenses, df$House.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Sex, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.active.defenses.scaled, df$CL.scaled),]
mod.full.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.ad.2.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + Hour.active.defenses.scaled + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.ad.2.adjust, mod.full.adjust)## Data: df.adjust
## Models:
## mod.ad.2.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + Hour.active.defenses.scaled + CL.scaled + (1 | Site) + (1 | ID)
## mod.full.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.active.defenses.scaled + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.2.adjust 13 2333.5 2396.3 -1153.7 2307.5
## mod.full.adjust 14 2335.5 2403.1 -1153.7 2307.5 0.017 1 0.8961summary(mod.ad.2)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + Hour.active.defenses.scaled + CL.scaled +
## (1 | Site) + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2333.5 2396.3 -1153.7 2307.5 913
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.33642 -0.57110 -0.07145 0.49996 2.75813
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.40048 0.6328
## Site (Intercept) 0.02543 0.1595
## Residual 0.35924 0.5994
## Number of obs: 926, groups: ID, 708; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.82199 0.07679 103.85805 10.705 < 2e-16
## House.200.scaled -0.13492 0.04905 25.28375 -2.751 0.01084
## Distance.channel.scaled -0.15084 0.04878 25.79642 -3.093 0.00472
## Vessels.activities.2019.scaled 0.19964 0.05366 19.20684 3.720 0.00143
## Order.trial.scaled 0.02305 0.03209 402.94796 0.718 0.47301
## SexM 0.30576 0.07236 681.21852 4.226 2.71e-05
## Tank.SiteT -0.16882 0.10623 230.63230 -1.589 0.11337
## Year2020 0.05897 0.07383 332.72104 0.799 0.42503
## Hour.active.defenses.scaled 0.01742 0.03496 612.02658 0.498 0.61843
## CL.scaled 0.11244 0.03838 634.15412 2.929 0.00352
##
## (Intercept) ***
## House.200.scaled *
## Distance.channel.scaled **
## Vessels.activities.2019.scaled **
## Order.trial.scaled
## SexM ***
## Tank.SiteT
## Year2020
## Hour.active.defenses.scaled
## CL.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020 Hr.c..
## Hs.200.scld 0.046
## Dstnc.chnn. 0.142 0.064
## Vssl..2019. -0.021 0.054 -0.450
## Ordr.trl.sc 0.211 -0.025 0.007 0.038
## SexM -0.632 0.018 -0.041 0.044 -0.033
## Tank.SiteT -0.136 0.083 -0.071 -0.035 0.081 -0.019
## Year2020 -0.418 -0.183 -0.203 0.071 -0.137 0.061 -0.183
## Hr.ctv.dfn. 0.073 0.016 0.056 -0.040 0.121 -0.003 -0.537 0.197
## CL.scaled -0.291 -0.067 -0.065 0.043 -0.040 0.373 -0.017 0.066 0.051We will delete the hour of the test (Hour.active.defenses.scaled)
mod.ad.3 <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$active.defenses, df$House.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Sex, df$Tank.Site, df$Year, df$Hour.active.defenses.scaled, df$CL.scaled),]
mod.full.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + Hour.active.defenses.scaled + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.ad.3.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.ad.3.adjust, mod.full.adjust)## Data: df.adjust
## Models:
## mod.ad.3.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + CL.scaled + (1 | Site) + (1 | ID)
## mod.full.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + Hour.active.defenses.scaled + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.3.adjust 12 2331.7 2389.7 -1153.9 2307.7
## mod.full.adjust 13 2333.5 2396.3 -1153.7 2307.5 0.2475 1 0.6189summary(mod.ad.3)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + CL.scaled + (1 | Site) + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2736.4 2796.3 -1356.2 2712.4 1079
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.28155 -0.62349 -0.07709 0.53823 3.10931
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.32853 0.5732
## Site (Intercept) 0.02815 0.1678
## Residual 0.42794 0.6542
## Number of obs: 1091, groups: ID, 714; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.83879 0.07397 92.98861 11.340 < 2e-16
## House.200.scaled -0.14391 0.04914 25.29328 -2.929 0.00711
## Distance.channel.scaled -0.13090 0.04905 25.64563 -2.668 0.01303
## Vessels.activities.2019.scaled 0.19467 0.05338 19.02046 3.647 0.00171
## Order.trial.scaled 0.04384 0.02735 700.47453 1.603 0.10937
## SexM 0.32307 0.06924 653.75623 4.666 3.73e-06
## Tank.SiteT -0.15989 0.07162 569.22734 -2.232 0.02598
## Year2020 0.05069 0.07093 354.53719 0.715 0.47536
## CL.scaled 0.10840 0.03687 635.31677 2.940 0.00340
##
## (Intercept) ***
## House.200.scaled **
## Distance.channel.scaled *
## Vessels.activities.2019.scaled **
## Order.trial.scaled
## SexM ***
## Tank.SiteT *
## Year2020
## CL.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020
## Hs.200.scld 0.054
## Dstnc.chnn. 0.128 0.064
## Vssl..2019. -0.029 0.066 -0.448
## Ordr.trl.sc 0.143 0.012 -0.005 -0.015
## SexM -0.625 0.021 -0.037 0.037 -0.042
## Tank.SiteT -0.047 0.074 -0.031 -0.038 0.355 -0.027
## Year2020 -0.451 -0.175 -0.201 0.060 -0.215 0.053 -0.109
## CL.scaled -0.284 -0.076 -0.066 0.043 -0.034 0.367 0.004 0.061We will delete the year (Year)
mod.ad.4 <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$active.defenses, df$House.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Sex, df$Tank.Site, df$Year, df$CL.scaled),]
mod.full.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.ad.4.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.ad.4.adjust, mod.full.adjust)## Data: df.adjust
## Models:
## mod.ad.4.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + CL.scaled + (1 | Site) + (1 | ID)
## mod.full.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.4.adjust 11 2734.9 2789.8 -1356.4 2712.9
## mod.full.adjust 12 2736.4 2796.3 -1356.2 2712.4 0.483 1 0.4871summary(mod.ad.4)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + CL.scaled + (1 | Site) + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2734.9 2789.8 -1356.4 2712.9 1080
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.28249 -0.61878 -0.07975 0.52423 3.06436
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.32964 0.5741
## Site (Intercept) 0.03064 0.1750
## Residual 0.42689 0.6534
## Number of obs: 1091, groups: ID, 714; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.86202 0.06694 67.25674 12.877 < 2e-16
## House.200.scaled -0.13743 0.04949 25.31289 -2.777 0.01019
## Distance.channel.scaled -0.12406 0.04916 25.13232 -2.524 0.01831
## Vessels.activities.2019.scaled 0.19254 0.05467 19.94307 3.522 0.00215
## Order.trial.scaled 0.04685 0.02675 697.30770 1.751 0.08037
## SexM 0.32092 0.06923 655.95931 4.636 4.3e-06
## Tank.SiteT -0.15784 0.07138 576.98001 -2.211 0.02742
## CL.scaled 0.10740 0.03687 641.22721 2.913 0.00370
##
## (Intercept) ***
## House.200.scaled *
## Distance.channel.scaled *
## Vessels.activities.2019.scaled **
## Order.trial.scaled .
## SexM ***
## Tank.SiteT *
## CL.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 Ordr.. SexM Tnk.ST
## Hs.200.scld -0.028
## Dstnc.chnn. 0.042 0.030
## Vssl..2019. -0.001 0.079 -0.447
## Ordr.trl.sc 0.052 -0.026 -0.048 -0.003
## SexM -0.665 0.030 -0.027 0.033 -0.032
## Tank.SiteT -0.106 0.055 -0.052 -0.032 0.346 -0.022
## CL.scaled -0.285 -0.065 -0.054 0.039 -0.023 0.365 0.010We will delete the order of the trial (Order.trial.scaled)
mod.ad.5 <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Sex + Tank.Site + CL.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$active.defenses, df$House.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Sex, df$Tank.Site, df$CL.scaled),]
mod.full.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.ad.5.adjust <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Sex + Tank.Site + CL.scaled + (1|Site) + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.ad.5.adjust, mod.full.adjust)## Data: df.adjust
## Models:
## mod.ad.5.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Sex + Tank.Site + CL.scaled + (1 | Site) + (1 | ID)
## mod.full.adjust: active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.ad.5.adjust 10 2735.7 2785.7 -1357.9 2715.7
## mod.full.adjust 11 2734.9 2789.8 -1356.4 2712.9 2.8421 1 0.09182 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod.ad.5)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Sex + Tank.Site + CL.scaled +
## (1 | Site) + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2735.7 2785.7 -1357.9 2715.7 1081
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.32255 -0.61721 -0.07972 0.53876 2.97925
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.33610 0.5797
## Site (Intercept) 0.03778 0.1944
## Residual 0.42286 0.6503
## Number of obs: 1091, groups: ID, 714; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.85488 0.06949 63.44592 12.303 < 2e-16
## House.200.scaled -0.13463 0.05259 25.70409 -2.560 0.01671
## Distance.channel.scaled -0.12077 0.05221 25.21780 -2.313 0.02914
## Vessels.activities.2019.scaled 0.19329 0.05851 20.70255 3.303 0.00343
## SexM 0.32590 0.06954 664.53628 4.686 3.37e-06
## Tank.SiteT -0.20887 0.06708 558.21213 -3.114 0.00194
## CL.scaled 0.11033 0.03708 657.68228 2.975 0.00303
##
## (Intercept) ***
## House.200.scaled *
## Distance.channel.scaled *
## Vessels.activities.2019.scaled **
## SexM ***
## Tank.SiteT **
## CL.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 SexM Tnk.ST
## Hs.200.scld -0.028
## Dstnc.chnn. 0.041 0.030
## Vssl..2019. 0.004 0.081 -0.448
## SexM -0.643 0.026 -0.027 0.031
## Tank.SiteT -0.126 0.065 -0.034 -0.031 -0.013
## CL.scaled -0.275 -0.062 -0.053 0.037 0.365 0.016Final model
Summary statistics
We need to fit the REML to TRUE to calculate the summary statistics of the model.
mod.ad.final <- lmer(active.defenses ~ House.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Sex + Tank.Site + CL.scaled + (1|Site) + (1|ID), data = df, REML = TRUE, na.action=na.exclude)
summary(mod.ad.final)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: active.defenses ~ House.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Sex + Tank.Site + CL.scaled +
## (1 | Site) + (1 | ID)
## Data: df
##
## REML criterion at convergence: 2743.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.30641 -0.60610 -0.08527 0.52494 2.97977
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.33808 0.5814
## Site (Intercept) 0.05099 0.2258
## Residual 0.42343 0.6507
## Number of obs: 1091, groups: ID, 714; Site, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.85355 0.07394 45.37461 11.544 4.25e-15
## House.200.scaled -0.13390 0.05781 20.46251 -2.316 0.03102
## Distance.channel.scaled -0.12190 0.05744 19.99195 -2.122 0.04650
## Vessels.activities.2019.scaled 0.19395 0.06495 17.02346 2.986 0.00829
## SexM 0.32747 0.06987 660.18254 4.687 3.38e-06
## Tank.SiteT -0.21805 0.06754 550.58716 -3.228 0.00132
## CL.scaled 0.11243 0.03731 653.49939 3.014 0.00268
##
## (Intercept) ***
## House.200.scaled *
## Distance.channel.scaled *
## Vessels.activities.2019.scaled **
## SexM ***
## Tank.SiteT **
## CL.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.200. Dstn.. V..201 SexM Tnk.ST
## Hs.200.scld -0.030
## Dstnc.chnn. 0.036 0.031
## Vssl..2019. 0.011 0.083 -0.449
## SexM -0.607 0.023 -0.025 0.027
## Tank.SiteT -0.118 0.059 -0.029 -0.029 -0.014
## CL.scaled -0.261 -0.058 -0.050 0.033 0.366 0.012Quick visualization of the preditor effects
plot(allEffects(mod.ad.final)) 
Calculation of the marginal and conditional variance explained by the final model
r.squaredGLMM(mod.ad.final)## R2m R2c
## [1,] 0.08179043 0.5214769marginal: only for the fixed effect. conditional: both fixed and random effect.
Calculation of the 95% confidence intervals
confint(mod.ad.final, level = 0.95, method = "Wald")## 2.5 % 97.5 %
## .sig01 NA NA
## .sig02 NA NA
## .sigma NA NA
## (Intercept) 0.70863757 0.99845973
## House.200.scaled -0.24721384 -0.02059445
## Distance.channel.scaled -0.23448490 -0.00932149
## Vessels.activities.2019.scaled 0.06664645 0.32124405
## SexM 0.19051818 0.46441761
## Tank.SiteT -0.35042096 -0.08567342
## CL.scaled 0.03930962 0.18555949Creation of the prediction figures for variables assessing human disturbance
Creation of the predictions of the final model for each variable
## Number of house with access to canal within 200 m
pred.mod.ad.house200 <- ggpredict(mod.ad.final, terms = "House.200.scaled")
pred.mod.ad.house200## # Predicted values of active.defenses
##
## House.200.scaled | Predicted | 95% CI
## -------------------------------------------
## -0.76 | 0.96 | [0.79, 1.13]
## 0.02 | 0.86 | [0.71, 1.00]
## 0.81 | 0.75 | [0.58, 0.92]
## 1.59 | 0.65 | [0.42, 0.87]
## 2.37 | 0.54 | [0.24, 0.84]
## 3.15 | 0.44 | [0.06, 0.82]
##
## Adjusted for:
## * Distance.channel.scaled = -0.01
## * Vessels.activities.2019.scaled = -0.01
## * Sex = F
## * Tank.Site = S
## * CL.scaled = 0.06
## * Site = 0 (population-level)
## * ID = 0 (population-level)## Shortest aquatic distance to the navigation channel
pred.mod.ad.channel <- ggpredict(mod.ad.final, terms = "Distance.channel.scaled")
pred.mod.ad.channel## # Predicted values of active.defenses
##
## Distance.channel.scaled | Predicted | 95% CI
## ---------------------------------------------------
## -1.00 | 0.98 | [ 0.80, 1.16]
## -0.50 | 0.92 | [ 0.77, 1.07]
## 0.50 | 0.80 | [ 0.64, 0.95]
## 1.50 | 0.68 | [ 0.45, 0.90]
## 2.00 | 0.62 | [ 0.34, 0.89]
## 2.50 | 0.55 | [ 0.23, 0.87]
## 3.50 | 0.43 | [ 0.01, 0.86]
## 5.00 | 0.25 | [-0.34, 0.84]
##
## Adjusted for:
## * House.200.scaled = 0.00
## * Vessels.activities.2019.scaled = -0.01
## * Sex = F
## * Tank.Site = S
## * CL.scaled = 0.06
## * Site = 0 (population-level)
## * ID = 0 (population-level)## Mean daily number of vessel crossings
pred.mod.ad.vessels <- ggpredict(mod.ad.final, terms = "Vessels.activities.2019.scaled")
pred.mod.ad.vessels## # Predicted values of active.defenses
##
## Vessels.activities.2019.scaled | Predicted | 95% CI
## ---------------------------------------------------------
## -0.84 | 0.70 | [0.52, 0.88]
## -0.71 | 0.72 | [0.55, 0.89]
## -0.63 | 0.74 | [0.58, 0.90]
## -0.61 | 0.74 | [0.58, 0.91]
## -0.56 | 0.75 | [0.59, 0.91]
## -0.54 | 0.76 | [0.60, 0.92]
## 0.03 | 0.87 | [0.72, 1.01]
## 2.56 | 1.36 | [1.00, 1.72]
##
## Adjusted for:
## * House.200.scaled = 0.00
## * Distance.channel.scaled = -0.01
## * Sex = F
## * Tank.Site = S
## * CL.scaled = 0.06
## * Site = 0 (population-level)
## * ID = 0 (population-level)## New dataset to only have complete observations for all the variables
df.ajust.final.model <- df[complete.cases(df$active.defenses, df$House.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Sex, df$Tank.Site, df$CL.scaled),]Figure for the number of houses with access to the canal within 200m
graph.house200 <- ggplot(data=pred.mod.ad.house200, aes(x, predicted)) + geom_point(data=df.ajust.final.model, aes(House.200.scaled, active.defenses), position = "jitter") +
geom_ribbon(data=pred.mod.ad.house200, aes(ymin=conf.low, ymax= conf.high), alpha=0.6) +
geom_line(data=pred.mod.ad.house200, color="black")+
theme_bw() + theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
theme(axis.text=element_text(size=15, colour="black",face = "bold")) +
ylab("Sum of active defensive behaviors") +
xlab("Number of houses within 200m")
graph.house200
Figure for the shortest aquatic distance to the navigation channel
graph.channel <- ggplot(data=pred.mod.ad.channel, aes(x, predicted)) + geom_point(data=df.ajust.final.model, aes(Distance.channel.scaled, active.defenses), position = "jitter") +
geom_ribbon(data=pred.mod.ad.channel, aes(ymin=conf.low, ymax= conf.high), alpha=0.6) +
geom_line(data=pred.mod.ad.channel, color="black")+
theme_bw() + theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
theme(axis.text=element_text(size=15, colour="black",face = "bold")) +
ylab("Sum of active defensive behaviors") +
xlab("Shortest distance to the navigation channel")
graph.channel
Figure for the mean daily number of vessel crossings
graph.vessels <- ggplot(data=pred.mod.ad.vessels, aes(x, predicted)) + geom_point(data=df.ajust.final.model, aes(Vessels.activities.2019.scaled, active.defenses), position = "jitter") +
geom_ribbon(data=pred.mod.ad.vessels, aes(ymin=conf.low, ymax= conf.high), alpha=0.6) +
geom_line(data=pred.mod.ad.vessels, color="black")+
theme_bw() + theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
theme(axis.text=element_text(size=15, colour="black",face = "bold")) +
ylab("Sum of active defensive behaviors") +
xlab("Mean daily number of vessel crossings")
graph.vessels
Escape latency
Random variable
First, we are testing if the identity of the turtle and the sampling site need to be added in the initial model. We are using likelihood ratio tests to test if the addition of these random variables make a significant effect on the initial model. We are using a dummy variable (same value for all the observations) to create a null mixed model to compared with the different combinations of mixed models.
REML is set at FALSE for the selection of random and fixed effects
see this link for more informations: https://stats.stackexchange.com/questions/272633/how-to-decide-whether-to-set-reml-to-true-or-false
Creation of the different mixed models
## Null mixed model
mod.escape.null <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|dummy), data = df, na.action=na.exclude, control=lmerControl(check.nlev.gtr.1="ignore"), REML = FALSE)
## only site identity as random variable
mod.escape.dummy.site <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|dummy) + (1|Site), data = df, control=lmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude, REML = FALSE)
## only turtle identity as random variable
mod.escape.dummy.ID <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|dummy) + (1|ID), data = df, control=lmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude, REML = FALSE)
## Turtle identity without the dummy variable
mod.escape.ID <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
## Turtle and site identity without the dummy variable
mod.escape.ID.site <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|Site) + (1|ID), data = df, na.action=na.exclude, REML = FALSE)Likelihood ratio tests between the mixed models
The less complex model is always first.
anova(mod.escape.null, mod.escape.dummy.site)## Data: df
## Models:
## mod.escape.null: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | dummy)
## mod.escape.dummy.site: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | dummy) + (1 | Site)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.null 16 2630.6 2706.8 -1299.3 2598.6
## mod.escape.dummy.site 17 2630.8 2711.7 -1298.4 2596.8 1.8206 1 0.1772anova(mod.escape.null, mod.escape.dummy.ID)## Data: df
## Models:
## mod.escape.null: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | dummy)
## mod.escape.dummy.ID: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | dummy) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.null 16 2630.6 2706.8 -1299.3 2598.6
## mod.escape.dummy.ID 17 2580.5 2661.4 -1273.2 2546.5 52.087 1 5.309e-13
##
## mod.escape.null
## mod.escape.dummy.ID ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod.escape.ID, mod.escape.ID.site)## Data: df
## Models:
## mod.escape.ID: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.escape.ID.site: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.ID 16 2578.5 2654.7 -1273.2 2546.5
## mod.escape.ID.site 17 2580.2 2661.2 -1273.1 2546.2 0.2413 1 0.6233Site identity does not have a significant effect on the initial model. We are only keeping turtle identity in our model.
Predictor variables
We are selecting the final model with the backward selection procedure. At each step, we deleted the variable with the higher p value. The deletion of each predictor variable is confirmed with a likelihood ratio test. We need to create a new dataset at each step to have only the rows with complete observations for all the preditors. It allow us to perform the likelihood ratio test that does not run between two models with not the same number of observations.
mod.escape.full <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
summary(mod.escape.full)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Sex +
## Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled +
## ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2578.5 2654.7 -1273.2 2546.5 848
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3594 -0.5406 -0.0574 0.4792 3.2932
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4736 0.6882
## Residual 0.7210 0.8491
## Number of obs: 864, groups: ID, 579
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.359760 0.094461 682.902375 35.568
## House.400.scaled 0.026837 0.061716 570.910245 0.435
## Urban.200.scaled 0.093267 0.063690 472.867253 1.464
## Distance.channel.scaled 0.016669 0.051795 602.037986 0.322
## Vessels.activities.2019.scaled -0.040090 0.054866 508.910133 -0.731
## Order.trial.scaled -0.255849 0.041091 736.238739 -6.226
## SexM 0.004524 0.091650 547.696647 0.049
## Tank.SiteT 0.341335 0.152755 854.698187 2.235
## Year2020 -0.099757 0.091444 695.706148 -1.091
## JulianDay.scaled 0.157908 0.073646 613.763768 2.144
## Hour.Platform.scaled -0.092431 0.046730 861.998943 -1.978
## ID.Temperature.scaled -0.141341 0.056543 826.634692 -2.500
## CL.scaled 0.167077 0.049090 564.224148 3.403
## Wind.scaled -0.130198 0.049839 850.478505 -2.612
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## House.400.scaled 0.663834
## Urban.200.scaled 0.143752
## Distance.channel.scaled 0.747695
## Vessels.activities.2019.scaled 0.465302
## Order.trial.scaled 8.02e-10 ***
## SexM 0.960647
## Tank.SiteT 0.025706 *
## Year2020 0.275692
## JulianDay.scaled 0.032413 *
## Hour.Platform.scaled 0.048247 *
## ID.Temperature.scaled 0.012622 *
## CL.scaled 0.000713 ***
## Wind.scaled 0.009150 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itWe will delete sex (Sex)
mod.escape.1 <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$House.400.scaled, df$Urban.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Sex, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.1.adjust <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.1.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.1.adjust: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.1.adjust 15 2576.5 2647.9 -1273.2 2546.5
## mod.full.adjust 16 2578.5 2654.7 -1273.2 2546.5 0.0024 1 0.9606summary(mod.escape.1)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled +
## Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site +
## Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled +
## CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2656.1 2727.9 -1313.0 2626.1 874
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3078 -0.5526 -0.0605 0.4765 3.2749
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4451 0.6672
## Residual 0.7483 0.8650
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.37373 0.07396 771.99286 45.618 < 2e-16
## House.400.scaled 0.04578 0.06002 594.97690 0.763 0.44590
## Urban.200.scaled 0.08225 0.06190 476.89894 1.329 0.18462
## Distance.channel.scaled 0.02503 0.05131 617.01907 0.488 0.62584
## Vessels.activities.2019.scaled -0.05112 0.05357 517.75618 -0.954 0.34041
## Order.trial.scaled -0.26135 0.04092 776.13527 -6.386 2.92e-10
## Tank.SiteT 0.29331 0.15160 882.12001 1.935 0.05334
## Year2020 -0.08493 0.08971 713.47763 -0.947 0.34408
## JulianDay.scaled 0.16223 0.07241 619.91443 2.240 0.02541
## Hour.Platform.scaled -0.08526 0.04660 887.94279 -1.830 0.06766
## ID.Temperature.scaled -0.15778 0.05560 844.78424 -2.838 0.00465
## CL.scaled 0.15343 0.04196 568.75964 3.656 0.00028
## Wind.scaled -0.14789 0.04950 878.08835 -2.988 0.00289
##
## (Intercept) ***
## House.400.scaled
## Urban.200.scaled
## Distance.channel.scaled
## Vessels.activities.2019.scaled
## Order.trial.scaled ***
## Tank.SiteT .
## Year2020
## JulianDay.scaled *
## Hour.Platform.scaled .
## ID.Temperature.scaled **
## CL.scaled ***
## Wind.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itWe will delete the distance to the navigation channel (Distance.channel.scaled)
mod.escape.2 <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$House.400.scaled, df$Urban.200.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.2.adjust <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.2.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.2.adjust: log.escape ~ House.400.scaled + Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ House.400.scaled + Urban.200.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.2.adjust 14 2654.3 2721.4 -1313.2 2626.3
## mod.full.adjust 15 2656.1 2727.9 -1313.0 2626.1 0.2368 1 0.6266summary(mod.escape.2)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## log.escape ~ House.400.scaled + Urban.200.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + Year + JulianDay.scaled +
## Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled +
## Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2654.3 2721.4 -1313.2 2626.3 875
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3164 -0.5501 -0.0611 0.4815 3.2641
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4476 0.6690
## Residual 0.7469 0.8642
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.36989 0.07357 774.87274 45.804 < 2e-16
## House.400.scaled 0.03198 0.05292 632.25994 0.604 0.545860
## Urban.200.scaled 0.08655 0.06130 486.87999 1.412 0.158634
## Vessels.activities.2019.scaled -0.03515 0.04237 558.71779 -0.830 0.407154
## Order.trial.scaled -0.25938 0.04072 771.49261 -6.370 3.24e-10
## Tank.SiteT 0.30520 0.14960 884.67570 2.040 0.041640
## Year2020 -0.08125 0.08946 722.69354 -0.908 0.364049
## JulianDay.scaled 0.14790 0.06628 604.34004 2.231 0.026025
## Hour.Platform.scaled -0.08428 0.04657 888.11020 -1.810 0.070688
## ID.Temperature.scaled -0.15366 0.05499 846.42851 -2.794 0.005318
## CL.scaled 0.15478 0.04190 568.53209 3.694 0.000242
## Wind.scaled -0.14656 0.04943 878.82156 -2.965 0.003107
##
## (Intercept) ***
## House.400.scaled
## Urban.200.scaled
## Vessels.activities.2019.scaled
## Order.trial.scaled ***
## Tank.SiteT *
## Year2020
## JulianDay.scaled *
## Hour.Platform.scaled .
## ID.Temperature.scaled **
## CL.scaled ***
## Wind.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.400. U.200. V..201 Ordr.. Tnk.ST Yr2020 JlnDy. Hr.Pl.
## Hs.400.scld -0.238
## Urbn.200.sc 0.013 -0.198
## Vssl..2019. 0.026 -0.287 0.091
## Ordr.trl.sc 0.264 -0.045 -0.231 0.005
## Tank.SiteT -0.385 0.217 0.019 -0.074 -0.003
## Year2020 -0.643 0.204 -0.063 0.013 -0.223 -0.128
## JulnDy.scld -0.195 0.232 0.294 -0.201 -0.251 -0.203 0.258
## Hr.Pltfrm.s 0.362 -0.145 -0.004 0.022 0.216 -0.534 -0.092 -0.062
## ID.Tmprtr.s -0.176 0.011 0.174 -0.207 -0.002 0.377 -0.186 -0.147 -0.155
## CL.scaled -0.089 -0.085 -0.031 0.015 -0.035 -0.014 0.086 0.140 0.015
## Wind.scaled -0.168 0.165 0.084 -0.002 0.001 0.558 -0.183 0.013 -0.069
## ID.Tm. CL.scl
## Hs.400.scld
## Urbn.200.sc
## Vssl..2019.
## Ordr.trl.sc
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.Pltfrm.s
## ID.Tmprtr.s
## CL.scaled 0.020
## Wind.scaled 0.143 -0.023We will delete the number of houses within 400m (House.400.scaled)
mod.escape.3 <- lmer(log.escape ~ Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$House.400.scaled, df$Urban.200.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ House.400.scaled + Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.3.adjust <- lmer(log.escape ~ Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.3.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.3.adjust: log.escape ~ Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ House.400.scaled + Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.3.adjust 13 2652.7 2714.9 -1313.3 2626.7
## mod.full.adjust 14 2654.3 2721.4 -1313.2 2626.3 0.365 1 0.5457summary(mod.escape.3)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: log.escape ~ Urban.200.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + Year + JulianDay.scaled +
## Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled +
## Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2652.7 2714.9 -1313.3 2626.7 876
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3200 -0.5538 -0.0604 0.4847 3.2628
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4472 0.6687
## Residual 0.7476 0.8646
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.38048 0.07146 752.34238 47.304 < 2e-16
## Urban.200.scaled 0.09392 0.06009 492.77974 1.563 0.118673
## Vessels.activities.2019.scaled -0.02778 0.04058 582.19727 -0.685 0.493899
## Order.trial.scaled -0.25827 0.04069 773.78365 -6.348 3.72e-10
## Tank.SiteT 0.28559 0.14606 888.80644 1.955 0.050862
## Year2020 -0.09229 0.08759 728.07789 -1.054 0.292396
## JulianDay.scaled 0.13861 0.06447 605.28452 2.150 0.031957
## Hour.Platform.scaled -0.08022 0.04609 888.63947 -1.741 0.082113
## ID.Temperature.scaled -0.15406 0.05500 845.30185 -2.801 0.005206
## CL.scaled 0.15693 0.04175 564.03520 3.759 0.000189
## Wind.scaled -0.15149 0.04876 880.39775 -3.107 0.001952
##
## (Intercept) ***
## Urban.200.scaled
## Vessels.activities.2019.scaled
## Order.trial.scaled ***
## Tank.SiteT .
## Year2020
## JulianDay.scaled *
## Hour.Platform.scaled .
## ID.Temperature.scaled **
## CL.scaled ***
## Wind.scaled **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) U.200. V..201 Ordr.. Tnk.ST Yr2020 JlnDy. Hr.Pl. ID.Tm.
## Urbn.200.sc -0.036
## Vssl..2019. -0.046 0.036
## Ordr.trl.sc 0.261 -0.245 -0.008
## Tank.SiteT -0.351 0.065 -0.012 0.006
## Year2020 -0.625 -0.024 0.076 -0.219 -0.180
## JulnDy.scld -0.148 0.357 -0.144 -0.248 -0.267 0.221
## Hr.Pltfrm.s 0.341 -0.034 -0.021 0.212 -0.521 -0.064 -0.030
## ID.Tmprtr.s -0.178 0.180 -0.213 -0.002 0.383 -0.192 -0.154 -0.155
## CL.scaled -0.113 -0.049 -0.010 -0.039 0.005 0.106 0.165 0.003 0.021
## Wind.scaled -0.134 0.121 0.048 0.008 0.542 -0.224 -0.026 -0.046 0.143
## CL.scl
## Urbn.200.sc
## Vssl..2019.
## Ordr.trl.sc
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.Pltfrm.s
## ID.Tmprtr.s
## CL.scaled
## Wind.scaled -0.009We will delete the mean daily number of vessel crossings (Vessel.activities.2019.scaled)
mod.escape.4 <- lmer(log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$Urban.200.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.4.adjust <- lmer(log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.4.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.4.adjust: log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ Urban.200.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.4.adjust 12 2651.2 2708.6 -1313.6 2627.2
## mod.full.adjust 13 2652.7 2714.9 -1313.3 2626.7 0.4685 1 0.4937summary(mod.escape.4)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site +
## Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled +
## CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2651.1 2708.6 -1313.6 2627.1 877
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3308 -0.5591 -0.0591 0.4905 3.2651
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4470 0.6686
## Residual 0.7482 0.8650
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.37824 0.07140 758.82702 47.312 < 2e-16 ***
## Urban.200.scaled 0.09540 0.06006 491.85232 1.589 0.112804
## Order.trial.scaled -0.25850 0.04070 773.75364 -6.352 3.63e-10 ***
## Tank.SiteT 0.28440 0.14609 888.86539 1.947 0.051875 .
## Year2020 -0.08775 0.08735 738.14739 -1.005 0.315443
## JulianDay.scaled 0.13227 0.06381 609.96917 2.073 0.038610 *
## Hour.Platform.scaled -0.08089 0.04609 888.45303 -1.755 0.079585 .
## ID.Temperature.scaled -0.16208 0.05375 830.54029 -3.016 0.002643 **
## CL.scaled 0.15665 0.04176 563.99610 3.751 0.000194 ***
## Wind.scaled -0.14989 0.04872 880.29271 -3.077 0.002158 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) U.200. Ordr.. Tnk.ST Yr2020 JlnDy. Hr.Pl. ID.Tm. CL.scl
## Urbn.200.sc -0.035
## Ordr.trl.sc 0.261 -0.245
## Tank.SiteT -0.352 0.066 0.006
## Year2020 -0.624 -0.026 -0.219 -0.180
## JulnDy.scld -0.156 0.366 -0.252 -0.272 0.235
## Hr.Pltfrm.s 0.341 -0.033 0.212 -0.521 -0.063 -0.033
## ID.Tmprtr.s -0.193 0.192 -0.004 0.390 -0.181 -0.190 -0.163
## CL.scaled -0.113 -0.049 -0.039 0.005 0.107 0.165 0.003 0.020
## Wind.scaled -0.132 0.119 0.009 0.544 -0.228 -0.019 -0.045 0.157 -0.009We will delete year (Year)
mod.escape.5 <- lmer(log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$Urban.200.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.5.adjust <- lmer(log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.5.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.5.adjust: log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.5.adjust 11 2650.2 2702.8 -1314.1 2628.2
## mod.full.adjust 12 2651.2 2708.6 -1313.6 2627.2 1.0058 1 0.3159summary(mod.escape.5)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site +
## JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled +
## CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2650.2 2702.8 -1314.1 2628.2 878
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3344 -0.5547 -0.0696 0.4788 3.2674
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4511 0.6716
## Residual 0.7466 0.8641
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.33356 0.05587 798.55711 59.668 < 2e-16 ***
## Urban.200.scaled 0.09373 0.06014 491.10802 1.559 0.119751
## Order.trial.scaled -0.26742 0.03971 759.10858 -6.734 3.27e-11 ***
## Tank.SiteT 0.25799 0.14379 888.98576 1.794 0.073126 .
## JulianDay.scaled 0.14718 0.06211 610.70883 2.370 0.018108 *
## Hour.Platform.scaled -0.08368 0.04602 887.83908 -1.818 0.069357 .
## ID.Temperature.scaled -0.17171 0.05291 814.07305 -3.245 0.001222 **
## CL.scaled 0.16114 0.04158 565.06718 3.875 0.000119 ***
## Wind.scaled -0.16100 0.04745 885.05546 -3.393 0.000723 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) U.200. Ordr.. Tnk.ST JlnDy. Hr.Pl. ID.Tm. CL.scl
## Urbn.200.sc -0.065
## Ordr.trl.sc 0.163 -0.257
## Tank.SiteT -0.603 0.062 -0.033
## JulnDy.scld -0.012 0.383 -0.211 -0.240
## Hr.Pltfrm.s 0.386 -0.035 0.203 -0.542 -0.019
## ID.Tmprtr.s -0.397 0.190 -0.045 0.369 -0.155 -0.178
## CL.scaled -0.060 -0.046 -0.016 0.024 0.145 0.010 0.039
## Wind.scaled -0.361 0.116 -0.044 0.525 0.036 -0.061 0.121 0.016We will delete the proportion of urban area within 200m (Urban.200.scaled)
mod.escape.6 <- lmer(log.escape ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$Urban.200.scaled, df$Order.trial.scaled, df$Tank.Site, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.6.adjust <- lmer(log.escape ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.6.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.6.adjust: log.escape ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ Urban.200.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.6.adjust 10 2650.6 2698.5 -1315.3 2630.6
## mod.full.adjust 11 2650.2 2702.8 -1314.1 2628.2 2.4176 1 0.12summary(mod.escape.6)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: log.escape ~ Order.trial.scaled + Tank.Site + JulianDay.scaled +
## Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled +
## Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2650.6 2698.5 -1315.3 2630.6 879
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3637 -0.5388 -0.0643 0.4778 3.2756
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4584 0.6770
## Residual 0.7444 0.8628
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.33919 0.05587 784.81329 59.767 < 2e-16 ***
## Order.trial.scaled -0.25158 0.03840 804.18070 -6.551 1.02e-10 ***
## Tank.SiteT 0.24428 0.14372 888.68953 1.700 0.089540 .
## JulianDay.scaled 0.10996 0.05754 595.28856 1.911 0.056506 .
## Hour.Platform.scaled -0.08101 0.04605 887.79102 -1.759 0.078918 .
## ID.Temperature.scaled -0.18708 0.05206 787.71270 -3.594 0.000346 ***
## CL.scaled 0.16416 0.04166 567.08104 3.940 9.16e-05 ***
## Wind.scaled -0.16939 0.04719 887.52895 -3.589 0.000350 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ordr.. Tnk.ST JlnDy. Hr.Pl. ID.Tm. CL.scl
## Ordr.trl.sc 0.152
## Tank.SiteT -0.601 -0.016
## JulnDy.scld 0.013 -0.127 -0.285
## Hr.Pltfrm.s 0.385 0.201 -0.542 -0.006
## ID.Tmprtr.s -0.393 0.004 0.365 -0.252 -0.174
## CL.scaled -0.063 -0.029 0.027 0.177 0.008 0.049
## Wind.scaled -0.356 -0.014 0.522 -0.009 -0.057 0.101 0.021We will delete the global environment of the test (Tank.Site)
mod.escape.7 <- lmer(log.escape ~ Order.trial.scaled + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$Order.trial.scaled, df$Tank.Site, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.7.adjust <- lmer(log.escape ~ Order.trial.scaled + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.7.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.7.adjust: log.escape ~ Order.trial.scaled + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.7.adjust 9 2651.4 2694.6 -1316.7 2633.4
## mod.full.adjust 10 2650.6 2698.5 -1315.3 2630.6 2.8841 1 0.08946 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod.escape.7)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## log.escape ~ Order.trial.scaled + JulianDay.scaled + Hour.Platform.scaled +
## ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2651.5 2694.6 -1316.7 2633.5 880
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3041 -0.5405 -0.0643 0.4724 3.2829
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4590 0.6775
## Residual 0.7474 0.8645
## Number of obs: 889, groups: ID, 594
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.39623 0.04473 528.39160 75.922 < 2e-16 ***
## Order.trial.scaled -0.25048 0.03847 811.48122 -6.512 1.30e-10 ***
## JulianDay.scaled 0.13786 0.05523 548.86303 2.496 0.012854 *
## Hour.Platform.scaled -0.03864 0.03878 877.89161 -0.996 0.319295
## ID.Temperature.scaled -0.21943 0.04854 807.58948 -4.521 7.09e-06 ***
## CL.scaled 0.16226 0.04171 565.67264 3.890 0.000112 ***
## Wind.scaled -0.21129 0.04032 888.46792 -5.241 2.00e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ordr.. JlnDy. Hr.Pl. ID.Tm. CL.scl
## Ordr.trl.sc 0.178
## JulnDy.scld -0.207 -0.138
## Hr.Pltfrm.s 0.089 0.228 -0.199
## ID.Tmprtr.s -0.233 0.011 -0.165 0.030
## CL.scaled -0.059 -0.029 0.192 0.027 0.042
## Wind.scaled -0.062 -0.006 0.171 0.314 -0.113 0.009We will delete the hour (Hour.Platform)
mod.escape.8 <- lmer(log.escape ~ Order.trial.scaled + JulianDay.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = FALSE)
df.adjust <- df[complete.cases(df$log.escape, df$Order.trial.scaled, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled, df$Wind.scaled),]
mod.full.adjust <- lmer(log.escape ~ Order.trial.scaled + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
mod.escape.8.adjust <- lmer(log.escape ~ Order.trial.scaled + JulianDay.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df.adjust, na.action=na.exclude, REML = FALSE)
anova(mod.escape.8.adjust, mod.full.adjust) ## Data: df.adjust
## Models:
## mod.escape.8.adjust: log.escape ~ Order.trial.scaled + JulianDay.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## mod.full.adjust: log.escape ~ Order.trial.scaled + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.escape.8.adjust 8 2650.4 2688.8 -1317.2 2634.4
## mod.full.adjust 9 2651.4 2694.6 -1316.7 2633.4 0.9897 1 0.3198summary(mod.escape.8)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## log.escape ~ Order.trial.scaled + JulianDay.scaled + ID.Temperature.scaled +
## CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 2732.7 2771.3 -1358.4 2716.7 911
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3095 -0.5513 -0.0602 0.4761 3.2045
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4515 0.6720
## Residual 0.7514 0.8669
## Number of obs: 919, groups: ID, 596
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.40111 0.04426 549.82575 76.849 < 2e-16 ***
## Order.trial.scaled -0.23129 0.03598 855.69664 -6.428 2.14e-10 ***
## JulianDay.scaled 0.14069 0.05293 525.67320 2.658 0.0081 **
## ID.Temperature.scaled -0.19946 0.04676 868.37804 -4.266 2.21e-05 ***
## CL.scaled 0.16694 0.04144 582.86640 4.029 6.35e-05 ***
## Wind.scaled -0.17599 0.03580 893.79900 -4.916 1.05e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ordr.. JlnDy. ID.Tm. CL.scl
## Ordr.trl.sc 0.150
## JulnDy.scld -0.204 -0.151
## ID.Tmprtr.s -0.241 0.038 -0.144
## CL.scaled -0.063 -0.036 0.207 0.033
## Wind.scaled -0.122 -0.195 0.205 -0.109 -0.001No variable related to human disturbance is significant.
Final model
Summary statistics
We need to fit the REML to TRUE to calculate the summary statistics of the model.
mod.escape.final <- lmer(log.escape ~ Order.trial.scaled + JulianDay.scaled + ID.Temperature.scaled + CL.scaled + Wind.scaled + (1|ID), data = df, na.action=na.exclude, REML = TRUE)
summary(mod.escape.final)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log.escape ~ Order.trial.scaled + JulianDay.scaled + ID.Temperature.scaled +
## CL.scaled + Wind.scaled + (1 | ID)
## Data: df
##
## REML criterion at convergence: 2743.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3047 -0.5472 -0.0598 0.4751 3.1972
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4576 0.6764
## Residual 0.7543 0.8685
## Number of obs: 919, groups: ID, 596
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.40101 0.04444 546.75913 76.525 < 2e-16 ***
## Order.trial.scaled -0.23150 0.03609 849.72025 -6.415 2.33e-10 ***
## JulianDay.scaled 0.14062 0.05316 522.77681 2.645 0.00841 **
## ID.Temperature.scaled -0.19920 0.04692 862.74430 -4.245 2.42e-05 ***
## CL.scaled 0.16696 0.04161 579.05456 4.013 6.79e-05 ***
## Wind.scaled -0.17591 0.03591 888.31626 -4.899 1.15e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ordr.. JlnDy. ID.Tm. CL.scl
## Ordr.trl.sc 0.151
## JulnDy.scld -0.204 -0.151
## ID.Tmprtr.s -0.241 0.038 -0.144
## CL.scaled -0.063 -0.036 0.207 0.033
## Wind.scaled -0.123 -0.195 0.205 -0.110 -0.001Quick visualization of the preditor effects
plot(allEffects(mod.escape.final)) 
Calculation of the marginal and conditional variance explained by the final model
r.squaredGLMM(mod.escape.final)## R2m R2c
## [1,] 0.104123 0.4423867marginal: only for the fixed effect. conditional: both fixed and random effect.
Calculation of the 95% confidence intervals
confint(mod.escape.final, level = 0.95, method = "Wald")## 2.5 % 97.5 %
## .sig01 NA NA
## .sigma NA NA
## (Intercept) 3.31390076 3.4881153
## Order.trial.scaled -0.30222514 -0.1607714
## JulianDay.scaled 0.03643371 0.2447981
## ID.Temperature.scaled -0.29116991 -0.1072327
## CL.scaled 0.08541190 0.2485132
## Wind.scaled -0.24629549 -0.1055340Emergence of the turtle after escaping
Random variable
First, we are testing if the identity of the turtle and the sampling site need to be added in the initial model. We are using likelihood ratio tests to test if the addition of these random variables make a significant effect on the initial model. We are using a dummy variable (same value for all the observations) to create a null mixed model to compared with the different combinations of mixed models.
For generalized mixed models, there is no REML option.
Creation of the different mixed models
## Null mixed model
mod.emergence.null <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|dummy), data = df, family = binomial, control=glmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude)
## only site identity as random variable
mod.emergence.dummy.site <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|dummy) + (1|Site), data = df, family = binomial, control=glmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude)
## only turtle identity as random variable
mod.emergence.dummy.ID <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|dummy) + (1|ID), data = df, family = binomial, control=glmerControl(check.nlev.gtr.1="ignore"), na.action=na.exclude)
## Turtle identity without the dummy variable
mod.emergence.ID <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
## Turtle and site identity without the dummy variable
mod.emergence.ID.site <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|Site) + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)Likelihood ratio tests between the mixed models
The less complex model is always first.
anova(mod.emergence.null, mod.emergence.dummy.site)## Data: df
## Models:
## mod.emergence.null: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy)
## mod.emergence.dummy.site: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy) + (1 | Site)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.null 13 1142.9 1206.8 -558.45 1116.9
## mod.emergence.dummy.site 14 1137.2 1206.0 -554.60 1109.2 7.6971 1
## Pr(>Chisq)
## mod.emergence.null
## mod.emergence.dummy.site 0.005531 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod.emergence.null, mod.emergence.dummy.ID)## Data: df
## Models:
## mod.emergence.null: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy)
## mod.emergence.dummy.ID: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | dummy) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.emergence.null 13 1142.9 1206.8 -558.45 1116.9
## mod.emergence.dummy.ID 14 1113.4 1182.1 -542.68 1085.4 31.55 1 1.944e-08
##
## mod.emergence.null
## mod.emergence.dummy.ID ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod.emergence.ID, mod.emergence.ID.site)## Data: df
## Models:
## mod.emergence.ID: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## mod.emergence.ID.site: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | Site) + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod.emergence.ID 13 1111.3 1175.2 -542.65 1085.3
## mod.emergence.ID.site 14 1113.1 1181.9 -542.56 1085.1 0.1714 1 0.6788When combined with turtle identity, site identity is not significant. Thus, site identity was not kept in the initial model.
Predictor variables
We are selecting the final model with the backward selection procedure. At each step, we deleted the variable with the higher p value. The deletion of each predictor variable is confirmed with a likelihood ratio test. We need to create a new dataset at each step to have only the rows with complete observations for all the preditors. It allow us to perform the likelihood ratio test that does not run between two models with not the same number of observations.
mod.emergence.full <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
summary(mod.emergence.full)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Sex + Tank.Site + Year + JulianDay.scaled +
## Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1111.3 1175.1 -542.6 1085.3 991
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0712 -0.4348 -0.2919 0.6191 4.2215
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 2.167 1.472
## Number of obs: 1004, groups: ID, 684
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.08832 0.25334 -4.296 1.74e-05 ***
## House.300.scaled -0.24078 0.13427 -1.793 0.072940 .
## Distance.channel.scaled -0.26125 0.13935 -1.875 0.060821 .
## Vessels.activities.2019.scaled 0.30458 0.12788 2.382 0.017224 *
## Order.trial.scaled 0.31098 0.11169 2.784 0.005363 **
## SexM -0.07172 0.24168 -0.297 0.766639
## Tank.SiteT -1.94050 0.41256 -4.704 2.56e-06 ***
## Year2020 0.26499 0.23502 1.128 0.259509
## JulianDay.scaled 0.57255 0.17220 3.325 0.000885 ***
## Hour.Platform.scaled 0.05927 0.12331 0.481 0.630738
## ID.Temperature.scaled 0.05879 0.13573 0.433 0.664916
## CL.scaled -0.34771 0.13113 -2.652 0.008009 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. Dstn.. V..201 Ordr.. SexM Tnk.ST Yr2020 JlnDy.
## Hs.300.scld 0.054
## Dstnc.chnn. 0.247 0.169
## Vssl..2019. -0.227 -0.156 -0.478
## Ordr.trl.sc 0.125 -0.159 -0.108 0.094
## SexM -0.594 0.086 -0.007 -0.002 -0.056
## Tank.SiteT 0.028 0.055 -0.121 -0.070 -0.060 0.000
## Year2020 -0.505 -0.075 -0.230 0.211 -0.180 0.050 -0.121
## JulnDy.scld -0.083 0.233 0.222 -0.199 -0.141 0.044 -0.501 0.041
## Hr.Pltfrm.s 0.160 -0.048 0.023 -0.020 0.221 -0.004 -0.479 -0.015 0.064
## ID.Tmprtr.s 0.001 0.138 -0.136 -0.112 0.037 0.022 0.329 -0.259 -0.420
## CL.scaled -0.084 -0.001 0.028 -0.067 -0.076 0.350 0.123 0.024 -0.003
## Hr.Pl. ID.Tm.
## Hs.300.scld
## Dstnc.chnn.
## Vssl..2019.
## Ordr.trl.sc
## SexM
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.Pltfrm.s
## ID.Tmprtr.s -0.074
## CL.scaled 0.054 -0.034We will delete sex (Sex)
mod.emergence.1 <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$House.300.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Sex, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Sex + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.1.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.1.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.1.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Sex + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.1.adjust 12 1109.4 1168.3 -542.69 1085.4
## mod.full.emergence.adjust 13 1111.3 1175.2 -542.65 1085.3 0.0886 1
## Pr(>Chisq)
## mod.emergence.1.adjust
## mod.full.emergence.adjust 0.7659summary(mod.emergence.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + Year + JulianDay.scaled +
## Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1136.4 1195.6 -556.2 1112.4 1017
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0741 -0.4311 -0.2848 0.6068 4.2163
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 2.298 1.516
## Number of obs: 1029, groups: ID, 699
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.16434 0.20830 -5.590 2.27e-08 ***
## House.300.scaled -0.24095 0.13407 -1.797 0.072293 .
## Distance.channel.scaled -0.25319 0.13951 -1.815 0.069548 .
## Vessels.activities.2019.scaled 0.30594 0.12882 2.375 0.017551 *
## Order.trial.scaled 0.29618 0.11165 2.653 0.007986 **
## Tank.SiteT -1.99644 0.41547 -4.805 1.55e-06 ***
## Year2020 0.29694 0.23583 1.259 0.207974
## JulianDay.scaled 0.58897 0.17388 3.387 0.000706 ***
## Hour.Platform.scaled 0.07558 0.12413 0.609 0.542617
## ID.Temperature.scaled 0.08207 0.13544 0.606 0.544519
## CL.scaled -0.31465 0.11541 -2.726 0.006405 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. Dstn.. V..201 Ordr.. Tnk.ST Yr2020 JlnDy. Hr.Pl.
## Hs.300.scld 0.136
## Dstnc.chnn. 0.300 0.179
## Vssl..2019. -0.287 -0.153 -0.483
## Ordr.trl.sc 0.110 -0.147 -0.109 0.087
## Tank.SiteT 0.062 0.059 -0.114 -0.086 -0.055
## Year2020 -0.596 -0.094 -0.233 0.213 -0.175 -0.136
## JulnDy.scld -0.083 0.234 0.227 -0.193 -0.137 -0.503 0.045
## Hr.Pltfrm.s 0.186 -0.048 0.020 -0.014 0.226 -0.488 -0.008 0.071
## ID.Tmprtr.s 0.007 0.130 -0.140 -0.104 0.041 0.324 -0.259 -0.413 -0.079
## CL.scaled 0.180 -0.010 0.025 -0.070 -0.068 0.164 0.031 -0.035 0.030
## ID.Tm.
## Hs.300.scld
## Dstnc.chnn.
## Vssl..2019.
## Ordr.trl.sc
## Tank.SiteT
## Year2020
## JulnDy.scld
## Hr.Pltfrm.s
## ID.Tmprtr.s
## CL.scaled -0.032We will delete the turtle temperature (ID.Temperature.scaled)
mod.emergence.2 <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$House.300.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$ID.Temperature.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.2.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.2.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.2.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + ID.Temperature.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.2.adjust 11 1134.7 1189.0 -556.36 1112.7
## mod.full.emergence.adjust 12 1136.4 1195.6 -556.18 1112.4 0.3658 1
## Pr(>Chisq)
## mod.emergence.2.adjust
## mod.full.emergence.adjust 0.5453summary(mod.emergence.2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + Year + JulianDay.scaled +
## Hour.Platform.scaled + CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1141.1 1195.5 -559.5 1119.1 1030
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1399 -0.4056 -0.2637 0.5823 4.0785
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 2.794 1.672
## Number of obs: 1041, groups: ID, 701
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.23341 0.22847 -5.399 6.71e-08 ***
## House.300.scaled -0.24829 0.13735 -1.808 0.07065 .
## Distance.channel.scaled -0.24398 0.14361 -1.699 0.08932 .
## Vessels.activities.2019.scaled 0.32562 0.13592 2.396 0.01659 *
## Order.trial.scaled 0.28855 0.11556 2.497 0.01253 *
## Tank.SiteT -2.17897 0.42230 -5.160 2.47e-07 ***
## Year2020 0.34784 0.23886 1.456 0.14532
## JulianDay.scaled 0.68081 0.16868 4.036 5.43e-05 ***
## Hour.Platform.scaled 0.08016 0.12783 0.627 0.53062
## CL.scaled -0.32716 0.12254 -2.670 0.00759 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. Dstn.. V..201 Ordr.. Tnk.ST Yr2020 JlnDy. Hr.Pl.
## Hs.300.scld 0.139
## Dstnc.chnn. 0.305 0.196
## Vssl..2019. -0.307 -0.137 -0.508
## Ordr.trl.sc 0.106 -0.141 -0.097 0.085
## Tank.SiteT 0.138 0.016 -0.059 -0.076 -0.056
## Year2020 -0.597 -0.066 -0.287 0.199 -0.174 -0.064
## JulnDy.scld -0.140 0.294 0.168 -0.228 -0.123 -0.473 -0.069
## Hr.Pltfrm.s 0.171 -0.029 0.015 -0.027 0.228 -0.480 -0.031 0.064
## CL.scaled 0.219 -0.004 0.025 -0.087 -0.064 0.213 0.014 -0.089 0.024We will delete the hour of the test (Hour.Platform.scaled)
mod.emergence.3 <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$House.300.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$Hour.Platform.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.3.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.3.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.3.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + Hour.Platform.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.3.adjust 10 1139.5 1189.0 -559.74 1119.5
## mod.full.emergence.adjust 11 1141.1 1195.5 -559.54 1119.1 0.3936 1
## Pr(>Chisq)
## mod.emergence.3.adjust
## mod.full.emergence.adjust 0.5304summary(mod.emergence.3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + Year + JulianDay.scaled +
## CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1187.9 1237.7 -584.0 1167.9 1061
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8830 -0.4206 -0.2687 0.6047 3.6654
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 2.872 1.695
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2293 0.2087 -5.889 3.88e-09 ***
## House.300.scaled -0.1979 0.1340 -1.477 0.1397
## Distance.channel.scaled -0.2311 0.1406 -1.644 0.1001
## Vessels.activities.2019.scaled 0.3024 0.1339 2.259 0.0239 *
## Order.trial.scaled 0.1968 0.1073 1.835 0.0665 .
## Tank.SiteT -2.0564 0.3685 -5.580 2.40e-08 ***
## Year2020 0.2985 0.2310 1.293 0.1961
## JulianDay.scaled 0.6907 0.1689 4.090 4.30e-05 ***
## CL.scaled -0.2914 0.1208 -2.412 0.0159 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. Dstn.. V..201 Ordr.. Tnk.ST Yr2020 JlnDy.
## Hs.300.scld 0.097
## Dstnc.chnn. 0.248 0.178
## Vssl..2019. -0.249 -0.112 -0.488
## Ordr.trl.sc 0.047 -0.145 -0.150 0.118
## Tank.SiteT 0.276 -0.006 -0.060 -0.095 0.142
## Year2020 -0.544 -0.034 -0.246 0.151 -0.142 -0.080
## JulnDy.scld -0.204 0.302 0.157 -0.223 -0.217 -0.503 -0.063
## CL.scaled 0.191 -0.019 0.004 -0.062 -0.050 0.238 0.049 -0.088We will delete year (Year)
mod.emergence.4 <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$House.300.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$Year, df$JulianDay.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.4.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.4.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.4.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + Year + JulianDay.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.4.adjust 9 1187.6 1232.4 -584.79 1169.6
## mod.full.emergence.adjust 10 1187.9 1237.7 -583.96 1167.9 1.6459 1
## Pr(>Chisq)
## mod.emergence.4.adjust
## mod.full.emergence.adjust 0.1995summary(mod.emergence.4)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled +
## (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1187.6 1232.4 -584.8 1169.6 1062
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8516 -0.4131 -0.2652 0.6022 3.3689
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.052 1.747
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1066 0.1835 -6.030 1.64e-09 ***
## House.300.scaled -0.1935 0.1351 -1.433 0.1520
## Distance.channel.scaled -0.1881 0.1381 -1.361 0.1734
## Vessels.activities.2019.scaled 0.2814 0.1346 2.090 0.0366 *
## Order.trial.scaled 0.2175 0.1076 2.022 0.0432 *
## Tank.SiteT -2.0603 0.3769 -5.466 4.60e-08 ***
## JulianDay.scaled 0.7197 0.1737 4.145 3.40e-05 ***
## CL.scaled -0.3067 0.1232 -2.490 0.0128 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. Dstn.. V..201 Ordr.. Tnk.ST JlnDy.
## Hs.300.scld 0.094
## Dstnc.chnn. 0.136 0.173
## Vssl..2019. -0.204 -0.108 -0.470
## Ordr.trl.sc -0.041 -0.148 -0.192 0.143
## Tank.SiteT 0.309 -0.006 -0.082 -0.084 0.129
## JulnDy.scld -0.315 0.291 0.142 -0.206 -0.221 -0.529
## CL.scaled 0.274 -0.014 0.016 -0.073 -0.046 0.253 -0.101We will delete the distance to the navigation channel (Distance.channel.scaled)
mod.emergence.5 <- glmer(emergence ~ House.300.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$House.300.scaled, df$Distance.channel.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$JulianDay.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.5.adjust <- glmer(emergence ~ House.300.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.5.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.5.adjust: emergence ~ House.300.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ House.300.scaled + Distance.channel.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.5.adjust 8 1187.5 1227.3 -585.75 1171.5
## mod.full.emergence.adjust 9 1187.6 1232.4 -584.79 1169.6 1.9262 1
## Pr(>Chisq)
## mod.emergence.5.adjust
## mod.full.emergence.adjust 0.1652summary(mod.emergence.5)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: emergence ~ House.300.scaled + Vessels.activities.2019.scaled +
## Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled +
## (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1187.5 1227.3 -585.7 1171.5 1063
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7916 -0.4051 -0.2658 0.5978 3.2281
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.167 1.78
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0994 0.1878 -5.853 4.83e-09 ***
## House.300.scaled -0.1601 0.1339 -1.196 0.2315
## Vessels.activities.2019.scaled 0.1958 0.1190 1.646 0.0998 .
## Order.trial.scaled 0.1901 0.1066 1.784 0.0745 .
## Tank.SiteT -2.1470 0.3806 -5.642 1.69e-08 ***
## JulianDay.scaled 0.7748 0.1753 4.420 9.86e-06 ***
## CL.scaled -0.3132 0.1251 -2.504 0.0123 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H.300. V..201 Ordr.. Tnk.ST JlnDy.
## Hs.300.scld 0.071
## Vssl..2019. -0.155 -0.022
## Ordr.trl.sc -0.023 -0.116 0.061
## Tank.SiteT 0.342 0.008 -0.131 0.112
## JulnDy.scld -0.362 0.267 -0.152 -0.192 -0.536
## CL.scaled 0.289 -0.014 -0.075 -0.048 0.265 -0.117We will delete the number of house with access to the canal within 300m (House.300.scaled)
mod.emergence.6 <- glmer(emergence ~ Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$House.300.scaled, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$JulianDay.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ House.300.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.6.adjust <- glmer(emergence ~ Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.6.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.6.adjust: emergence ~ Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ House.300.scaled + Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.6.adjust 7 1186.9 1221.8 -586.47 1172.9
## mod.full.emergence.adjust 8 1187.5 1227.3 -585.75 1171.5 1.4458 1
## Pr(>Chisq)
## mod.emergence.6.adjust
## mod.full.emergence.adjust 0.2292summary(mod.emergence.6)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: emergence ~ Vessels.activities.2019.scaled + Order.trial.scaled +
## Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1186.9 1221.8 -586.5 1172.9 1064
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7912 -0.4000 -0.2686 0.5881 3.1179
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.307 1.819
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1053 0.1941 -5.694 1.24e-08 ***
## Vessels.activities.2019.scaled 0.1959 0.1209 1.620 0.1051
## Order.trial.scaled 0.1764 0.1066 1.654 0.0981 .
## Tank.SiteT -2.1788 0.3902 -5.584 2.35e-08 ***
## JulianDay.scaled 0.8474 0.1731 4.896 9.79e-07 ***
## CL.scaled -0.3218 0.1268 -2.538 0.0112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) V..201 Ordr.. Tnk.ST JlnDy.
## Vssl..2019. -0.161
## Ordr.trl.sc -0.026 0.062
## Tank.SiteT 0.368 -0.138 0.106
## JulnDy.scld -0.416 -0.143 -0.157 -0.571
## CL.scaled 0.300 -0.079 -0.055 0.273 -0.128We will delete the mean daily number of vessel crossings (Vessels.activities.2019.scaled)
mod.emergence.7 <- glmer(emergence ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$Vessels.activities.2019.scaled, df$Order.trial.scaled, df$Tank.Site, df$JulianDay.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.7.adjust <- glmer(emergence ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.7.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.7.adjust: emergence ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ Vessels.activities.2019.scaled + Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.7.adjust 6 1187.7 1217.5 -587.83 1175.7
## mod.full.emergence.adjust 7 1186.9 1221.8 -586.47 1172.9 2.7133 1
## Pr(>Chisq)
## mod.emergence.7.adjust
## mod.full.emergence.adjust 0.09951 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod.emergence.7)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: emergence ~ Order.trial.scaled + Tank.Site + JulianDay.scaled +
## CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1187.7 1217.5 -587.8 1175.7 1065
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7853 -0.4036 -0.2656 0.5732 2.9896
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.34 1.827
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1039 0.1942 -5.685 1.31e-08 ***
## Order.trial.scaled 0.1674 0.1067 1.570 0.1165
## Tank.SiteT -2.1776 0.3909 -5.571 2.53e-08 ***
## JulianDay.scaled 0.9245 0.1733 5.335 9.55e-08 ***
## CL.scaled -0.3214 0.1268 -2.535 0.0112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ordr.. Tnk.ST JlnDy.
## Ordr.trl.sc -0.022
## Tank.SiteT 0.369 0.112
## JulnDy.scld -0.460 -0.145 -0.611
## CL.scaled 0.295 -0.053 0.269 -0.147We will delete the order of the trial (Order.trial.scaled)
mod.emergence.8 <- glmer(emergence ~ Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
df.adjust <- df[complete.cases(df$emergence, df$Order.trial.scaled, df$Tank.Site, df$JulianDay.scaled, df$CL.scaled),]
mod.full.emergence.adjust <- glmer(emergence ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
mod.emergence.8.adjust <- glmer(emergence ~ Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df.adjust, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
anova(mod.emergence.8.adjust, mod.full.emergence.adjust) # less complex first## Data: df.adjust
## Models:
## mod.emergence.8.adjust: emergence ~ Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## mod.full.emergence.adjust: emergence ~ Order.trial.scaled + Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## npar AIC BIC logLik deviance Chisq Df
## mod.emergence.8.adjust 5 1188.1 1213.0 -589.04 1178.1
## mod.full.emergence.adjust 6 1187.7 1217.5 -587.83 1175.7 2.4288 1
## Pr(>Chisq)
## mod.emergence.8.adjust
## mod.full.emergence.adjust 0.1191summary(mod.emergence.8)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: emergence ~ Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1188.1 1213.0 -589.0 1178.1 1066
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4438 -0.4014 -0.2571 0.5682 2.9466
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.563 1.888
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1353 0.2078 -5.464 4.67e-08 ***
## Tank.SiteT -2.3175 0.4170 -5.558 2.73e-08 ***
## JulianDay.scaled 0.9938 0.1816 5.473 4.42e-08 ***
## CL.scaled -0.3230 0.1306 -2.473 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tnk.ST JlnDy.
## Tank.SiteT 0.446
## JulnDy.scld -0.514 -0.641
## CL.scaled 0.325 0.311 -0.192Final model
Summary statistics
mod.emergence.final <- glmer(emergence ~ Tank.Site + JulianDay.scaled + CL.scaled + (1|ID), data = df, family = binomial, control = glmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')), na.action=na.exclude)
summary(mod.emergence.final)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: emergence ~ Tank.Site + JulianDay.scaled + CL.scaled + (1 | ID)
## Data: df
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
##
## AIC BIC logLik deviance df.resid
## 1188.1 1213.0 -589.0 1178.1 1066
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4438 -0.4014 -0.2571 0.5682 2.9466
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.563 1.888
## Number of obs: 1071, groups: ID, 704
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1353 0.2078 -5.464 4.67e-08 ***
## Tank.SiteT -2.3175 0.4170 -5.558 2.73e-08 ***
## JulianDay.scaled 0.9938 0.1816 5.473 4.42e-08 ***
## CL.scaled -0.3230 0.1306 -2.473 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tnk.ST JlnDy.
## Tank.SiteT 0.446
## JulnDy.scld -0.514 -0.641
## CL.scaled 0.325 0.311 -0.192Quick visualization of the preditor effects
plot(allEffects(mod.emergence.final)) 
Calculation of the marginal and conditional variance explained by the final model
r.squaredGLMM(mod.emergence.final)## R2m R2c
## theoretical 0.1440800 0.5891187
## delta 0.1206006 0.4931154marginal: only for the fixed effect. conditional: both fixed and random effect.
The delta method can be used with all distributions and link functions.
Calculation of the 95% confidence intervals
confint(mod.emergence.final, level = 0.95, method = "Wald")## 2.5 % 97.5 %
## .sig01 NA NA
## (Intercept) -1.5426261 -0.7280630
## Tank.SiteT -3.1347722 -1.5001995
## JulianDay.scaled 0.6379106 1.3496770
## CL.scaled -0.5790026 -0.0670086