Supplement 2.

R script from Tardin et al. Modeling the habitat use of Bryde's whales Balaenoptera edeni off southeastern Brazil



##Getting libraries library(spdep) library(ggplot2) library(car) library(maptools) library(hier.part) library(MuMIn) library(bbmle) library(boot) library(AER)


##########################################################################

#	1 km models	##

##########################################################################

#Getting data

bryde1km <- readShapePoly("Fishnet_1km_RT-KM_Brydeteste3.shp")






#############################################################

#	Data exploration	####

#############################################################

#Visual inspection to include or not a polynomial term depth=bryde1km$DEPTH occbryde=bryde1km$OCC_BRYDE_ plot(depth,occbryde)

#Cubic




distcoast=bryde1km$DISTCOAST

plot(distcoast,occbryde)

#Linear




sstsd=bryde1km$SST_SD_1 plot(sstsd,occbryde)
#Curvilinear




sstmean=bryde1km$SST_MEAN_1 plot(sstmean,occbryde)
#Linear




sstmin=bryde1km$SST_MIN_1 plot(sstmin,occbryde)
#Curvilinear




sstmax=bryde1km$SST_MAX_1 plot(sstmax,occbryde)
#Curvilinear




chlormean=bryde1km$CHLOR_ME_1 plot(chlormean,occbryde)
#Linear




chlormin=bryde1km$CHLOR_MI_1 plot(chlormin,occbryde)
#linear





##########################################################################

##	GLM modelling	##

##########################################################################

##Overdispersion check dispersiontest(glmbase_p) dispersiontest(glmbase_o) dispersiontest(glmbase_f)


# Physiographic model

glmbase_pquasipoly2 <- glm(OCC_BRYDE_~ DEPTH+I(DEPTH^2) + DISTCOAST, data=bryde1km, offset=log(RT_KM),family="quasipoisson",maxit=100)






#Oceanographic model

glmbase_oquasipoly <- glm(OCC_BRYDE_~ SST_SD+I(SST_SD_1^2)+SST_MEAN+SST_MIN_1+I(SST_MIN_1^2)+SST_MAX_1+I(SS T_MAX_1^2)+CHLOR_MEAN+CHLOR_MIN, data=bryde1km, offset=log(RT_KM),family="quasipoisson",maxit=100)



#Full model

glmbase_fquasipoly <- glm(OCC_BRYDE_~ DEPTH+I(DEPTH^2) + DISTCOAST+SST_SD_1+I(SST_SD_1^2)+SST_MEAN+SST_MIN_1+I(SST_MIN_1^2)+SS T_MAX_1+I(SST_MAX_1^2)+CHLOR_MEAN+CHLOR_MIN, data=bryde1km, offset=log(RT_KM),family="quasipoisson",maxit=100)



# GLM Model Selection x.quasipoisson <- function(...) { res <- quasipoisson(...)
res$aic <- poisson(...)$aic

res

}

dfun <- function(object) {

with(object,sum((weights * residuals^2)[weights > 0])/df.residual)

}




qaicglm_pquasipoly2 <- update(glmbase_pquasipoly,family="x.quasipoisson", na.action=na.fail)
qaicglm_oquasipoly <- update(glmbase_oquasipoly,family="x.quasipoisson", na.action=na.fail)
qaicglm_fquasipoly <- update(glmbase_fquasipoly,family="x.quasipoisson", na.action=na.fail)


#qAIC with dispersion argument qAIC(qaicglm_pquasipoly2,dispersion=dfun(qaicglm_pquasipoly))
#560.9




qAIC(qaicglm_oquasipoly,dispersion=dfun(qaicglm_oquasipoly))

#563.9




qAIC(qaicglm_fquasipoly,dispersion=dfun(qaicglm_fquasipoly))

#548.4




#Best GLM summary summary(glmbase_fquasipoly2)



###########################################################################

#	SEV-GLM modelling	##

###########################################################################





##Neighbourhood with Queen specification bryde.nb.q1 <- poly2nb(bryde1km, queen=TRUE) bryde.lw.q1 <- nb2listw(bryde.nb.q1, style="W")


##Moran Test check for spatial autocorrelation moran.test(bryde1km$OCC_BRYDE_,nb2listw(bryde.nb.q1, style="W"))
#Existence of spatial autocorrelation




##Spatial filtering with ME

# Moran Eigenvectors




#Physiographic

errcol2.pquasipoly <- ME(OCC_BRYDE_ ~ DEPTH+(DEPTH^2) + DISTCOAST, data=bryde1km, family="quasipoisson", listw=bryde.lw.q1,offset=log(RT_KM), alpha=0.1, verbose=TRUE)






#Oceanographic

errcol2.oquasipoly <- ME(OCC_BRYDE_ ~ SST_MEAN+SST_SD_1+I(SST_SD_1^2)+SST_MIN_1+I(SST_MIN_1^2)+SST_MAX_1+I(S ST_MAX_1^2)+CHLOR_ME_1+CHLOR_MI_1, data=bryde1km, family="quasipoisson", listw=bryde.lw.q1,offset=log(RT_KM), alpha=0.1, verbose=TRUE)



#Full

errcol2.q_fquasipoly <- ME(OCC_BRYDE_ ~ DEPTH+(DEPTH^2) + DISTCOAST+ SST_MEAN+SST_SD_1+I(SST_SD_1^2)+SST_MIN_1+I(SST_MIN_1^2)+SST_MAX_1+I(S ST_MAX_1^2)+CHLOR_ME_1+CHLOR_MI_1, data=bryde1km, family="quasipoisson", listw=bryde.lw.q1,offset=log(RT_KM), alpha=0.1, verbose=TRUE)



## SEV-GLM


#Physiographic

sevglm_pquasipoly <- glm(OCC_BRYDE_ ~ DEPTH+I(DEPTH^2) + DISTCOAST+fitted(errcol2.pquasipoly), data=bryde1km, family="quasipoisson", offset=log(RT_KM))



#Oceanographic

sevglm_oquasipoly<- glm(OCC_BRYDE_ ~ SST_MEAN+SST_SD_1+I(SST_SD_1^2)+SST_MIN_1+I(SST_MIN_1^2)+SST_MAX_1+I(S ST_MAX_1^2)+CHLOR_ME_1+CHLOR_MI_1+fitted(errcol2.oquasipoly), data=bryde1km, family="quasipoisson", offset=log(RT_KM))



#Full

sevglm_fquasipoly<- glm(OCC_BRYDE_ ~ DEPTH+(DEPTH^2) + DISTCOAST+ SST_MEAN+SST_SD_1+I(SST_SD_1^2)+SST_MIN_1+I(SST_MIN_1^2)+SST_MAX_1+I(S ST_MAX_1^2)+CHLOR_ME_1+CHLOR_MI_1+fitted(errcol2.q_fquasipoly),
data=bryde1km, family="quasipoisson", offset=log(RT_KM))




## Model Selection

#QAIC

x.quasipoisson <- function(...) {

res <- quasipoisson(...) res$aic <- poisson(...)$aic res
}

dfun <- function(object) {

with(object,sum((weights * residuals^2)[weights > 0])/df.residual)

}




qaicsevglm_pquasipoly <- update(sevglm_pquasipoly,family="x.quasipoisson", na.action=na.fail)
qaicsevglm_oquasipoly <- update(sevglm_oquasipoly,family="x.quasipoisson",

na.action=na.fail)

qaicsevglm_fquasipoly <- update(sevglm_fquasipoly,family="x.quasipoisson", na.action=na.fail)
#qAIC with dispersion argument qAIC(qaicsevglm_pquasipoly,dispersion=dfun(qaicsevglm_pquasipoly))
#506.3 qAIC(qaicsevglm_oquasipoly,dispersion=dfun(qaicsevglm_oquasipoly))
#521.0 qAIC(qaicsevglm_fquasipoly,dispersion=dfun(qaicsevglm_fquasipoly))
#513.9






##Summary for best model summary(sevglm_pquasipoly)





##Hierarchical partitioning analysis names(bryde1km) depth1=as.data.frame(depth) distcoast1=as.data.frame(distcoast) env<-c(depth1,distcoast1)
env1<-as.data.frame(env)

hier.part(bryde1km$OCC_BRYDE_, env1)




## test if the fitted values explain the observed value

glmME <- glm(bryde1km$OCC_BRYDE_ ~ fitted(sevglm_pquasipoly), family="quasipoisson")

anova(glmME, test="Chisq")

# Test if SEVGLM is best than GLM

ln.lr <- -2*(logLik(glmbase_pquasipoly)[1]-logLik(sevglm_pquasipoly)[1])

1-pchisq(ln.lr, df=4)




#PSeudoR

psd.r21apoly <- lm(bryde1km$OCC_BRYDE_~ fitted(sevglm_pquasipoly))

summary(psd.r21apoly)$r.square




#############################################################################
######################

##	Spatially predicting B. edeni habitat use	##

#############################################################################
######################




#GLM predbrydeglmquasipoly33=predict.glm(glmbase_pquasipoly,type='response') which(predbrydeglmquasipoly33>2) write.csv(predbrydeglmquasipoly33,file="predbryglmpoly33.csv")



#SEV-GLM predbrydesevglmquasipoly33=predict.glm(sevglm_pquasipoly,type="response") which(predbrydesevglmquasipoly33>1) write.csv(predbrydesevglmquasipoly33,file="brydepredsevglmquasipoly33.csv")


#Residualas mapping

#SEV-GLM residuals

resglm <- residuals.glm(glmbase_pquasipoly, type="pearson")

write.csv(resglm,file="residualsBrydePHYSglm1.csv")

#GLM residuals

ressevglm <- residuals.glm(sevglm_pquasipoly, type="pearson")

write.csv(res,file="residualsBrydePHYS1.csv")




#############################################################################
#####################

##	Eigenvector Mapping	##

#############################################################################
#####################

beta <- matrix(coefficients(sevglm_pquasipoly)[5:9])

x<-as.matrix(fitted(errcol2.pquasipoly)) sf <- x%*%beta write.csv(sf,file="spatialfilterBryde.csv")


##End of coding