# Scripts included in the paper # Rinde et al. 2014. The influence of physical factors on kelp and sea urchin distribution # in previously and still grazed areas in the NE Atlantic. Plos One.DOI: 10.1371/journal.pone.0100222 # GAMM analysis of the temperature atlas data # set the path to the working directory where you place the temp data file setwd(".......") library(gdata) # library which support import of excelfiles data_mlongres = read.xls("TempAtlasdata.xlsx") # reads the file with temperature atlas data for 10 m depth within the study area page(data_mlongres,method="print") names(data_mlongres) # [1] "depth" "lat" "lon" "temp" "longresv" "Month" "Year" dim(data_mlongres) # [1] 3325 7 # head(data_mlongres) # depth lat lon temp longresv Month Year #1 10 70.00 20.0 7.221 0.8319050 2 2001 #2 10 65.00 11.0 7.006 0.2220310 2 2007 #3 10 65.00 11.5 6.985 0.5794880 2 2007 #4 10 65.33 11.0 6.980 -0.0683441 2 2007 #5 10 66.00 12.5 6.975 0.3694680 2 1992 #6 10 65.67 11.0 6.967 -0.3679440 2 2007 colMeans(data_mlongres) summary(data_mlongres) detach("package:gdata") library(mgcv) library(MuMIn) # Adds a Site variable tab <- xtabs(~lat + lon, data_mlongres) # creates a frequency table for the observations across latitude and longitude # lon #lat 11 11.5 12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 20 # 65 65 64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # 65.33 72 72 69 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # 65.67 72 72 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # 66 0 72 71 69 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # 66.33 0 0 72 72 69 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # 66.67 0 0 0 72 0 71 0 0 0 0 0 0 0 0 0 0 0 0 0 # 67 0 0 0 0 72 72 71 0 0 0 0 0 0 0 0 0 0 0 0 # 67.33 0 0 0 0 72 72 72 0 0 0 0 0 0 0 0 0 0 0 0 # 67.67 0 0 0 0 0 72 72 71 71 0 0 0 0 0 0 0 0 0 0 # 68 0 0 0 0 0 0 72 71 71 0 0 0 0 0 0 0 0 0 0 # 68.33 0 0 0 0 0 0 0 0 71 71 70 0 0 0 0 0 0 0 0 # 69 0 0 0 0 0 0 0 0 0 0 72 0 72 0 0 0 0 0 0 # 69.33 0 0 0 0 0 0 0 0 0 0 0 72 72 0 16 16 0 0 0 # 69.67 0 0 0 0 0 0 0 0 0 0 0 0 72 72 72 0 16 0 0 # 70 0 0 0 0 0 0 0 0 0 0 0 0 0 72 72 72 0 0 72 # 70.33 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 72 72 72 72 tab <- as.data.frame.table(tab) # converts the table to a dataframe # lat lon Freq #1 65 11 65 #2 65.33 11 72 #3 65.67 11 72 #4 66 11 0 #5 66.33 11 0 #6 66.67 11 0 #7 67 11 0 #8 67.33 11 0 tab <- tab[tab$Freq>0,] # identifies positions with observations page(tab,method="print") dim(tab) # 49 positions / sites tab$lat <- as.numeric(as.character(tab$lat)) tab$lon <- as.numeric(as.character(tab$lon)) data_mlongres$Site <- NA # include a Site variable with NA for all observations for (i in 1:nrow(tab)) { # then include a unik site number from 1-49 to the data sel <- data_mlongres$lat==tab$lat[i] & data_mlongres$lon==tab$lon[i] data_mlongres$Site[sel] <- rep(i, sum(sel)) } # adds a Time variable data_mlongres$Year.fraction <- data_mlongres$Year + (data_mlongres$Month-2)/12 data_mlongres$Time <- 4*(data_mlongres$Year.fraction - 1990) head(data_mlongres) page(data_mlongres,method="print") # Analysis including autocorrelasjon from sampling to sampling within each site temp10m.k3.sampautoc <- gamm(temp ~ s(lat, bs="cr", k=3) + s(longresv, bs="cr", k=3) + s(Year, bs = "cs") + s(Month, bs="cc", k=4), family = gaussian, knots=list(Month=c(0,12)), data=data_mlongres, correlation=corAR1(form=~Time|Site)) # summary of the gam part of the model summary(temp10m.k3.sampautoc$gam) # #Family: gaussian #Link function: identity # #Formula: #temp ~ s(lat, bs = "cr", k = 3) + s(longresv, bs = "cr", k = 3) + # s(Year, bs = "cs") + s(Month, bs = "cc", k = 4) # #Parametric coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) 7.25275 0.01596 454.5 <2e-16 *** #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # #Approximate significance of smooth terms: # edf Ref.df F p-value #s(lat) 1.995 1.995 209.35 < 2e-16 *** #s(longresv) 1.943 1.943 12.50 5.09e-06 *** #s(Year) 8.435 8.435 76.22 < 2e-16 *** #s(Month) 2.000 2.000 11595.25 < 2e-16 *** #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # #R-sq.(adj) = 0.794 Scale est. = 1.6465 n = 3325 # Partial response plots of the mixed GAM model plot(temp10m.k3.sampautoc$gam,all.terms=T, pages=1,scale=0) # scripts to creat the plots in Figure 8 # creates a dataset for predicting latitude versus years for late June (month=6.5) and late July (month=7.5) # for longres values = 0.1. The newgriddata below is for late July, substitute 7.5 by 6.5 in Month, to achieve the same # plot for late June, after performing the predictions below based on the new dataset, and creating a new matrix for plotting. lat.grid <- seq(63, 70.14, by=0.42) # får n antall punkter for en serie ved å ta (max-min)/(n-1) longres.grid <- seq(-1,1.04, by=0.12) year.grid <- seq(1990, 2007, by=1) newgriddata <- data.frame( lat=rep(lat.grid, times=length(year.grid)), Year=rep(year.grid, each=length(lat.grid)), longresv=rep(0.1,18), Month=rep(6.5,18)) names(newgriddata) #"lat" "Year" "longresv" "Month" # predicts the model to the dataset in newgriddata p.newdatatemp10m <- predict.gam(temp10m.k3.sampautoc$gam, newdata=newgriddata, type="response", se.fit = T) # Creates a matrix of these predcitions m.p.newdatatemp10m <- matrix(p.newdatatemp10m$fit, ncol=length(lat.grid), dimnames=list(lat.grid,year.grid)) # plots one of the figures in Figure 8 filled.contour(year.grid, lat.grid, t(m.p.newdatatemp10m) , cex.axis = 2, cex.lab= 1.5, cex.sub=1.5, color = terrain.colors, xlab="Year", ylab="Latitude", plot.axes={ axis(1); axis(2); contour(x = year.grid, y = lat.grid, z = t(m.p.newdatatemp10m), vfont=c("sans serif", "bold"), labcex=2, add=T) })