# script by E. Garcia-Berthou
#
# Supplementary to Magellan & García-Berthou
# Influences of size and sex on invasive species aggression and native species vulnerability: a case for modern regression techniques
# Reviews in Fish Biology and Fisheries
#
# This R script reproduces Fig. 2 top and first row of Table 3
# Warning! ADAPT to your working directory
setwd("C:/Folder_with_data")
FILE<-read.delim("Table S1 R1.txt",header=T, sep="\t")
# Only need to install the packages once
install.packages("VGAM")
install.packages("pscl")
library(pscl)
library(VGAM)
library(bbmle)
# How to cite packages used
citation("VGAM")
citation("pscl")
citation("bbmle")
names(FILE)
summary(FILE)
# Replace the second part of the equalities with names and labels of variables to be analyzed
FILE$y=FILE$Attack.G
LABEL_Y="Gambusia attacks"
FILE$x=FILE$Length.G
LABEL_X="Gambusia length (mm)"
FILE <- FILE[order(FILE$x), ]
FILE$x_centered=FILE$x-mean(FILE$x)
FILE$x_sq=FILE$x_centered*FILE$x_centered
# How to check the (useful) help pages?
# Type e.g. "?vgam" in the Console
?vgam
plot(FILE$y ~ FILE$x,ylab=LABEL_Y,xlab=LABEL_X, pch=19, axes = FALSE)
axis(1);axis(2);box( bty = "l", lwd = 1 )
#####################################
### Quadratic linear regression ###
#####################################
fit_QLR <- lm(log10(y+1) ~ x_centered+x_sq,data=FILE);summary(fit_QLR)
FILE$f <- (10**fitted(fit_QLR))-1;# FILE$f
lines( FILE$x, FILE$f, col="black",lty = 2 )
################################################
### (Conventional GLM with Poisson errors) ###
################################################
fit_GLM <- glm(y ~ x_centered+x_sq,family=poisson,data=FILE);summary(fit_GLM)
# Should be around 1; > 1 => overdispersion
fit_GLM$deviance/fit_GLM$df.residual
# % explained deviance
100*(1-fit_GLM$deviance/fit_GLM$null.deviance)
# (Approximate!) Goodness of fit test
1 - pchisq(summary(fit_GLM)$deviance,summary(fit_GLM)$df.residual)
anova(fit_GLM, test = "Chisq")
FILE$M3cpred <- predict(fit_GLM, se = F, type = "response")
lines(FILE$x, FILE$M3cpred, col="blue", lwd=3,lty = 3)
################################################
### Zero-inflated regression #################
################################################
fit_ZIP <- zeroinfl( y ~ x_centered + x_sq , dist = "poisson",data=FILE);summary(fit_ZIP)
################################################
### Negative-binomial model #################
################################################
fit_NB <- glm.nb(y ~ x_centered+x_sq,data=FILE);summary(fit_NB)
# % explained deviance
100*(1-fit_NB$deviance/fit_NB$null.deviance)
# (Approximate!) Goodness of fit test
1 - pchisq(summary(fit_NB)$deviance, summary(fit_NB)$df.residual)
anova(fit_NB, test = "Chisq")
FILE$fit_NB <- predict(fit_NB, se = F, type = "response")
lines(FILE$x, FILE$fit_NB, col="green", lwd=3,lty=4)
################################################
### NBM with zero-inflation #################
################################################
fit_ZIP_NB <- zeroinfl( y ~ x_centered + x_sq , dist = "negbin",data=FILE);summary(fit_ZIP_NB)
################################################
### Expectile Poisson regression ############
################################################
fit_EPR = vgam(y ~ s(x), fam = amlpoisson(w.aml = c(0.04, 3)),data=FILE)
fit_EPR@extra
summary(fit_EPR)
matlines(FILE$x, fitted(fit_EPR), col = c(6,2), lwd = 2,lty = 1)
################################################
### Comparison of 4 models with AIC #########
################################################
AICtab(fit_GLM,fit_ZIP,fit_NB,fit_ZIP_NB,base=TRUE,logLik=TRUE,weights=TRUE)