Chapter Two Workflow
Fiona Spooner
19 October 2017
Europe
Have ignored the bootstrapping of the models and the use of max, min and precipitation data, and the dredging!
LPI<-read.csv("LPI_pops_20160523_edited.csv")
LPI_pop<-subset(LPI, Specific_location==1 & System !="Marine" & Class != "Actinopterygii"& Class != "Cephalaspidomorphi"& ID != 4438)
ID<-LPI_pop$ID
pop_data<- LPI_pop[,c(1,65:127)]
pop_datab <- (pop_data [,2:64] !="NULL")
points_per_pop1950_2012 = rowSums(pop_datab)
length_id <- data.frame(ID,points_per_pop1950_2012)
LPI_EU<-merge(length_id, LPI_pop, by = "ID")
LPI_EU<-subset(LPI_EU, points_per_pop1950_2012 >=2)
LPI_EU2<-LPI_EU[,c(1:3,66:128)] #just pop trends
LPI_EU2[LPI_EU2 == 'NULL'] = NAlibrary(taRifx)
library(mgcv)
library(zoo)
doFit = function(sp_name) {
spid2 = subset(LPI_EU2, ID == sp_name) #subsetting the population data by each population
spid = spid2[,4:66] #subsetting only the dates
colnames(spid)<-1950:2012 #renaming the date column names as R doesn't like numbered column names
name<-spid2$Binomial
id<-spid2$ID
points<-spid2$points_per_pop1950_2012
name_id<-paste(name, id, sep="_") #creating id for naming files of plots
Date<-as.numeric(colnames(spid))
spidt<-destring(t(spid))
time<-length(min(which(!is.na(spidt))):max(which(!is.na(spidt))))
missing<-time-points
Year<-Date[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
Population<-spidt[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
Population[Population == 0] <- mean(Population, na.rm=TRUE)*0.01
df<-data.frame(Year,Population)
if (points >=6) { ###should I be trying GAMs for populations with less than six points and if that doesn't fit then use a linear model - here I am automatically fitting a linear model if there are less than six points
PopN = log10(Population)
if (length(na.omit(PopN)) >=6) {
SmoothParm = round(length(na.omit(PopN))/2) ####added na.omit in as was getting " A term has fewer unique covariate combinations than specified maximum degrees of freedom" error
} else{
SmoothParm=3
}
mg2<-mgcv:::gam(PopN ~ s(Year, k=SmoothParm), fx=TRUE)
pv2 <- predict(mg2,df,type="response",se=TRUE)
R_sq2<-summary(mg2)$r.sq
model<-1
pv2$fit[pv2$fit <= 0] <- NA
lambda2<-diff(pv2$fit)
lambda_sum2<-sum(lambda2, na.rm=TRUE)
lambda_mean2<-mean(lambda2, na.rm=TRUE)
} else {
SmoothParm<-NA
PopN = log10(Population)
ml2<-lm(PopN~df$Year)
R_sq2<-summary(ml2)$r.sq
model<-0
Pop_interp2<-na.approx(PopN)
Pop_interp2[Pop_interp2<=0] <- NA
lambda2<-diff(Pop_interp2)
lambda_sum2<-sum(lambda2, na.rm=TRUE)
lambda_mean2<-mean(lambda2, na.rm=TRUE)
}
res_df = data.frame(sp_name=sp_name, points=points, SmoothParm=SmoothParm, r_sq=R_sq2, model=model,lambda_sum=lambda_sum2,lambda_mean=lambda_mean2,time=time, missing=missing)
#print(id)
return(res_df)
}
all_df_list <- lapply(unique(LPI_EU2$ID), doFit)all_matrix <- matrix(unlist(all_df_list, use.names =FALSE), ncol=9, byrow=TRUE)
all_df <- data.frame(all_matrix)
colnames(all_df) <- c("ID", "points","SmoothParm", "r_sq", "model", "lambda_sum","lambda_mean", "length_time", "missing_years")
wd<-getwd()write.csv(all_df, paste(wd,"/Chapter_Two/Europe/", "Global_Population_Lambdas.csv", sep=""))Download climate data from here - http://www.ecad.eu/download/ensembles/download.php
rmean<-brick("C:/Users/Fiona/Desktop/PhD/Climate/Europe/MeanT/tg_0.25deg_reg_v11.0.nc", varname = "tg")
tm <- seq(as.Date('1950-01-01'), as.Date('2014-12-31'), 'day')
r2 <- setZ(rmean, tm, 'days')
xmean <- zApply(r2, by=as.yearmon, fun=mean, name='months') #this bit can take ages (hours) - it is summarasing the daily data into monthly data
writeRaster(xmean, filename="Europe_mean_mon.grd", bandorder='BIL', overwrite=TRUE)library(raster)## Warning: package 'raster' was built under R version 3.4.2## Loading required package: sp## Warning: package 'sp' was built under R version 3.4.2##
## Attaching package: 'raster'## The following object is masked from 'package:nlme':
##
## getData## The following object is masked from 'package:taRifx':
##
## shiftxy<-cbind(LPI_EU$Longitude, LPI_EU$Latitude)
id<-LPI_EU$ID
binomial<-LPI_EU$Binomial
xmean<- brick(paste(wd,"/Chapter_Two/Europe/", "Europe_mean_mon.grd", sep=""))all_EU_mean_data<-data.frame(id=numeric(0), Binomial=character(0),xy=numeric(0), year=numeric(0), month=numeric(0), rasterex=numeric(0))
for (i in 1:nlayers(xmean)) {
rasterex<- extract(xmean[[i]], xy, buffer=50000, fun=mean, na.rm=TRUE)
date<-names(xmean[[i]])
date2 <- as.yearmon(date, "%b.%Y")
dataex<-data.frame(id, binomial,xy, format(date2, "%Y"),format(date2, "%b"), rasterex)
# print(date2)
all_EU_mean_data = rbind(all_EU_mean_data, dataex)
}colnames(all_EU_mean_data)[c(2:7)]<-c("Binomial","Longitude","Latitude", "Year", "Month", "Mean_Temp")
all_EU_mean_data_na_omit<-na.omit(all_EU_mean_data)
write.csv(all_EU_mean_data, paste(wd,"/Chapter_Two/Europe/" ,"Mean_Temp_E-OBS_50k_buff.csv", sep=""))
write.csv(all_EU_mean_data_na_omit, paste(wd,"/Chapter_Two/Europe/" ,"Mean_Temp_E-OBS_50k_buff_na_omitted.csv", sep=""))all_EU_mean_data_na_omit<-read.csv(paste(wd,"/Chapter_Two/Europe/" ,"Mean_Temp_E-OBS_50k_buff_na_omitted.csv", sep=""))library(plyr)## Warning: package 'plyr' was built under R version 3.4.2library(broom)
LPI_EU2<-LPI_EU[LPI_EU$ID %in% all_EU_mean_data_na_omit$id & LPI_EU$ID!=18236 & LPI_EU$ID!=18239 ,]
nrow(LPI_EU2)## [1] 1118doMean = function(sp_name) {
spid2 = subset(LPI_EU2, ID == sp_name) #subsetting the population data by each population
spid = spid2[,66:128] #subsetting only the dates
colnames(spid)<-1950:2012 #renaming the date column names as R doesn't like numbered column names
climid=subset(all_EU_mean_data_na_omit, id == sp_name) #subsetting the climate data by each population
year_temp <- ddply(climid, "Year", summarise, #calculating the annual mean for max temp, min temp and precipitation
mean_mean = mean(na.omit(Mean_Temp)))
lt_avg <- ddply(climid, "id", summarise, #calculating the long term mean (1950-2014) of the max temp, min temp and precipitation - for each population
lt_avg_mean = mean(na.omit(Mean_Temp)))
year_temp$anom_mean = year_temp$mean_mean - lt_avg$lt_avg_mean #calculating anomalies 1950-2014
name<-spid2$Binomial
id<-spid2$ID
points<-spid2$points_per_pop1950_2012
name_id<-paste(name, id, sep="_") #creating id for naming files of plots
Date<-as.numeric(colnames(spid))
spidt<-destring(t(spid))
Year<-Date[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
Population<-spidt[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
Mean_mon<-climid[climid$Year %in% Year, ]$Mean_Temp
Mean_anom<-year_temp$anom_mean[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
Mean<-year_temp$mean_mean[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
if (sum(is.nan(Mean))!=length(Mean)){
lm_mean<-lm(Mean~Year)
lm_mean_df<-tidy(lm_mean)[2,]
mean_df<-cbind(id,lm_mean_df)
} else{
mean_df<-matrix(c(id,NA,NA,NA,NA,NA), nrow=1, ncol=6)
colnames(mean_df)<-c("id", "term", "estimate", "std.error", "statistic", "p.value")
mean_df<-data.frame(mean_df)
}
#print(id)
return(mean_df)
}
all_df_list <- lapply(LPI_EU2$ID, doMean)all_matrix <- matrix(unlist(all_df_list), ncol=6, byrow=TRUE)
mean_df <- data.frame(all_matrix)
colnames(mean_df) <- c("ID", "Term","Estimate","SE","Statistic","p.val")
write.csv(mean_df, paste(wd,"/Chapter_Two/Europe/" ,"Rate_Mean_Temp_Change.csv", sep=""))Land use - HILDA dataset from http://www.wageningenur.nl/en/Expertise-Services/Chair-groups/Environmental-Sciences/Laboratory-of-Geoinformation-Science-and-Remote-Sensing/Models/Hilda.htm
library(maptools)## Warning: package 'maptools' was built under R version 3.4.2## Checking rgeos availability: TRUEr_1950 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch1950.asc", sep="")))
r_1960 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch1960.asc", sep="")))
r_1970 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch1970.asc", sep="")))
r_1980 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch1980.asc", sep="")))
r_1990 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch1990.asc", sep="")))
r_2000 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch2000.asc", sep="")))
r_2010 <-raster(readAsciiGrid( paste(wd,"/Chapter_Two/Europe/" ,"eu27ch2010.asc", sep="")))
r_2010b<-crop(r_2010, r_1950) #the 2010 raster is larger than the other years so has to be ropped to the same size so that they can be stacked
r<-stack(r_1950,r_1960,r_1970,r_1980,r_1990,r_2000,r_2010b)
#writeRaster(r, paste(wd,"/Chapter_Two/Europe/" ,"Land_Use_1950_2010.grd", sep=""))
names(r)<-seq(1950,2010, by=10)
crs.geo <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
proj4string(r) <- crs.geolibrary(rgdal)## Warning: package 'rgdal' was built under R version 3.4.2## rgdal: version: 1.2-13, (SVN revision 686)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
## Path to GDAL shared files: C:/Users/Fiona/R/win-library/3.4/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/Fiona/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-5xy<-cbind(LPI_EU$Longitude, LPI_EU$Latitude)
id<-LPI_EU$ID
xy_eu<-project(xy,"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
check_xy<-data.frame(id,xy_eu) #subsetting the populations so that only those within the landuse raster are included - means that non-EU27 populations are lost
check_xy$check<-extract(r,check_xy[,(2:3)])
check_xy2<-na.omit(check_xy)
#check_xy2<-subset(check_xy, check != "NA")
crop_id<-check_xy2[,1]
selectedRows <- (LPI_EU$ID %in% check_xy2$id)
LPI_EU_sel <- LPI_EU[selectedRows,]
points <- SpatialPoints(cbind(check_xy2[,2], check_xy2[,3]))
x<-check_xy2[,2]
y<-check_xy2[,3]
xy<-cbind(x,y) #creating the grid around the population - 5km either side gives a grid of 121km^2
xmin<- colFromX(r[[1]], points[1:length(points),]) - 5
xmax<- colFromX(r[[1]], points[1:length(points),]) + 5
ymin<- rowFromY(r[[1]], points[1:length(points),]) - 5
ymax<- rowFromY(r[[1]], points[1:length(points),]) + 5 #changed from 5 to 12
grid_crop<-data.frame(ID=crop_id,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) #setting extents of grid to extract
grid_crop2<-merge(grid_crop,LPI_EU_sel, by="ID")
grid_crop2<-subset(grid_crop2, Specific_location ==1) result <- data.frame() #empty result dataframe
for (i in 1:length(grid_crop2$ID)) {
ID<-grid_crop2[i,1]
Binomial<-as.character(grid_crop2[i,7])
spid = grid_crop2[i,70:132] #subsetting only the dates
colnames(spid)<-1950:2012
Date<-as.numeric(colnames(spid))
spidt<-destring(t(spid))
Year<-Date[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
Population<-spidt[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
LUmin<-min(round_any(Year, 10))
LUmax<-max(round_any(Year, 10))
crop_check<-crop(r, extent(r, grid_crop2[i,4],grid_crop2[i,5],grid_crop2[i,2],grid_crop2[i,3]))
crop_df<-data.frame(as.matrix(crop_check))
#crop_df<-na.omit(crop_df)
colnames(crop_df)<-seq(1950,2010, by=10)
col<-as.numeric(colnames(crop_df[,1:7]))
min_yr<-min(which(col>=LUmin))
max_yr<-max(which(col<=LUmax))
if (min_yr != max_yr) {
crop_df$check<-rowSums(crop_df[,min_yr] == crop_df[,c(min_yr:max_yr)])
}else{
crop_df$check<-1
}
crop_df$change<- crop_df$check != length(seq(min_yr,max_yr, by=1))
change_rate<-sum(na.omit(crop_df$change=="TRUE"))/length(na.omit(crop_df$change))
#print(change_rate)
final<-cbind(ID,Binomial,change_rate)
result<-rbind(final,result)
}write.csv(result, paste(wd,"/Chapter_Two/Europe/" ,"LPI_LandUse_121_18_09_2015.csv", sep="" ))
LPI_LUC<-merge(LPI_EU_sel, result[,c(1,3)], by="ID",all = TRUE)bm<-read.csv("Amniote_Database_Aug_2015.csv")
bm$Binomial<-paste(bm$genus, bm$species, sep="_")
lpibm<-merge(LPI, bm, by="Binomial")
lpibm<-data.frame(lpibm$Binomial, lpibm$ID, lpibm$adult_body_mass_g)
head(lpibm)
colnames(lpibm)<-c("Binomial", "ID", "Body_mass_g")
lpibm$Log_Body_Mass_g<-log10(lpibm$Body_mass_g)
lpibm2<-lpibm[lpibm$Body_mass_g !=-999,]
lpibm2<-unique(lpibm2)
lpibm2$Log_Body_Mass_g<-log10(lpibm2$Body_mass_g)
write.csv(lpibm2, paste(wd, "/Chapter_Two/Europe/" ,"LPI_BodyMass_Amniote_Database.csv", sep=""))lpi<-read.csv("LPI_pops_20160523_edited.csv")
lpi<-lpi[,c(1,15,32,31)]
pop<-read.csv(paste(wd,"/Chapter_Two/Europe/", "Global_Population_Lambdas.csv", sep=""))
pop<-pop[,c(2,5,8,9)]
clim<-read.csv(paste(wd,"/Chapter_Two/Europe/" ,"Rate_Mean_Temp_Change.csv", sep=""))
clim<-clim[,c(2,4)]
luc<- read.csv(paste(wd,"/Chapter_Two/Europe/" ,"LPI_LandUse_121_18_09_2015.csv", sep="" ))
luc<-luc[,c(2,4)]
bm<- read.csv(paste(wd, "/Chapter_Two/Europe/", "LPI_BodyMass_Amniote_Database.csv", sep=""))
bm<-bm[,-1]
df<-merge(merge(pop, luc,by="ID", all=T), merge(clim, bm, by="ID", all=T), by="ID", all=T)
df<-merge(df, lpi, by="ID")
head(df)## ID r_sq lambda_mean length_time change_rate Estimate
## 1 1 NA NA NA NA NA
## 2 2 0.9311001 0.03987404 30 NA NA
## 3 3 NA NA NA NA NA
## 4 4 0.9030387 -0.02074499 15 NA NA
## 5 5 NA NA NA NA NA
## 6 6 0.8831535 -0.04365634 21 NA NA
## Binomial Body_mass_g Log_Body_Mass_g Class Longitude
## 1 Balearica_regulorum 3614.5 3.558048 Aves 31.13306
## 2 Acrocephalus_sechellensis 15.9 1.201397 Aves 55.66667
## 3 Copsychus_sechellarum 71.0 1.851258 Aves 55.66667
## 4 Copsychus_sechellarum 71.0 1.851258 Aves 55.93333
## 5 Falco_punctatus 160.0 2.204120 Aves 57.58333
## 6 Pternistis_ochropectus 809.0 2.907949 Aves 42.65806
## Latitude
## 1 -0.07889
## 2 -4.33333
## 3 -4.58333
## 4 -4.58333
## 5 -20.30000
## 6 11.76667Mammals
df_mammal<-subset(df, !is.na(lambda_mean)& !is.na(change_rate) & !is.na(Estimate) & !is.na(Binomial) & length_time > 5 & Class=="Mammalia")
nrow(df_mammal)## [1] 83library(data.table)
library(plyr)
sp_dups<-data.frame(ddply(df_mammal,.(Longitude,Latitude),nrow))
sp_dups$loc_id<-1:length(sp_dups$Longitude)
df_mammal<-merge(sp_dups, df_mammal, by=c("Longitude","Latitude"))
parm_df<-df_mammal[,c("ID","change_rate", "Estimate", "Log_Body_Mass_g")] ##ID, land use, and climate
parm_mat<-as.matrix(parm_df)
parm_scale_m<-scale(parm_mat[,c("change_rate", "Estimate", "Log_Body_Mass_g")]) #use the scaling factors at the bottom of these to scale the rasters
parm_id<-parm_mat[,"ID"]
parm_df_scale<-data.frame(parm_id,parm_scale_m)
colnames(parm_df_scale)<-c("ID","change_rate_scale","mean_slope_scale", "body_mass_scale" )
sp_df_scale<-merge(df_mammal, parm_df_scale, by="ID")
dt_mammal<-data.table(sp_df_scale)source("rsquaredglmm.R")
library(lme4)
library(MuMIn)
mm0<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+change_rate_scale:mean_slope_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm0a<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm0b<-lmer(lambda_mean ~ change_rate_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm0c<-lmer(lambda_mean ~ mean_slope_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm0d<-lmer(lambda_mean ~ body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm1<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+change_rate_scale:mean_slope_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm1a<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm1b<-lmer(lambda_mean ~ change_rate_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mm1c<-lmer(lambda_mean ~ mean_slope_scale+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mmnull<-lmer(lambda_mean ~ 1+(1|Binomial)+(1|loc_id),data=dt_mammal, REML=F)
mmsAICc <- model.sel(mm0,mm0a,mm0b,mm0c,mm0d,mm1,mm1a,mm1b,mm1c,mmnull)
mmsAICc$model<-rownames(mmsAICc)
mmsAICc<-data.frame(mmsAICc)
mmsAICc## X.Intercept. body_mass_scale change_rate_scale mean_slope_scale
## mm0b 0.01524545 0.01557146 -0.01877577 NA
## mm1b 0.01167096 NA -0.01877543 NA
## mm0a 0.01555543 0.01369557 -0.01891198 0.0031686026
## mm1a 0.01266659 NA -0.01894963 0.0040690305
## mm1 0.01257949 NA -0.02323793 0.0090951017
## mm0 0.01543953 0.01323369 -0.01968707 0.0040951884
## mm0d 0.01904809 0.01115632 NA NA
## mmnull 0.01661278 NA NA NA
## mm0c 0.01896941 0.01266544 NA -0.0034826605
## mm1c 0.01666464 NA NA 0.0005387514
## change_rate_scale.mean_slope_scale df logLik AICc delta
## mm0b NA 6 148.3299 -283.5546 0.000000
## mm1b NA 5 146.5883 -282.3974 1.157156
## mm0a NA 7 148.7748 -282.0562 1.498340
## mm1a NA 6 147.4821 -281.8590 1.695599
## mm1 0.009018657 7 147.7699 -280.0464 3.508128
## mm0 0.001628487 8 148.7828 -279.6196 3.935003
## mm0d NA 5 138.9615 -267.1437 16.410852
## mmnull NA 4 137.4392 -266.3655 17.189077
## mm0c NA 6 139.1490 -265.1928 18.361767
## mm1c NA 5 137.4441 -264.1090 19.445594
## weight model
## mm0b 3.603206e-01 mm0b
## mm1b 2.020300e-01 mm1b
## mm0a 1.703447e-01 mm0a
## mm1a 1.543456e-01 mm1a
## mm1 6.236039e-02 mm1
## mm0 5.037488e-02 mm0
## mm0d 9.842779e-05 mm0d
## mmnull 6.670040e-05 mmnull
## mm0c 3.710924e-05 mm0c
## mm1c 2.158400e-05 mm1c AIC(mm0,mm0a,mm0b,mm0c,mm0d,mm1,mm1a,mm1b,mm1c,mmnull)## df AIC
## mm0 8 -281.5655
## mm0a 7 -283.5496
## mm0b 6 -284.6598
## mm0c 6 -266.2981
## mm0d 5 -267.9229
## mm1 7 -281.5398
## mm1a 6 -282.9642
## mm1b 5 -283.1766
## mm1c 5 -264.8882
## mmnull 4 -266.8783 mmodels_list<-list(mm0,mm0a,mm0b,mm0c,mm0d,mm1,mm1a,mm1b,mm1c,mmnull)
mmodelsR<-lapply(mmodels_list,rsquared.glmm)
mmodelsRsq <- matrix(unlist(mmodelsR), ncol=6, byrow=T)
rownames(mmodelsRsq)<-c("m0","m0a","m0b","m0c","m0d","m1","m1a","m1b","m1c","mnull")
mmodelsRsq## [,1] [,2] [,3] [,4] [,5] [,6]
## m0 1 1 1 0.1564913380 0.9999483 -281.5655
## m0a 1 1 1 0.1570750292 0.9999483 -283.5496
## m0b 1 1 1 0.1607273449 0.9999165 -284.6598
## m0c 1 1 1 0.0574337584 0.2543548 -266.2981
## m0d 1 1 1 0.0533727107 0.2636722 -267.9229
## m1 1 1 1 0.1113336396 0.9999512 -281.5398
## m1a 1 1 1 0.1036501787 0.9999498 -282.9642
## m1b 1 1 1 0.1016993754 0.9999155 -283.1766
## m1c 1 1 1 0.0001235974 0.3063563 -264.8882
## mnull 1 1 1 0.0000000000 0.2990910 -266.8783 library(MuMIn)
mvar_imp<-summary(model.avg(mmodels_list))
#mav<-model.avg(models_list, subset = cumsum(weight) <= 0.95)
mmav<-model.avg(mmodels_list, subset = delta <= 6)
msmav<-summary(mmav)
mcoef_av<-msmav$coefmat.subset[1:5,"Estimate"]
mcoef_df<-data.frame(mcoef_av)
mcoef_df$lowCI<-confint(mmav)[1:5,1]
mcoef_df$highCI<-confint(mmav)[1:5,2]
mcoef_df## mcoef_av lowCI highCI
## (Intercept) 0.014022223 -0.004759333 0.03280378
## change_rate_scale -0.019149988 -0.023615432 -0.01468454
## body_mass_scale 0.014818061 -0.001952104 0.03158823
## mean_slope_scale 0.004437057 -0.005221473 0.01409559
## change_rate_scale:mean_slope_scale 0.005716417 -0.019824273 0.03125711 mcoef_pcnt<-data.frame(((10^mcoef_df) - 1)*100)
mcoef_pcnt## mcoef_av lowCI highCI
## (Intercept) 3.281425 -1.0898940 7.845935
## change_rate_scale -4.313645 -5.2924575 -3.324716
## body_mass_scale 3.470861 -0.4484799 7.544505
## mean_slope_scale 1.026907 -1.1950900 3.298874
## change_rate_scale:mean_slope_scale 1.324954 -4.4620920 7.462541Birds
df_birds<-subset(df, !is.na(lambda_mean)& !is.na(change_rate) & !is.na(Estimate) & !is.na(Binomial) & length_time > 5 & Class=="Aves")
nrow(df_birds)## [1] 409library(data.table)
sp_dups<-data.frame(ddply(df_birds,.(Longitude,Latitude),nrow))
sp_dups$loc_id<-1:length(sp_dups$Longitude)
df_birds<-merge(sp_dups, df_birds, by=c("Longitude","Latitude"))
parm_df<-df_birds[,c("ID","change_rate", "Estimate", "Log_Body_Mass_g")] ##ID, land use, and climate
parm_mat<-as.matrix(parm_df)
parm_scale_b<-scale(parm_mat[,c("change_rate", "Estimate", "Log_Body_Mass_g")]) #use the scaling factors at the bottom of these to scale the rasters
parm_id<-parm_mat[,"ID"]
parm_df_scale<-data.frame(parm_id,parm_scale_b)
colnames(parm_df_scale)<-c("ID","change_rate_scale","mean_slope_scale", "body_mass_scale" )
sp_df_scale<-merge(df_birds, parm_df_scale, by="ID")
dt_birds<-data.table(sp_df_scale)source("rsquaredglmm.R")
library(lme4)
library(MuMIn)
bm0<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+change_rate_scale:mean_slope_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm0a<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm0b<-lmer(lambda_mean ~ change_rate_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm0c<-lmer(lambda_mean ~ mean_slope_scale+body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm0d<-lmer(lambda_mean ~ body_mass_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm1<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+change_rate_scale:mean_slope_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm1a<-lmer(lambda_mean ~ change_rate_scale+mean_slope_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm1b<-lmer(lambda_mean ~ change_rate_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bm1c<-lmer(lambda_mean ~ mean_slope_scale+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bmnull<-lmer(lambda_mean ~ 1+(1|Binomial)+(1|loc_id),data=dt_birds, REML=F)
bmsAICc <- model.sel(bm0,bm0a,bm0b,bm0c,bm0d,bm1,bm1a,bm1b,bm1c,bmnull)
bmsAICc$model<-rownames(bmsAICc)
bmsAICc<-data.frame(bmsAICc)
bmsAICc## X.Intercept. body_mass_scale change_rate_scale mean_slope_scale
## bm0d 0.0012163996 0.010485119 NA NA
## bm0c 0.0020197094 0.010919352 NA -0.003405315
## bm0b 0.0008508011 0.009955871 -0.001792717 NA
## bm0a 0.0016608104 0.010356137 -0.001907385 -0.003447511
## bmnull 0.0043160134 NA NA NA
## bm1b 0.0024175506 NA -0.006537363 NA
## bm0 0.0041861363 0.009393373 0.002609988 0.001987205
## bm1c 0.0048013764 NA NA -0.001974737
## bm1a 0.0030156109 NA -0.006718286 -0.002388365
## bm1 0.0072265945 NA 0.001977744 0.006921863
## change_rate_scale.mean_slope_scale df logLik AICc delta
## bm0d NA 5 527.5369 -1044.925 0.000000
## bm0c NA 6 527.8588 -1043.509 1.416298
## bm0b NA 6 527.5854 -1042.962 1.963125
## bm0a NA 7 527.9142 -1041.549 3.375741
## bmnull NA 4 524.3593 -1040.620 4.305310
## bm1b NA 5 525.1149 -1040.081 4.843926
## bm0 0.009171693 8 528.1878 -1040.016 4.909370
## bm1c NA 5 524.4700 -1038.791 6.133870
## bm1a NA 6 525.2735 -1038.338 6.586919
## bm1 0.016008974 7 526.1804 -1038.081 6.843484
## weight model
## bm0d 0.40661149 bm0d
## bm0c 0.20027851 bm0c
## bm0b 0.15236754 bm0b
## bm0a 0.07518771 bm0a
## bmnull 0.04723820 bmnull
## bm1b 0.03608564 bm1b
## bm0 0.03492396 bm0
## bm1c 0.01893331 bm1c
## bm1a 0.01509553 bm1a
## bm1 0.01327810 bm1 AIC(bm0,bm0a,bm0b,bm0c,bm0d,bm1,bm1a,bm1b,bm1c,bmnull)## df AIC
## bm0 8 -1040.376
## bm0a 7 -1041.828
## bm0b 6 -1043.171
## bm0c 6 -1043.718
## bm0d 5 -1045.074
## bm1 7 -1038.361
## bm1a 6 -1038.547
## bm1b 5 -1040.230
## bm1c 5 -1038.940
## bmnull 4 -1040.719 bmodels_list<-list(bm0,bm0a,bm0b,bm0c,bm0d,bm1,bm1a,bm1b,bm1c,bmnull)
bmodelsR<-lapply(bmodels_list,rsquared.glmm)
bmodelsRsq <- matrix(unlist(bmodelsR), ncol=6, byrow=T)
rownames(bmodelsRsq)<-c("m0","m0a","m0b","m0c","m0d","m1","m1a","m1b","m1c","mnull")
bmodelsRsq## [,1] [,2] [,3] [,4] [,5] [,6]
## m0 1 1 1 0.0266402581 0.2894048 -1040.376
## m0a 1 1 1 0.0238195371 0.2874339 -1041.828
## m0b 1 1 1 0.0235873521 0.2920669 -1043.171
## m0c 1 1 1 0.0227498644 0.2875633 -1043.718
## m0d 1 1 1 0.0221258867 0.2914887 -1045.074
## m1 1 1 1 0.0164147169 0.2729675 -1038.361
## m1a 1 1 1 0.0084387334 0.2663729 -1038.547
## m1b 1 1 1 0.0087079173 0.2707932 -1040.230
## m1c 1 1 1 0.0008019437 0.2609907 -1038.940
## mnull 1 1 1 0.0000000000 0.2631177 -1040.719 library(MuMIn)
bvar_imp<-summary(model.avg(bmodels_list))
#mav<-model.avg(models_list, subset = cumsum(weight) <= 0.95)
bmav<-model.avg(bmodels_list, subset = delta <= 6)
bsmav<-summary(bmav)
nums<-c("1", "2", "3", "12", "13", "23", "123", "234", "1234", "(Null)")
mods<-c("bm0d", "bm1b", "bm1c", "bm0b", "bm0c", "bm1a", "bm0a", "bm1", "bm0", "bmnull")
mod_num<-cbind(nums, mods)
b_models<-rownames(bvar_imp$msTable)
b_weights<-Weights(bmav)
b_mod_weights<-cbind(b_models, b_weights)
b_model_weights<-merge(b_mod_weights, mod_num, by.x="b_models", by.y="nums")
bcoef_av<-bsmav$coefmat.subset[1:5,"Estimate"]
bcoef_df<-data.frame(bcoef_av)
bcoef_df$lowCI<-confint(bmav)[1:5,1]
bcoef_df$highCI<-confint(bmav)[1:5,2]
bcoef_df## bcoef_av lowCI highCI
## (Intercept) 0.001823313 -0.01000010 0.013646727
## body_mass_scale 0.010437385 0.00213598 0.018738790
## mean_slope_scale -0.002385688 -0.01313264 0.008361269
## change_rate_scale -0.001946767 -0.01495797 0.011064440
## change_rate_scale:mean_slope_scale 0.011055142 -0.01371498 0.035825260 bcoef_pcnt<-data.frame(((10^bcoef_df) - 1)*100)
bcoef_pcnt## bcoef_av lowCI highCI
## (Intercept) 0.4207159 -2.276300 3.192166
## body_mass_scale 2.4324086 0.493039 4.409205
## mean_slope_scale -0.5478189 -2.978641 1.943906
## change_rate_scale -0.4472564 -3.385563 2.580412
## change_rate_scale:mean_slope_scale 2.5782160 -3.108646 8.598858bcoef_pcnt$Class<-"Birds"
mcoef_pcnt$Class<-"Mammals"
bcoef_pcnt$Var<-rownames(bcoef_pcnt)
mcoef_pcnt$Var<-rownames(mcoef_pcnt)
colnames(bcoef_pcnt)[1]<-"ave"
colnames(mcoef_pcnt)[1]<-"ave"
coef_pcnt<-rbind(bcoef_pcnt, mcoef_pcnt)
library(ggplot2)
p1<-ggplot(coef_pcnt, aes(colour=Class))
p1<- p1 + geom_hline(yintercept = 0, colour=gray(1/2), lty=2)
p1<- p1 + geom_linerange(aes(x=Var, ymin=lowCI, ymax=highCI), lwd=2.5, position = position_dodge(width=2/3))
p1<- p1 + geom_pointrange(aes(x= Var, y=ave, ymin=lowCI, ymax=highCI), lwd=2, position=position_dodge(width=2/3), shape=21, fill="White")
#p1<- p1 + scale_y_continuous(breaks=seq(-8, 6, 2), limits=(c(-9,5))) +theme_bw() + labs(y = "Population Change (%)", x = "") + theme(legend.position="none",text=element_text(size=20),axis.text.x=element_text(size=8) , axis.title.x = element_text(margin = unit(c(5, 0, 0, 0), "mm")))
#p1<- p1 + scale_color_manual(values=c("black", "black"))
p1<-p1+coord_flip()
print(p1)
Spatial Predictions - Europe
mean_t<- brick(paste(wd,"/Chapter_Two/Europe/" ,"Europe_mean_mon.grd", sep=""))
mean_t2<-trim(mean_t)
x<-seq(-38.625, 75.125, by=0.25)
y<-seq(25.375, 71.875, by=0.25)
xy<-expand.grid(x,y)
colnames(xy)<-c("Lon", "Lat")
time<-1:nlayers(mean_t)
## add 1 for a model with an intercept
X <- cbind(1, time)
## pre-computing constant part of least squares
invXtX <- solve(t(X) %*% X) %*% t(X)
## much reduced regression model; [2] is to get the slope
quickfun <- function(y) (invXtX %*% y)[2]
x4 <- calc(mean_t, quickfun)writeRaster(x4, filename=paste(wd,"/Chapter_Two/Europe/" ,"climate_change_slope.tif", sep=""), format="GTiff", overwrite=TRUE)r<-brick(paste(wd,"/Chapter_Two/Europe/" ,"Land_Use_1950_2010.grd", sep=""))
crs.geo <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
proj4string(r) <- crs.geo
names(r)<-seq(1950,2010, by=10)
plot(r)
r_df<-data.frame(as.matrix(r))
r_df$a<-as.numeric(r_df[,1] != r_df[,2]) #checking for each decade if there is a change in land use category
r_df$b<-as.numeric(r_df[,2] != r_df[,3])
r_df$c<-as.numeric(r_df[,3] != r_df[,4])
r_df$d<-as.numeric(r_df[,4] != r_df[,5])
r_df$e<-as.numeric(r_df[,5] != r_df[,6])
r_df$f<-as.numeric(r_df[,6] != r_df[,7])
r_df$total_change<-rowSums(r_df[,c(8:13)]) #total number of category changes - a way of quantifying how unstable a landscape is?
r_rast<-raster(matrix(r_df$total_change, nrow=4075,ncol=4011, byrow=T))
plot(r_rast)
extent(r_rast)<-c(2554419, 6565419, 1260409, 5335409)
crs(r_rast)<-"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
r_rast
n <- c(0,0.99,0,1,6,1)
m <- matrix(n, ncol=3, byrow=TRUE)
rc <- reclassify(r_rast, m) #reclassifying to binary rather than number of changes
f<-rep.int(1, 121)
m<-matrix(f, nrow=11)
my_fun2 <- function(x) {
result <- sum(na.omit(x))/length(na.omit(x))
return(result)
}
rm <- focal(rc, w=matrix(1,nrow=11,ncol=11), fun=my_fun2) #looking at amount of change in the surrounding area
plot(rm)writeRaster(rm, paste(wd,"/Chapter_Two/Europe/" ,"land_use_change.tif", sep=""), format="GTiff", overwrite=TRUE )Spatial Prediction - Europe Birds
library(raster)
LU<-raster( paste(wd,"/Chapter_Two/Europe/" ,"land_use_change.tif", sep=""))
TM<-raster(paste(wd,"/Chapter_Two/Europe/" ,"climate_change_slope.tif", sep=""))
crs_laea<-"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
extent(LU)<-c(2554419, 6565419, 1260409, 5335409)
crs(LU)<-crs_laea
TM<-projectRaster(TM, crs = crs_laea)
#plot(TM)
centre_temp_b<-attr(parm_scale_b, 'scaled:center')[2]
centre_luc_b<-attr(parm_scale_b, 'scaled:center')[1]
scale_temp_b<-attr(parm_scale_b, 'scaled:scale')[2]
scale_luc_b<-attr(parm_scale_b, 'scaled:scale')[1]
LU_b<-(LU - centre_luc_b )/scale_luc_b
TM_b<-(TM -centre_temp_b)/ scale_temp_b
LU_bt<-trim(LU_b)
TM_b2<-resample(TM_b, LU_bt)
#LU_bt2<-resample(LU_ust,TM_us)
#bm<-reclassify(TM_b, c(-1,max(na.omit(values(TM_b))),0))
bird_stack<-stack(TM_b2,LU_bt)
names(bird_stack)<-c("mean_slope_scale","change_rate_scale")
# bird_average<-function(model, weight){
# pred<-predict(bird_stack,model, re.form=NA)
# pred_weight<-pred*weight
# return(pred_weight)
# }
#
# b_mod_weights$b_weights<-as.numeric(as.character(b_model_weights$b_weights))
#
# mapply(bird_average, model=bmodels_list, weight=b_model_weights$b_weights)
pred<-predict(bird_stack,bm1, re.form=NA)
plot(trim(10^pred))
xy <- data.frame(lon=dt_birds$Longitude, lat=dt_birds$Latitude)
coordinates(xy) <- c("lon", "lat")
proj4string(xy) <- CRS("+init=epsg:4326") # WGS 84 UTM 35S
CRS.new <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
xy_new<- spTransform(xy, CRS.new)
points(xy_new)
Spatial Prediction - Europe Mammals
library(raster)
LU<-raster( paste(wd,"/Chapter_Two/Europe/" ,"land_use_change.tif", sep=""))
TM<-raster(paste(wd,"/Chapter_Two/Europe/" ,"climate_change_slope.tif", sep=""))
crs_laea<-"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
extent(LU)<-c(2554419, 6565419, 1260409, 5335409)
crs(LU)<-crs_laea
TM<-projectRaster(TM, crs = crs_laea)
#plot(TM)
centre_temp_m<-attr(parm_scale_m, 'scaled:center')[2]
centre_luc_m<-attr(parm_scale_m, 'scaled:center')[1]
scale_temp_m<-attr(parm_scale_m, 'scaled:scale')[2]
scale_luc_m<-attr(parm_scale_m, 'scaled:scale')[1]
LU_m<-(LU - centre_luc_m )/scale_luc_m
TM_m<-(TM -centre_temp_m)/ scale_temp_m
LU_mt<-trim(LU_m)
TM_m2<-resample(TM_m, LU_mt)
mammal_stack<-stack(TM_m2,LU_mt)
names(mammal_stack)<-c("mean_slope_scale","change_rate_scale")
pred<-predict(mammal_stack,mm1, re.form=NA)
plot(trim(10^pred))
xy <- data.frame(lon=dt_mammal$Longitude, lat=dt_mammal$Latitude)
coordinates(xy) <- c("lon", "lat")
proj4string(xy) <- CRS("+init=epsg:4326") # WGS 84 UTM 35S
CRS.new <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
xy_new<- spTransform(xy, CRS.new)
points(xy_new)