prey_capture_Gasteracantha
Fabian C. Salgado-Roa
22/19/2023
Introduction
Load the libraries and standard error function
library(tidyverse)
library(glmmTMB)
library(car)
library(DHARMa)Standard error function
se <- function(x) sd(x)/sqrt(length(x))Load the datasets
Total dataset
tmp_noNA<-read_csv("data_gasteracantha.csv",col_names=T)Treat time as a categorical variable
tmp_noNA$Time<-as.character(tmp_noNA$Time)Let’s create subset excluding the empty web observations. This will facilitate the comparisons between colour morphs
tmp_ind<-tmp_noNA %>% filter(colour!="empty")Because the measures are repeated for each individual, lets create a subset with a single observation to explore association between colour morph and the variables measured
unique_ind <- tmp_noNA %>% distinct(code, .keep_all=T)
unique_noEmpty<-unique_ind %>% filter(colour!="empty")Association between the predictors
Check the normality of the variables per colour category (including empty web when possible)
#Web height
hist(unique_ind$Web_height)
qqPlot(unique_ind$Web_height,groups=factor(unique_ind$colour))
#opisthosoma width
hist(unique_noEmpty$opisto_w)
qqPlot(unique_noEmpty$opisto_w)
#web Area
hist(log(unique_ind$web_area+10))
qqPlot(log(unique_ind$web_area+10),groups=factor(unique_ind$colour))
It seems that all the variables are following a normal distribution in every colour category
Linear models
Web height
Run the linear model
mod_webHeight<-lm(Web_height~colour,data=unique_ind)
#plot(mod_webHeight)Check the significance of the predictor
summary(mod_webHeight)##
## Call:
## lm(formula = Web_height ~ colour, data = unique_ind)
##
## Residuals:
## Min 1Q Median 3Q Max
## -90.400 -29.238 0.762 30.560 93.790
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 130.400 8.241 15.824 <2e-16 ***
## colourempty -14.352 12.197 -1.177 0.241
## colourwhite -4.190 9.762 -0.429 0.668
## colouryellow -1.162 10.408 -0.112 0.911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.2 on 146 degrees of freedom
## Multiple R-squared: 0.01192, Adjusted R-squared: -0.008382
## F-statistic: 0.5871 on 3 and 146 DF, p-value: 0.6244anova(mod_webHeight)## Analysis of Variance Table
##
## Response: Web_height
## Df Sum Sq Mean Sq F value Pr(>F)
## colour 3 2990 996.83 0.5871 0.6244
## Residuals 146 247871 1697.75Opisthosoma width
Run the linear model
mod_opistoW<-lm(opisto_w~colour,data=unique_noEmpty)
#plot(mod_opistoW)Check the significance of the predictor
summary(mod_opistoW)##
## Call:
## lm(formula = opisto_w ~ colour, data = unique_noEmpty)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2229 -0.8017 0.0826 0.8767 2.2771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.6229 0.2640 36.444 <2e-16 ***
## colourwhite 0.5916 0.3245 1.823 0.0712 .
## colouryellow 0.8945 0.3428 2.609 0.0104 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.294 on 103 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.0625, Adjusted R-squared: 0.04429
## F-statistic: 3.433 on 2 and 103 DF, p-value: 0.03603anova(mod_opistoW)## Analysis of Variance Table
##
## Response: opisto_w
## Df Sum Sq Mean Sq F value Pr(>F)
## colour 2 11.489 5.7445 3.4331 0.03603 *
## Residuals 103 172.350 1.6733
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Opisthosoma width seems to differ between colour morphs
Web area
Run the linear model
mod_webArea<-lm(sqrt(web_area+10)~colour,data=unique_ind)Check the assumptions of the linear model
plot(mod_webArea)
shapiro.test(mod_webArea$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_webArea$residuals
## W = 0.98884, p-value = 0.2767everything seems ok
Check the significance of the predictor
summary(mod_webArea)##
## Call:
## lm(formula = sqrt(web_area + 10) ~ colour, data = unique_ind)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.756 -2.785 -0.096 3.053 9.800
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.34362 0.81604 23.704 <2e-16 ***
## colourempty -0.09935 1.20775 -0.082 0.935
## colourwhite -0.63583 0.96666 -0.658 0.512
## colouryellow 0.29297 1.03068 0.284 0.777
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.08 on 146 degrees of freedom
## Multiple R-squared: 0.009452, Adjusted R-squared: -0.0109
## F-statistic: 0.4644 on 3 and 146 DF, p-value: 0.7076anova(mod_webArea)## Analysis of Variance Table
##
## Response: sqrt(web_area + 10)
## Df Sum Sq Mean Sq F value Pr(>F)
## colour 3 23.19 7.7307 0.4644 0.7076
## Residuals 146 2430.59 16.6478get the summary value per colour morph
unique_ind %>% group_by(colour) %>% summarize(num(mean(web_area),digits=2),se(web_area))## # A tibble: 4 x 3
## colour `num(mean(web_area), digits = 2)` `se(web_area)`
## <chr> <num:.2!> <dbl>
## 1 black 372.07 22.1
## 2 empty 375.58 34.7
## 3 white 360.56 22.2
## 4 yellow 390.76 23.6unique_ind %>% group_by(colour) %>% summarize(num(mean(Web_height),digits=2),se(Web_height))## # A tibble: 4 x 3
## colour `num(mean(Web_height), digits = 2)` `se(Web_height)`
## <chr> <num:.2!> <dbl>
## 1 black 130.40 9.33
## 2 empty 116.05 9.49
## 3 white 126.21 5.45
## 4 yellow 129.24 5.15unique_ind %>% filter(colour!="empty") %>% drop_na(opisto_w) %>% group_by(colour) %>% summarize(num(mean(opisto_w),digits=2),se(opisto_w))## # A tibble: 3 x 3
## colour `num(mean(opisto_w), digits = 2)` `se(opisto_w)`
## <chr> <num:.2!> <dbl>
## 1 black 9.62 0.333
## 2 white 10.21 0.170
## 3 yellow 10.52 0.203Because we found that the colour morphs differ in their opisthosoma width, we explored if this predictor is associated with prey capture
model with negative binomial distribution
Run the model
capture_opistoW<-glmmTMB(formula=total_preys ~ opisto_w + (1|code) + (1|check),
family="nbinom2",REML=T,data=tmp_ind)Check the uniformity, dispersion and zero inflation of the model
simulationOutput <- simulateResiduals(fittedModel = capture_opistoW, n = 1000)
testDispersion(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.82239, p-value = 0.528
## alternative hypothesis: two.sidedtestUniformity(simulationOutput = simulationOutput, alternative ="two.sided")
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.023038, p-value = 0.3784
## alternative hypothesis: two-sidedtestZeroInflation(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0024, p-value = 0.946
## alternative hypothesis: two.sidedeverything seems ok
Check the significance of the predictor
Anova(capture_opistoW)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: total_preys
## Chisq Df Pr(>Chisq)
## opisto_w 0.328 1 0.5668Because width is not related with prey capture, which is the variable we are interested in, we decided to remove it for the subsequent analyses
Web area and web height present collinearity
Run the model
mod_webArea_Heigth<-lm(sqrt(web_area+10)~Web_height,data=unique_ind)Check the assumptions of the linear model
plot(mod_webArea_Heigth)
shapiro.test(mod_webArea_Heigth$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_webArea_Heigth$residuals
## W = 0.99232, p-value = 0.6011everything seems ok
Check the significance of the predictor
summary(mod_webArea_Heigth)##
## Call:
## lm(formula = sqrt(web_area + 10) ~ Web_height, data = unique_ind)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6118 -2.7954 -0.0652 2.7971 9.7547
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.982515 1.017724 14.722 < 2e-16 ***
## Web_height 0.032980 0.007664 4.303 3.05e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.839 on 148 degrees of freedom
## Multiple R-squared: 0.1112, Adjusted R-squared: 0.1052
## F-statistic: 18.52 on 1 and 148 DF, p-value: 3.048e-05Anova(mod_webArea_Heigth)## Anova Table (Type II tests)
##
## Response: sqrt(web_area + 10)
## Sum Sq Df F value Pr(>F)
## Web_height 272.85 1 18.516 3.048e-05 ***
## Residuals 2180.93 148
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Light differences between sides and colour morphs
Light measures are bimodal, in consequences we decided to this variable as binary (high or low). This categorization can be explore
paleta<-c("#636363","#f0f0f0","#ffeda0")
ggplot(data=tmp_noNA, aes(x=log(luxes*100), group=colour, fill=colour)) +
geom_density(adjust=1.5, alpha=.4) +
scale_fill_manual(values=paleta)+
scale_x_continuous(limits=c(0,12))+
theme_classic()+labs(y="Density",x="Log(lux*100)")## Warning: Removed 274 rows containing non-finite values (`stat_density()`).
We explore if there were differences between sides in the exposure to a certain light condition
Run the model
light_side<-glmmTMB(formula=b_luxes ~ Side+ (1|code)+(1|day),family="binomial",data=tmp_ind)Check the significance of the predictor
summary(light_side)## Family: binomial ( logit )
## Formula: b_luxes ~ Side + (1 | code) + (1 | day)
## Data: tmp_ind
##
## AIC BIC logLik deviance df.resid
## 1516.6 1538.5 -754.3 1508.6 1774
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## code (Intercept) 3.9175 1.9793
## day (Intercept) 0.4532 0.6732
## Number of obs: 1778, groups: code, 129; day, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.8531 0.3745 -7.618 2.57e-14 ***
## Sideventral 1.5409 0.1504 10.245 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(light_side)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: b_luxes
## Chisq Df Pr(>Chisq)
## Side 104.95 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We found that the ventral side is more often expose to conditions of high luminosity. Let’s visualize the results
#Predict the values based on the model
tmp_ind$light_side<-predict(light_side,type="response")
#Colours
paleta_side<-c("#B84C7A","#BCA7D2")
#Plot
ggplot(data = tmp_ind,aes(x=Side,y=light_side, fill=Side))+
scale_y_continuous(limits=c(0,1))+
scale_fill_manual(values=paleta_side)+
scale_colour_manual(values=paleta_side)+
stat_summary(fun = mean,aes(color = Side,group=Side),fun.min = function(x) mean(x) - (2*se(x)),fun.max = function(x) mean(x)+(2*se(x)),geom = "pointrange", position = position_jitterdodge(jitter.width=0.15), size=1.5)+
theme_classic()+labs(x="Spider side",y="Light environment")
We also tested if the colour morphs differ in which side is expose to the different light environments
Run the model with yellow and ventral as fixed levels
tmp_ind$colour<-relevel(as.factor(tmp_ind$colour),"yellow")
tmp_ind$Side<-relevel(as.factor(tmp_ind$Side),"ventral")
colour_light1<-glmmTMB(formula=b_luxes ~ Side*colour + (1|code) + (1|day) ,family="binomial",data=tmp_ind)Check the significance of the predictor
summary(colour_light1)## Family: binomial ( logit )
## Formula: b_luxes ~ Side * colour + (1 | code) + (1 | day)
## Data: tmp_ind
##
## AIC BIC logLik deviance df.resid
## 1485.9 1529.7 -734.9 1469.9 1770
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## code (Intercept) 3.9748 1.994
## day (Intercept) 0.4597 0.678
## Number of obs: 1778, groups: code, 129; day, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5977 0.4540 -1.317 0.18800
## Sidedorsal -2.7857 0.2873 -9.696 < 2e-16 ***
## colourblack -1.4535 0.6207 -2.342 0.01919 *
## colourwhite -0.9977 0.4570 -2.183 0.02904 *
## Sidedorsal:colourblack 1.4252 0.4776 2.984 0.00284 **
## Sidedorsal:colourwhite 1.9959 0.3488 5.722 1.05e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(colour_light1,type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: b_luxes
## Chisq Df Pr(>Chisq)
## (Intercept) 1.7332 1 0.188
## Side 94.0183 1 < 2.2e-16 ***
## colour 7.1512 2 0.028 *
## Side:colour 32.8057 2 7.522e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Run the model with black and dorsal as fixed levels
tmp_ind$colour<-relevel(as.factor(tmp_ind$colour),"black")
tmp_ind$Side<-relevel(as.factor(tmp_ind$Side),"dorsal")
colour_light2<-glmmTMB(formula=b_luxes ~ Side*colour + (1|code) + (1|day),
family="binomial",data=tmp_ind)Check the significance of the predictors
summary(colour_light2)## Family: binomial ( logit )
## Formula: b_luxes ~ Side * colour + (1 | code) + (1 | day)
## Data: tmp_ind
##
## AIC BIC logLik deviance df.resid
## 1485.9 1529.7 -734.9 1469.9 1770
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## code (Intercept) 3.9749 1.994
## day (Intercept) 0.4597 0.678
## Number of obs: 1778, groups: code, 129; day, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.41165 0.63944 -5.335 9.53e-08 ***
## Sideventral 1.36046 0.38597 3.525 0.000424 ***
## colouryellow 0.02826 0.68811 0.041 0.967236
## colourwhite 1.02645 0.63868 1.607 0.108024
## Sideventral:colouryellow 1.42520 0.47755 2.984 0.002841 **
## Sideventral:colourwhite -0.57068 0.43453 -1.313 0.189072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(colour_light2,type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: b_luxes
## Chisq Df Pr(>Chisq)
## (Intercept) 28.4665 1 9.534e-08 ***
## Side 12.4243 1 0.0004238 ***
## colour 5.1592 2 0.0758037 .
## Side:colour 32.8061 2 7.520e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Let’s visualize the results
#Predict the values based on the model
tmp_ind$predicted_light<-predict(colour_light1,type="response")
#Colours
tmp_ind$colour<-relevel(as.factor(tmp_ind$colour),"white")
tmp_ind$Side<-relevel(as.factor(tmp_ind$Side),"dorsal")
paleta2<-c("#f0f0f0","#636363","#ffeda0")
#Plot
ggplot(data = tmp_ind,aes(x=Side,y=predicted_light, fill=colour))+
scale_y_continuous(limits=c(0,1))+labs(x="Spider side",y="Light environment")+
scale_fill_manual(values=paleta2)+
scale_colour_manual(values=paleta2)+
stat_summary(fun = mean,aes(color = colour,group=colour),fun.min = function(x) mean(x) - (2*se(x)),fun.max = function(x) mean(x)+(2*se(x)),geom = "pointrange", shape=22,col="black",position = position_jitterdodge(jitter.width=0.15), size=1)+
theme_classic()
Side filtering and analyses
we compared dorsal and ventral light measures of each individual, where the higher measure was considered as open space and the lower measure as vegetation. We did this for everyday of observation. The following code will create a new dataset.
unique_ind2<-data.frame()
for(d in unique(tmp_ind$day)){
sub<-tmp_ind %>% filter(day==d & Time==6.3)
unique_ind2<-rbind(unique_ind2,sub)
}
side_data<-tibble(code=character(),side=character(),open_space=numeric(),day=numeric())
for(d in unique(unique_ind2$day)){
filtered_day<-unique_ind2 %>% filter(day==d)
tmp<-tibble(code=character(),side=character(),open_space=numeric(),day=numeric())
for(ind in unique(filtered_day$code)) {
dorsal<-filtered_day %>% filter(code==ind & Side=="dorsal")
ventral<-filtered_day %>% filter(code==ind & Side=="ventral")
if(nrow(ventral)==0 | nrow(dorsal)==0){next}
else if(ventral$luxes>dorsal$luxes){
ventral<-tibble(code=ind,side="ventral",open_space=1,day=d)
dorsal<-tibble(code=ind,side="dorsal",open_space=0,day=d)
tmp<-rbind(tmp,ventral,dorsal)
} else if(ventral$luxes<dorsal$luxes){
ventral<-tibble(code=ind,side="ventral",open_space=0,day=d)
dorsal<-tibble(code=ind,side="dorsal",open_space=1,day=d)
tmp<-rbind(tmp,ventral,dorsal)
}
}
side_data<-rbind(side_data,tmp)
}Let’s evaluate if the new categorization is independent of how the spider position its body towards the different light conditions (High or Low luminosity)
test_association<-chisq.test(factor(unique_ind2$t_luxes),factor(unique_ind2$site_exp))
test_association##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: factor(unique_ind2$t_luxes) and factor(unique_ind2$site_exp)
## X-squared = 2.1355, df = 1, p-value = 0.1439Because we did not find association, let’s run an independent model
open_side<-glmmTMB(formula=open_space ~ side+ (1|day) + (1|code),family="binomial",data=side_data) #test which side is mostCheck the model’s fit
simulationOutput <- simulateResiduals(fittedModel = open_side, n = 1000)
testDispersion(simulationOutput = simulationOutput, alternative ="two.sided")
testUniformity(simulationOutput = simulationOutput, alternative ="two.sided")
Check the significance of the predictor
summary(open_side)## Family: binomial ( logit )
## Formula: open_space ~ side + (1 | day) + (1 | code)
## Data: side_data
##
## AIC BIC logLik deviance df.resid
## 601.9 618.4 -297.0 593.9 450
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## day (Intercept) 1.099e-11 3.315e-06
## code (Intercept) 4.651e-10 2.157e-05
## Number of obs: 454, groups: day, 6; code, 129
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5700 0.1382 -4.125 3.70e-05 ***
## sideventral 1.1400 0.1954 5.834 5.41e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(open_side)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: open_space
## Chisq Df Pr(>Chisq)
## side 34.038 1 5.406e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Capture models
Effect of the presence of the colour morphs on the web
We explored the effect of the presence of the presence of the spiders on the web. We also included the presence of silk decorations as a predictor
tmp_noNA$colour<-relevel(as.factor(tmp_noNA$colour),"empty")
tmp_noNA$Time<-relevel(as.factor(tmp_noNA$Time),"6.3")
with_control_model<-glmmTMB(formula=total_preys ~ colour+sqrt(web_area+10)+presence_decorations + check + (1|code) + (1|day/Time),family="poisson",REML=F,data=tmp_noNA) Check the uniformity, dispersion and zero inflation of the model
simulationOutput <- simulateResiduals(fittedModel = with_control_model, n = 10000)
testDispersion(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1036, p-value = 0.5726
## alternative hypothesis: two.sidedtestUniformity(simulationOutput = simulationOutput, alternative ="two.sided")
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055171, p-value = 7.514e-06
## alternative hypothesis: two-sidedtestZeroInflation(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0595, p-value = 0.1614
## alternative hypothesis: two.sidedBecause there is not uniformity of the residuals when using a poisson distribution, we run the model again with a negative binomial distribution
tmp_noNA$colour<-relevel(as.factor(tmp_noNA$colour),"empty")
tmp_noNA$Time<-relevel(as.factor(tmp_noNA$Time),"6.3")
with_control_model<-glmmTMB(formula=total_preys ~ colour+sqrt(web_area+10)+presence_decorations + check + (1|code) + (1|day/Time),family="nbinom2",REML=F,data=tmp_noNA) Check the uniformity, dispersion and zero inflation of the model
simulationOutput <- simulateResiduals(fittedModel = with_control_model, n = 10000)
testDispersion(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.57703, p-value = 0.4786
## alternative hypothesis: two.sidedtestUniformity(simulationOutput = simulationOutput, alternative ="two.sided")
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.02864, p-value = 0.06903
## alternative hypothesis: two-sidedtestZeroInflation(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0122, p-value = 0.7122
## alternative hypothesis: two.sidedeverything seems ok
Check the significance of the predictors
summary(with_control_model)## Family: nbinom2 ( log )
## Formula:
## total_preys ~ colour + sqrt(web_area + 10) + presence_decorations +
## check + (1 | code) + (1 | day/Time)
## Data: tmp_noNA
##
## AIC BIC logLik deviance df.resid
## 2135.1 2202.6 -1055.6 2111.1 2040
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## code (Intercept) 0.1192 0.3452
## Time:day (Intercept) 0.1922 0.4384
## day (Intercept) 0.2696 0.5193
## Number of obs: 2052, groups: code, 150; Time:day, 23; day, 6
##
## Dispersion parameter for nbinom2 family (): 0.339
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.45213 0.48582 -2.989 0.0028 **
## colourblack -1.19340 0.26087 -4.575 4.77e-06 ***
## colourwhite -1.11115 0.21595 -5.145 2.67e-07 ***
## colouryellow -1.18243 0.22258 -5.312 1.08e-07 ***
## sqrt(web_area + 10) 0.03853 0.01995 1.931 0.0534 .
## presence_decorationsyes -0.10909 0.16922 -0.645 0.5192
## checkF 0.00459 0.19537 0.023 0.9813
## checkM -0.18729 0.18275 -1.025 0.3054
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(with_control_model)## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: total_preys
## Chisq Df Pr(>Chisq)
## colour 36.9717 3 4.665e-08 ***
## sqrt(web_area + 10) 3.7301 1 0.05344 .
## presence_decorations 0.4155 1 0.51917
## check 1.4183 2 0.49207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1the presence of any of the colour morphs on the web seems to reduce the capture of prey
paleta2<-c("#648ace","#636363","#f0f0f0","#ffeda0")
#pdf("dorsal.pdf",height=8,width=10)
tmp_noNA$total_response<-predict(with_control_model,type="response")
ggplot(data = tmp_noNA,aes(x=colour,y=total_response))+
geom_point(data = tmp_noNA,aes(x=colour,y=total_preys, fill=colour),shape = 21,size=3, position = position_jitterdodge(jitter.width=0.2,jitter.height=0.1),alpha=0.5)+
scale_fill_manual(values=paleta2)+
geom_violin(aes(x=colour,y=total_response,fill=colour),alpha=0.3,size=0.4)+scale_color_manual(values=paleta2)+
stat_summary(fun = mean,aes(group=colour,col=colour),fun.min = function(x) mean(x) - sd(x),fun.max = function(x) mean(x)+sd(x),geom = "pointrange",linewidth = 2,size=1.5,shape=22,col="black",fill=paleta2)+
theme_classic()+labs(x="Colour",y="Number of prey capture per hour")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Explore the differences on prey capture considering side and colour morph
Run the model with white and ventral as fixed levels
tmp_ind$colour<-relevel(as.factor(tmp_ind$colour),"white")
tmp_ind$Side<-relevel(as.factor(tmp_ind$Side),"ventral")
tmp_ind$Time<-relevel(as.factor(tmp_ind$Time),"6.3")
total_model_poisson<-glmmTMB::glmmTMB(formula=total_preys ~ Side+colour+t_luxes+sqrt(web_area+10)+ check + (1|code) + (1|day/Time),family="poisson",REML=F,data=tmp_ind)Poisson model
total_model_poisson<-glmmTMB::glmmTMB(formula=total_preys ~ Side+colour+t_luxes+sqrt(web_area+10)+ check + (1|code) + (1|day/Time),family="poisson",REML=F,data=tmp_ind)Test the fit of the model and zero inflation
simulationOutput <- simulateResiduals(fittedModel = total_model_poisson, n = 10000)
testDispersion(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.2804, p-value = 0.3496
## alternative hypothesis: two.sidedtestUniformity(simulationOutput = simulationOutput, alternative ="two.sided")
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.047735, p-value = 0.0006053
## alternative hypothesis: two-sidedtestZeroInflation(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.04, p-value = 0.2346
## alternative hypothesis: two.sidedBecause this model has no uniformity; its is necessary to try other model that takes this into account
Negative binomial model
total_model_nb<-glmmTMB::glmmTMB(formula=total_preys ~ colour+Side+t_luxes+sqrt(web_area+10)+ check + (1|code) + (1|day/Time),family="nbinom2",REML=F,data=tmp_ind)Test the fit of the model and zero inflation
simulationOutput <- simulateResiduals(fittedModel = total_model_nb, n = 10000)
testDispersion(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.7817, p-value = 0.8058
## alternative hypothesis: two.sidedtestUniformity(simulationOutput = simulationOutput, alternative ="two.sided")
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.018165, p-value = 0.6004
## alternative hypothesis: two-sidedtestZeroInflation(simulationOutput = simulationOutput, alternative ="two.sided")
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0086, p-value = 0.7794
## alternative hypothesis: two.sidedThis model seems to be ok, let’s now check the significance of the predictors
summary(total_model_nb)## Family: nbinom2 ( log )
## Formula:
## total_preys ~ colour + Side + t_luxes + sqrt(web_area + 10) +
## check + (1 | code) + (1 | day/Time)
## Data: tmp_ind
##
## AIC BIC logLik deviance df.resid
## 1625.4 1691.2 -800.7 1601.4 1766
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## code (Intercept) 0.2037 0.4513
## Time:day (Intercept) 0.2596 0.5095
## day (Intercept) 0.2066 0.4546
## Number of obs: 1778, groups: code, 129; Time:day, 23; day, 6
##
## Dispersion parameter for nbinom2 family (): 0.387
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.10971 0.56554 -5.499 3.83e-08 ***
## colourblack -0.17275 0.25510 -0.677 0.498281
## colouryellow -0.12137 0.20569 -0.590 0.555169
## Sidedorsal 0.18636 0.15710 1.186 0.235528
## t_luxeslow 0.72967 0.22143 3.295 0.000983 ***
## sqrt(web_area + 10) 0.03035 0.02287 1.327 0.184353
## checkF 0.04343 0.23224 0.187 0.851646
## checkM -0.35026 0.20964 -1.671 0.094773 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(total_model_nb,type = 3)## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: total_preys
## Chisq Df Pr(>Chisq)
## (Intercept) 30.2355 1 3.826e-08 ***
## colour 0.6111 2 0.7367028
## Side 1.4072 1 0.2355283
## t_luxes 10.8590 1 0.0009831 ***
## sqrt(web_area + 10) 1.7622 1 0.1843529
## check 3.9086 2 0.1416656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To create the plot for the results, we have to predict the values under the model
#Predict values to plot
tmp_ind$response_nb<-predict(total_model_nb,type="response")With this, we can create two plots, one per side in order to make easy to visualize the results
###Plot dorsal
Create subset
plot_dorsal <- tmp_ind %>% filter(Side=="dorsal")
plot_dorsal$colour<-relevel(as.factor(plot_dorsal$colour),"black")Plot the results
#Colours
paleta2<-c("#636363","#f0f0f0","#ffeda0")
#plot
#pdf("dorsal.pdf",height=8,width=10)
ggplot(data = plot_dorsal,aes(x=t_luxes,y=response_nb, fill=colour))+
geom_point(data = plot_dorsal,aes(x=t_luxes,y=total_preys, fill=colour),shape = 21,size=4, position = position_jitterdodge(jitter.width=0.15,jitter.height=0.1),alpha=1.0)+
scale_fill_manual(values=paleta2)+
#geom_violin(aes(fill=colour),alpha=0.3,size=0.4)+scale_color_manual(values=paleta2)+
stat_summary(fun = mean,aes(color = colour,group=colour),fun.min = function(x) mean(x) - sd(x),fun.max = function(x) mean(x)+sd(x),geom = "pointrange", position = position_jitterdodge(jitter.width=0.01), linewidth = 2,size=1.5,shape=22, col="black")+
theme_classic()+labs(x="Light environment",y="Number of prey capture per hour")
###Plot ventral
Create subset
plot_ventral <- tmp_ind %>% filter(Side=="ventral")
plot_ventral$colour<-relevel(as.factor(plot_ventral$colour),"black")Plot the results
#Colours
paleta2<-c("#636363","#f0f0f0","#ffeda0")
#plot
ggplot(data = plot_ventral,aes(x=t_luxes,y=response_nb, fill=colour))+
geom_point(data = plot_ventral,aes(x=t_luxes,y=total_preys, fill=colour),shape = 21,size=4, position = position_jitterdodge(jitter.width=0.15,jitter.height=0.1),alpha=1.0)+
scale_fill_manual(values=paleta2)+
#geom_violin(aes(fill=colour),alpha=0.3,size=0.4)+scale_color_manual(values=paleta2)+
stat_summary(fun = mean,aes(color = colour,group=colour),fun.min = function(x) mean(x) - sd(x),fun.max = function(x) mean(x)+sd(x),geom = "pointrange", position = position_jitterdodge(jitter.width=0.01), linewidth = 2,size=1.5,shape=22, col="black")+
theme_classic()+labs(x="Light environment",y="Number of prey capture per hour")
###Plot by light condition
#Colours
paleta2<-c("#c5dd8f","#5db23a")
#pdf("light_prey.pdf",height=8,width=10)
#plot
ggplot(data = plot_dorsal,aes(x=t_luxes,y=response_nb, fill=t_luxes))+
geom_point(data = plot_dorsal,aes(x=t_luxes,y=total_preys, fill=t_luxes),shape = 21,size=4, position = position_jitterdodge(jitter.width=0.15,jitter.height=0.1),alpha=0.2)+
scale_fill_manual(values=paleta2)+
stat_summary(fun = mean,aes(color = t_luxes,group=t_luxes),fun.min = function(x) mean(x) - sd(x),fun.max = function(x) mean(x)+sd(x),geom = "pointrange", size=1.5, shape=22, col="black")+
theme_classic()+labs(x="Light environment",y="Number of prey capture per hour")
some supplementary plots
tmp_ind %>% filter(total_preys>0) %>% drop_na(Preys_ID) %>% pull(Preys_ID) %>% table() %>% prop.table()## .
## coleoptera diptera hemiptera homoptera hymenoptera lepidoptera
## 0.01775148 0.82840237 0.02958580 0.01183432 0.08875740 0.01775148
## neuroptera
## 0.00591716preys<-tmp_ind %>% filter(total_preys>0) %>% drop_na(Preys_ID) %>% select(colour,Preys_ID,total_preys) %>% uncount(total_preys)
prop.table(table(preys$Preys_ID))*100##
## coleoptera diptera hemiptera homoptera hymenoptera lepidoptera
## 1.5094340 86.4150943 2.2641509 0.7547170 7.1698113 1.5094340
## neuroptera
## 0.3773585paleta2<-c("#f0f0f0","#636363","#ffeda0")
preys %>% ggplot(aes(x=Preys_ID,fill=colour)) + geom_bar(aes(y = after_stat((count/sum(count))*100)),colour="black") + scale_fill_manual(values=paleta2)+theme_classic()+scale_y_continuous(breaks=seq(0,100,10))+labs(y="Percentage of idenfied",x="Prey order")