Dune Vegetation
Data Exploration
and
Spatial-temporal multi-species joint model for vegetation dynamics
Data Source: Miller Lab: http://www.bio.fsu.edu/~miller/StGeorge/index.htm
Code Contact: John: jmh09r@my.fsu.edu
date() ## [1] "Thu Nov 16 13:48:25 2017"Options
suppressMessages(library(knitr))
suppressMessages(library(rgl))
opts_knit$set(verbose = TRUE)
knit_hooks$set(webgl = hook_webgl)Load Packages:
suppressMessages(library(INLA))
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(fields))
suppressMessages(library(FactoMineR))
suppressMessages(library(adegenet))
suppressMessages(library(factoextra))
suppressMessages(library(corrplot))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(Rmisc))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(ggmap))
suppressMessages(library(ggthemes))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(xtable))
suppressMessages(library(PresenceAbsence))
suppressMessages(library(dismo))
suppressMessages(library(dplyr))Set Directory
setwd("D:/Miller_Data_100917")Vegetation Data
DVeg = read.csv("./Occurence.csv",
header = TRUE, sep=",")
dim(DVeg)## [1] 5537 78head(DVeg)## year habitat rep row col ALGALMAT ANDGLO ANDGYR ARIPUR BACSPP BUCAME
## 1 1999 B 1 A 1 0 0 0 0 0 0
## 2 1999 F 1 A 1 0 0 0 0 0 0
## 3 1999 I 1 A 1 0 0 0 0 0 0
## 4 1999 B 2 A 1 0 0 0 0 0 0
## 5 1999 F 2 A 1 0 0 0 0 0 0
## 6 1999 I 2 A 1 0 0 0 0 0 0
## BUCFLO BULCIL CAKCON CENINC CETASI CHAMAC CNISTI COMERE CONCAN CROPUN
## 1 0 0 0 0 0 0 0 0 1 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0
## CYAANG CYPCRO CYPLEC DICCOL DIOTER ELEFLA EREHIE ERALUG EUSPET FIMSPP
## 1 0 0 0 0 0 0 0 0 0 2
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 0 0 0 0 0
## FUISCI GALHIS HEDUNI HETSUB HYDBON HYPGEN IPOIMP IPOSPA IVAIMB JUNEFF
## 1 0 0 4 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 32 0 0 0 0 0 0
## 4 0 0 0 0 3 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 8 3 0 0 0 0 0
## JUNMEG JUNROE JUNSCI LICHEN LINCAR LINMED LIPMIC LUDMAR MUHCAP MYRCER
## 1 0 0 0 0 0 0 0 0 15 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 5 0 0 0 0 0 0 0 0 0
## OENHUM OPUHUM PANACI PANAMA PARERE PASDIS PHSANG PHYNOD POLPRO POLPUN
## 1 1 0 0 2 0 0 0 0 20 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 15 0 0 0 0 0 0 0 0 0
## 4 1 0 0 0 0 0 0 0 2 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 30 0 0 0 0
## RHESPP SABSTE SAGGRA SCHMAR SCHTAB SCIAME SCLSPP SESMAR SETSPP SMIAUR
## 1 0 2 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 23 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0
## SPAPAT SPOVIR STISET TRIPUR UNIPAN XYRCAR comments
## 1 0 0 0 0 0 0 spapat to muchcap
## 2 0 0 0 0 30 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 5 0
## 5 0 0 0 0 45 0
## 6 0 0 0 0 0 0 is triplasis PASDIS??Dropping New Plots
Replicate = 3
DVeg = DVeg %>%
filter(rep != 3)
dim(DVeg)## [1] 4704 78Unique Location Identifier
DVeg$Plot = paste0(DVeg$habitat, DVeg$rep, DVeg$row, DVeg$col, sep = "")
length(levels(as.factor(DVeg$Plot)))## [1] 294Get Plot Coordinates
DV.coords = read.csv("./Plot_coords.csv",
header = TRUE, sep=",")
DV.coords$Plot = paste0(DV.coords$HABITAT, DV.coords$REP, DV.coords$ROW, DV.coords$COL, sep = "")
dim(DV.coords)## [1] 294 10head(DV.coords)## IDX X_COORD Y_COORD REP ROW HABITAT COL Long Lat Plot
## 1 1 -37.195 -5.599 1 A B 1 332831 640361 B1A1
## 2 2 -27.341 -6.258 1 A B 2 332841 640360 B1A2
## 3 3 -17.123 -6.940 1 A B 3 332851 640360 B1A3
## 4 4 -7.056 -7.647 1 A B 4 332861 640359 B1A4
## 5 5 3.012 -8.293 1 A B 5 332871 640358 B1A5
## 6 6 12.983 -9.993 1 A B 6 332881 640357 B1A6Get Shape File of SGI
SGI = readOGR(dsn = "./ArcMap/Island",
layer = "SGIsland",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "./ArcMap/Island", layer: "SGIsland"
## with 1 features
## It has 1 fieldsdd = DV.coords[,8:9]
p.pnts = SpatialPointsDataFrame(dd, DV.coords)
proj4string(p.pnts) = proj4string(SGI)Quick Map
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
SGI2 = spTransform(SGI, LL84)
p.pnts2 = spTransform(p.pnts, LL84)
p.pnts2.df = p.pnts2@data %>%
mutate(nLong = p.pnts2@coords[,1],
nLat = p.pnts2@coords[,2])
colScale = scale_colour_manual(name = "Dune", values = c("red", "green", "orange"),
labels = c("Backdune", "Foredune", "Interdune"))
map = get_map(location = bbox(p.pnts2), maptype = "satellite", zoom = 13)## Warning: bounding box given to google - spatial extent only approximate.## converting bounding box to center/zoom specification. (experimental)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=29.767778,-84.694176&zoom=13&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=falseFar.plt = ggmap(map) +
geom_point(aes(x = nLong, y = nLat, col = factor(HABITAT)),
data = p.pnts2.df, pch =19) +
colScale +
xlab("Longitude") +
ylab("Latitude") +
theme(legend.title = element_text(size=16, face="bold"),
axis.title.x = element_text(size=18, face="bold"),
axis.title.y = element_text(size=18, face="bold"),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
plot.title = element_blank())
map = get_map(location = bbox(p.pnts2), maptype = "satellite", zoom = 17)## Warning: bounding box given to google - spatial extent only approximate.## converting bounding box to center/zoom specification. (experimental)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=29.767778,-84.694176&zoom=17&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=falsezoom.plt = ggmap(map) +
geom_point(aes(x = nLong, y = nLat, col = factor(HABITAT)),
data = p.pnts2.df, pch =19) +
colScale +
xlab("Longitude") +
ylab("Latitude") +
theme(legend.title = element_text(size=16, face="bold"),
axis.title.x = element_text(size=18, face="bold"),
axis.title.y = element_text(size=18, face="bold"),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
plot.title = element_blank())
grid.arrange(Far.plt, zoom.plt, ncol=2)
Changing to presence/absence and selecting 30 most abundant species
DVeg2 = DVeg[,6:77]
DVeg2[DVeg2>=1] = 1
tmp.vec = as.data.frame(colSums(DVeg2))
names(tmp.vec) = "Count"
tmp.vec$Names = rownames(tmp.vec)
tmp.vec = arrange(tmp.vec, desc(Count))
Keepers = as.character(tmp.vec$Names[1:30])
DVeg2 = DVeg2[,Keepers]
DVeg2 = DVeg2 %>%
mutate(Plot = DVeg$Plot,
Year = DVeg$year,
Habitat = DVeg$habitat)View Counts
ggplot(data=tmp.vec) +
geom_bar(aes(x=factor(Names), y=Count), stat="identity") +
ylab("Count") +
xlab("30 Most Abundant Species") +
theme_classic() +
guides(col = guide_legend(ncol = 16)) +
theme(axis.title.x=element_text(size=10, face="bold"),
axis.text.x=element_text(size=8, face="bold", vjust=1, hjust=1, angle=90),
axis.text.y=element_text(size=8, face="bold"),
axis.ticks.x=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=10, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
legend.title = element_text(size=10, face="bold"),
legend.position="bottom")
Calculating Richness by Year and Habitat
DTypes = c("F", "I", "B")
DVeg2$Year = as.character(DVeg2$Year)
SurvYears = levels(factor(as.character(DVeg2$Year)))
for(i in 1:3){
Tmp.df = DVeg2 %>%
filter(Habitat == DTypes[i])
Tmp.df$pRich = rowSums(Tmp.df[,1:30])
for(j in 1:length(SurvYears)){
Yr.Tmp.df = Tmp.df %>%
filter(Year == SurvYears[j])
tmp.CI = CI(Yr.Tmp.df$pRich , ci = 0.95)
if(j == 1){out.h = tmp.CI}
else(out.h = rbind(out.h, tmp.CI))}
if(i == 1){out.df = out.h}
else{out.df = cbind(out.df, out.h)}
}
out.df = as.data.frame(out.df)
rownames(out.df) = 1:16
names(out.df) = c(paste("F_", c("up", "m", "low"), sep=""),
paste("I_", c("up", "m", "low"), sep=""),
paste("B_", c("up", "m", "low"), sep=""))
out.df$Year = SurvYearsMean Plot Richness
ggplot(out.df, aes(factor(Year), F_m)) +
geom_point(size=2, pch=1, col = "black") +
geom_linerange(aes(ymin=F_low, ymax=F_up),
colour="black") +
geom_line(group=1, lty = "solid") +
geom_text(label="Foredune", x= 3, y=1.75, cex=8) +
geom_point(data = out.df, aes(factor(Year), I_m),
size=2, pch=1, col = "black") +
geom_linerange(data = out.df, aes(ymin=I_low, ymax=I_up),
colour="black") +
geom_line(data = out.df, aes(factor(Year), I_m),
group=1, lty = "dotted") +
geom_text(label="Interdune", x= 12, y=4, cex=8) +
geom_point(data = out.df, aes(factor(Year), B_m),
size=2, pch=1, col = "black") +
geom_linerange(data = out.df, aes(ymin=B_low, ymax=B_up),
colour="black") +
geom_line(data = out.df, aes(factor(Year), B_m),
group=1, lty = "dashed") +
geom_text(label="Backdune", x= 6, y=5, cex=8) +
theme_classic() +
xlab("Year") +
ylab("Mean Plot Richness") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=16),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14, vjust=1, hjust=1, angle=45))
Species Prevelance
Percent of plots with species.
for(i in 1:3){
Tmp.df = DVeg2 %>%
filter(Habitat == DTypes[i])
for(j in 1:length(SurvYears)){
Yr.Tmp.df = Tmp.df %>%
filter(Year == SurvYears[j])
Plt_N = nrow(Yr.Tmp.df)
Sp_total = colSums(Yr.Tmp.df[,1:30])
Perc.occ = as.data.frame(Sp_total/Plt_N)
names(Perc.occ) = "Perc"
Perc.occ$Dune = DTypes[i]
Perc.occ$Year = SurvYears[j]
Perc.occ$Species = rownames(Perc.occ)
if(j == 1){out.h = Perc.occ}
else(out.h = rbind(out.h, Perc.occ))}
if(i == 1){Pocc.df = out.h}
else{Pocc.df = rbind(Pocc.df, out.h)}
}
Pocc.df$Idx = 1:nrow(Pocc.df)Plot Mean Occupancy
#Foredune
DuneFU.sub = Pocc.df %>%
filter(Dune == "F",
Species == "UNIPAN")
DuneFO.sub = Pocc.df %>%
filter(Dune == "F",
Species == "OENHUM")
DuneFH.sub = Pocc.df %>%
filter(Dune == "F",
Species == "HETSUB")
DF.plt = ggplot(DuneFU.sub, aes(factor(Year), Perc)) +
geom_point(size=2, pch=1, col = "black") +
geom_line(group=1, lty = "dotted", col="black") +
geom_text(label="UNIPAN", x= 3, y=0.6, cex=6) +
ylim(0,1) +
geom_point(data = DuneFO.sub, aes(factor(Year), Perc),
size=2, pch=1, col = "black") +
geom_line(data = DuneFO.sub, aes(factor(Year), Perc),
group=1, lty = "dashed") +
geom_text(label="HETSUB", x= 7, y=0.30, cex=6) +
geom_point(data = DuneFH.sub, aes(factor(Year), Perc),
size=2, pch=1, col = "black") +
geom_line(data = DuneFH.sub, aes(factor(Year), Perc),
group=1, lty = "solid") +
geom_text(label="OENHUM", x= 11, y=0.15, cex=6) +
theme_classic() +
xlab("Year") +
ylab(NULL) +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=16),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14,
vjust=1, hjust=1, angle=45))
#Interdune
DuneFU.sub = Pocc.df %>%
filter(Dune == "I",
Species == "PASDIS")
DuneFO.sub = Pocc.df %>%
filter(Dune == "I",
Species == "PHYNOD")
DuneFH.sub = Pocc.df %>%
filter(Dune == "I",
Species == "HYDBON")
DI.plt = ggplot(DuneFU.sub, aes(factor(Year), Perc)) +
geom_point(size=2, pch=1, col = "black") +
geom_line(group=1, lty = "dotted", col="black") +
geom_text(label="PASSPP", x= 3, y=0.5, cex=6) +
ylim(0,1) +
geom_point(data = DuneFO.sub, aes(factor(Year), Perc),
size=2, pch=1, col = "black") +
geom_line(data = DuneFO.sub, aes(factor(Year), Perc),
group=1, lty = "dashed") +
geom_text(label="PHYNOD", x= 7, y=0.30, cex=6) +
geom_point(data = DuneFH.sub, aes(factor(Year), Perc),
size=2, pch=1, col = "black") +
geom_line(data = DuneFH.sub, aes(factor(Year), Perc),
group=1, lty = "solid") +
geom_text(label="HYDBON", x= 12, y=0.11, cex=6) +
theme_classic() +
xlab("Year") +
ylab("Percent of Plots Occupied") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=16),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14,
vjust=1, hjust=1, angle=45))
#Backdune
DuneFU.sub = Pocc.df %>%
filter(Dune == "B",
Species == "SCHMAR")
DuneFO.sub = Pocc.df %>%
filter(Dune == "B",
Species == "ERALUG")
DuneFH.sub = Pocc.df %>%
filter(Dune == "B",
Species == "HETSUB")
DB.plt = ggplot(DuneFU.sub, aes(factor(Year), Perc)) +
geom_point(size=2, pch=1, col = "black") +
geom_line(group=1, lty = "dotted", col="black") +
geom_text(label="SCHMAR", x= 3, y=0.15, cex=6) +
ylim(0,1) +
geom_point(data = DuneFO.sub, aes(factor(Year), Perc),
size=2, pch=1, col = "black") +
geom_line(data = DuneFO.sub, aes(factor(Year), Perc),
group=1, lty = "dashed") +
geom_text(label="ERALUG", x= 8, y=0.65, cex=6) +
geom_point(data = DuneFH.sub, aes(factor(Year), Perc),
size=2, pch=1, col = "black") +
geom_line(data = DuneFH.sub, aes(factor(Year), Perc),
group=1, lty = "solid") +
geom_text(label="HETSUB", x= 8, y=0.3, cex=6) +
theme_classic() +
xlab("Year") +
ylab(NULL) +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=16),
axis.text.y = element_text(face="bold", size=12),
axis.text.x = element_text(face="bold", size=14,
vjust=1, hjust=1, angle=45))
grid.arrange(DF.plt, DI.plt, DB.plt, ncol=1)
Organize Data
Also, looking-up coordinates by plot ID
DVeg3 = DVeg2
DVeg3$pRich = rowSums(DVeg3[,1:30])
#DVeg3 = DVeg3 %>%
# filter(pRich !=0)
LuTable = as.data.frame(cbind(p.pnts$Long, p.pnts$Plot))
DVeg3$nLong = with(LuTable,
V1[match(DVeg3$Plot,
V2)])
LuTable = as.data.frame(cbind(p.pnts$Lat, p.pnts$Plot))
DVeg3$nLat = with(LuTable,
V1[match(DVeg3$Plot,
V2)])
LuTable = as.data.frame(cbind(1:16, SurvYears))
DVeg3$Time = with(LuTable,
V1[match(DVeg3$Year,
SurvYears)])
DVeg3$nLong = as.numeric(as.character(DVeg3$nLong))
DVeg3$nLat = as.numeric(as.character(DVeg3$nLat))
DVeg3$Time = as.integer(as.character(DVeg3$Time))Get Terrain Variables
DEM = raster("./DEM/DEM.tif")
tpi = terrain(DEM,
neighbors = 8,
opt = 'tpi')
slope = terrain(DEM,
neighbors = 8,
opt = 'slope')
tri = terrain(DEM,
neighbors = 8,
opt = 'tri')
trr = terrain(DEM,
neighbors = 8,
opt = 'roughness')
Fdir = terrain(DEM,
neighbors = 8,
opt = 'flowdir')
aspect = terrain(DEM,
neighbors = 8,
opt = 'aspect')
Terra.stk = stack(DEM, tpi, tri, trr, slope, Fdir, aspect)
names(Terra.stk) = c("DEM", "tpi", "tri", "trr", "slope", "Fdir", "aspect")Observations to Points
DVeg3p = SpatialPointsDataFrame(cbind(DVeg3$nLong, DVeg3$nLat), DVeg3)
proj4string(DVeg3p) = proj4string(SGI)
Terra.df = as.data.frame(extract(Terra.stk, DVeg3p, method="simple"))
TCor = cor(Terra.df)
corrplot(TCor)
DVeg3p@data = cbind(DVeg3p@data, Terra.df)
DVeg3 = cbind(DVeg3, Terra.df)Compare Terrain Variables
Bxp.df = Terra.df
Bxp.df$Hab = DVeg3$Habitat
melt.bx = melt(Bxp.df, "Hab")
ggplot(melt.bx, aes(x=Hab, y=value, fill = Hab)) +
geom_boxplot() +
ylab("Terrain Factors") +
facet_wrap(~variable, scales = "free", ncol = 4) +
theme_classic() +
guides(col = guide_legend(ncol = 16)) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_text(size=12, face="bold"),
axis.ticks.x=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=12, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
legend.title = element_blank(),
legend.position="bottom")
Climate (PRISM)
files = list.files(path="./PRISM/PRISM_ppt_stable_4kmM3_198101_201704_bil",
pattern=".bil$", full.names=TRUE)
for(i in 1:length(SurvYears)){
y.files = mean(stack(files[grep(SurvYears[i], files)]))
if(i==1){mPrec.stk = y.files}
else{mPrec.stk = stack(mPrec.stk, y.files)}
}
names(mPrec.stk) = SurvYears
files = list.files(path="./PRISM/PRISM_tmean_stable_4kmM2_198101_201704_bil",
pattern=".bil$", full.names=TRUE)
for(i in 1:length(SurvYears)){
y.files = mean(stack(files[grep(SurvYears[i], files)]))
if(i==1){mT.stk = y.files}
else{mT.stk = stack(mT.stk, y.files)}
}
names(mT.stk) = SurvYearsPRISM to Observations
By year
Point.tmp = subset(DVeg3p, Plot=="B1A1" & Year=="1999")
Point.tmp = spTransform(Point.tmp, proj4string(mPrec.stk))
mPrec.tab = as.numeric(extract(mPrec.stk, Point.tmp, method = "bilinear"))
LuTableP = as.data.frame(cbind(mPrec.tab, SurvYears))
DVeg3$mPrec = with(LuTableP,
mPrec.tab[match(DVeg3$Year,
SurvYears)])
mT.tab = as.numeric(extract(mT.stk, Point.tmp, method = "bilinear"))
LuTableT = as.data.frame(cbind(mT.tab, SurvYears))
DVeg3$mTemp = with(LuTableT,
mT.tab[match(DVeg3$Year,
SurvYears)])
DVeg3$mPrec = as.numeric(as.character(DVeg3$mPrec))
DVeg3$mTemp = as.numeric(as.character(DVeg3$mTemp))Distance from Coastline
CoastLine = readOGR(dsn = "./ArcMap/Coastline",
layer = "Coast_line",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "./ArcMap/Coastline", layer: "Coast_line"
## with 1 features
## It has 1 fieldsplot(DEM)
plot(CoastLine, add=T, col="red", lwd=3)
suppressMessages(library(spatstat))
C.psp = as.psp(CoastLine)
P.pp = as.ppp(DVeg3p)
DVeg3p$CDist = nncross(P.pp, C.psp)[,"dist"]
DVeg3$CDist = DVeg3p$CDist
DVeg3 %>%
group_by(Habitat) %>%
summarize(Mean_Dist = mean(CDist))## # A tibble: 3 x 2
## Habitat Mean_Dist
## <fctr> <dbl>
## 1 B 425.1030
## 2 F 212.8868
## 3 I 319.0181Correlations
Cor.pa = cor(DVeg3[,c(1:30, 38:47)])
corrplot(Cor.pa)
Discriminant Analysis
Optimize alpha
Estimating number of components to include in the DAPC.
dapc.df = DVeg3[,c(1:30, 38:46)]
dapc1 = dapc(dapc.df,
DVeg2$Habitat,
n.da=15, n.pca=39)
temp = a.score(dapc1)
temp$pop.score## B F I
## 0.4726403 0.4775510 0.4120536temp$mean## [1] 0.4540816optim.a.score(dapc1, plot=FALSE)## $pop.score
## $pop.score$`39`
## B F I
## 0.4547832 0.4974490 0.4107781
##
## $pop.score$`1`
## B F I
## 0.14942602 -0.03852041 0.22008929
##
## $pop.score$`5`
## B F I
## 0.2290816 0.2140306 0.4089286
##
## $pop.score$`10`
## B F I
## 0.3436224 0.2447066 0.4733418
##
## $pop.score$`15`
## B F I
## 0.3837372 0.4065051 0.3894770
##
## $pop.score$`20`
## B F I
## 0.4114158 0.4320153 0.3704082
##
## $pop.score$`25`
## B F I
## 0.4294643 0.4772959 0.3666454
##
## $pop.score$`30`
## B F I
## 0.4100128 0.5051020 0.3970026
##
## $pop.score$`35`
## B F I
## 0.4514031 0.4957270 0.4105230
##
##
## $mean
## 39 1 5 10 15 20 25
## 0.4543367 0.1103316 0.2840136 0.3538903 0.3932398 0.4046131 0.4244685
## 30 35
## 0.4373724 0.4525510
##
## $pred
## $pred$x
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
##
## $pred$y
## [1] 0.1103316 0.1601070 0.2073428 0.2494935 0.2840136 0.3091641 0.3264327
## [8] 0.3381138 0.3465015 0.3538903 0.3620712 0.3708220 0.3794173 0.3871317
## [15] 0.3932398 0.3972587 0.3996757 0.4012205 0.4026231 0.4046131 0.4077248
## [22] 0.4117100 0.4161249 0.4205257 0.4244685 0.4276344 0.4302032 0.4324794
## [29] 0.4347676 0.4373724 0.4404925 0.4439020 0.4472694 0.4502630 0.4525510
## [36] 0.4539007 0.4544742 0.4545325 0.4543367
##
##
## $best
## [1] 38Conduct DAPC
DAPC.1 = dapc(dapc.df,
DVeg2$Habitat,
n.da=15, n.pca=39,
center = TRUE, scale = TRUE)
summary(DAPC.1)## $n.dim
## [1] 2
##
## $n.pop
## [1] 3
##
## $assign.prop
## [1] 0.8278061
##
## $assign.per.pop
## B F I
## 0.8233418 0.8743622 0.7857143
##
## $prior.grp.size
##
## B F I
## 1568 1568 1568
##
## $post.grp.size
##
## B F I
## 1385 1888 1431PCA Scatter
DAPCA.df = as.data.frame(DAPC.1$tab)
colnames(DAPCA.df) = paste("PC", seq(1,30,1), sep="")
DAPCA.df$Hab = as.factor(DVeg2$Habitat)
myShapes = c(16:18)
names(myShapes) = levels(DAPCA.df$Hab)
shpScale = scale_shape_manual(values = myShapes)
myCol = c("black", "red", "cyan")
names(myCol) = levels(DAPCA.df$Hab)
colScale = scale_colour_manual(name = "Habitat", values = myCol,
labels = c("Backdune", "Foredune", "Interdune"))
Labels = c("Backdune", "Foredune", "Interdune")
ggplot(DAPCA.df , aes(PC1, PC2,
col= DAPCA.df$Hab,
pch= DAPCA.df$Hab)) +
geom_point(size=2) +
ylab("PC2") +
xlab("PC1") +
scale_shape_manual(values = myShapes, labels = Labels) +
scale_colour_manual(values = myCol, labels = Labels) +
stat_ellipse(geom="polygon", show.legend = FALSE, alpha=0.2) +
theme_classic() +
guides(col = guide_legend(ncol = 16)) +
labs(shape = "Position", color = "Position", fill = "Position") +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size=12),
legend.position="bottom")
Linear Discriminant Scatter
DAPCA.df = as.data.frame(DAPC.1$ind.coord)
DAPCA.df$Hab = as.factor(DVeg2$Habitat)
myShapes = c(16:18)
names(myShapes) = levels(DAPCA.df$Hab)
shpScale = scale_shape_manual(values = myShapes)
myCol = c("black", "red", "cyan")
names(myCol) = levels(DAPCA.df$Hab)
colScale = scale_colour_manual(name = "Habitat", values = myCol,
labels = c("Backdune", "Foredune", "Interdune"))
Labels = c("Backdune", "Foredune", "Interdune")
ggplot(DAPCA.df , aes(LD1, LD2,
col= DAPCA.df$Hab,
pch= DAPCA.df$Hab)) +
geom_point(size=2) +
ylab("LD2") +
xlab("LD1") +
scale_shape_manual(values = myShapes, labels = Labels) +
scale_colour_manual(values = myCol, labels = Labels) +
stat_ellipse(geom="polygon", show.legend = FALSE, alpha=0.2) +
theme_classic() +
guides(col = guide_legend(ncol = 16)) +
labs(shape = "Position", color = "Position", fill = "Position") +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size=12),
legend.position="bottom")
Correct Assignment
temp = as.data.frame(summary(dapc(dapc.df, DVeg2$Habitat, n.da=15, n.pca=31))$assign.per.pop*100)
names(temp) = "Value"
myCol = c("black", "red", "cyan")
temp$Labels = c("Backdune", "Foredune", "Interdune")
#barplot(temp, xlab="% of correct reassignment to actual species",
# horiz=TRUE, las=1)
Reass.plt = ggplot(temp, aes(factor(Labels), Value, fill = factor(Labels), alpha=0.75)) +
geom_bar(stat="identity") +
scale_fill_manual(values = myCol) +
coord_flip() +
ylab("% Correct Plot Reassignment") +
theme_classic() +
guides(col = guide_legend(ncol = 3)) +
theme(axis.title.x=element_text(size=16, face="bold"),
axis.text.x=element_text(size=16, face="bold"),
axis.text.y=element_text(size=16, face="bold"),
axis.ticks.x=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=16, face="bold"),
axis.title.y = element_blank(),
legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size=12),
panel.grid.major.x = element_line(size = 1, color = "darkgray"),
panel.grid.minor.x = element_line(size = 0.5, color = "darkgray"),
legend.position="none")
Reass.plt
Probabilty of Group Assignment
Comp.df = as.data.frame(DAPC.1$posterior)
Comp.df$Hab = DVeg2$Habitat
Comp.df = arrange(Comp.df, Hab)
Comp.df$Idx = 1:nrow(Comp.df)
Comp.df = Comp.df[,c(1:3,5)]
Comp.dfm = melt(Comp.df, "Idx")
names(myCol) = levels(DVeg2$Habitat)
colScale = scale_fill_manual(name = "Habitat", values = myCol,
labels = c("Backdune", "Foredune", "Interdune"))
Comp.plt = ggplot(Comp.dfm, aes(factor(Idx), value, fill = variable)) +colScale+
geom_bar(stat="identity") +
ylab("Membership Probability") +
xlab("Plot") +
theme_classic() +
guides(fill = guide_legend(ncol = 8)) +
theme(axis.title.x=element_text(size=16, face="bold"),
axis.text.x=element_blank(),
axis.text.y=element_text(size=16, face="bold"),
axis.ticks.x=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size=12),
legend.position="bottom")
Comp.plt
ID Plots Below 50%
temp = which(apply(DAPC.1$posterior,1, function(e) all(e<0.5)))
length(temp)## [1] 91Comp.df = as.data.frame(DAPC.1$posterior)
Comp.df = Comp.df[temp,]
Comp.df$Idx = 1:nrow(Comp.df)
Comp.dfm = melt(Comp.df, "Idx")
Comp.plt = ggplot(Comp.dfm, aes(factor(Idx), value, fill = variable)) +
colScale +
geom_bar(stat="identity") +
ylab("Membership Probability") +
xlab("Plot") +
theme_classic() +
guides(col = guide_legend(ncol = 8)) +
theme(axis.title.x=element_text(size=12, face="bold"),
axis.text.x=element_blank(),
axis.text.y=element_text(size=12, face="bold"),
axis.ticks.x=element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=12, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
legend.title = element_text(size=12, face="bold"),
legend.position="bottom")
Comp.plt
Clustering by Plot
suppressMessages(library(phytools))
suppressMessages(library(phangorn))
Clust.df = DVeg3 %>%
filter(Year == "2015")
rownames(Clust.df) = Clust.df$Plot
Clust.df1 = Clust.df[,1:30]
Out = rep(1, 30)
Clust.df1 = rbind(Out, Clust.df1)
#rownames(Clust.df[1,]) = "Out"
Clust.mat = as.matrix(Clust.df1)
Clust2 = phyDat(Clust.mat, type="USER", levels=c("0","1"))
dm = dist.hamming(Clust2)
treeUPGMA = upgma(dm)Plot Composition Tree
suppressMessages(library(ggtree))
seq = treeUPGMA$tip.label #Create vector of tip labels
Tlab = data.frame(seq)
Tlab$Dune = substring(treeUPGMA$tip.label, 1, 2)
myCol = c("black", "red", "red4", "blue", "blue4", "green", "green4")
names(myCol) = levels(factor(Tlab$Dune))
colScale = scale_color_manual(name = "Dune", values = myCol,
labels = c("Outgroup", "Backdune 1", "Backdune 2", "Foredune 1",
"Foredune 2", "Interdune 1", "Interdune 2"))
p = ggtree(treeUPGMA, layout="circular")
p %<+% Tlab + geom_tippoint(aes(color = Dune),
alpha = 1, size = 4) + colScale +
geom_text2(aes(label = label,
subset = !isTip), size = 1.5, hjust = 1.2, vjust = -0.5) +
theme(legend.position = "right")
Spatio-temporal Modeling
Create Mesh
nProjKM = "+proj=aea +lat_1=24 +lat_2=31.5 +lat_0=24 +lon_0=-84 +x_0=400000
+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs"
SGI_km = spTransform(SGI, nProjKM)
p.pntsKM = spTransform(p.pnts, nProjKM)
bdry = inla.sp2segment(SGI_km)
MaxEdge = 0.05
mesh0 = inla.mesh.2d(#loc = p.pntsKM,
boundary = bdry,
cutoff = 0.01,
max.edge = MaxEdge,
min.angle = 25)
mesh1 = inla.mesh.2d(loc=mesh0$loc,
max.edge = MaxEdge*2,
offset = MaxEdge*2,
min.angle = 25)
mesh2 = inla.mesh.2d(loc=mesh1$loc,
max.edge = MaxEdge*4,
offset = MaxEdge*4,
min.angle = 25)
mesh2$n## [1] 2529plot(mesh2, rgl = F, col = "white", main = " ")
Collecting Mesh Verticies
mdd = as.data.frame(cbind(mesh2$loc[,1],
mesh2$loc[,2]))
names(mdd) = c("Long", "Lat")
v.pnts = SpatialPointsDataFrame(mdd, mdd)
proj4string(v.pnts) = proj4string(p.pntsKM)Distance from Coastline and Terrain
Terra.df2 = extract(Terra.stk, spTransform(v.pnts, proj4string(Terra.stk)), method="simple")
Terra.df2[is.na(Terra.df2)] = 0
v.pnts@data = cbind(v.pnts@data, Terra.df2)
#CDist
v.pnts.m = spTransform(v.pnts, proj4string(SGI))
P.pp = as.ppp(v.pnts.m)
v.pnts$CDist = nncross(P.pp, C.psp)[,"dist"]
#v.pnts$Ocean = !is.na(over(v.pnts, SGI_km))[,1]
#v.pnts = subset(v.pnts, Ocean == TRUE )Build Quadrature
Copy of quadrature for each year.
V1999 = as.data.frame(matrix(nrow = nrow(v.pnts@data), ncol = ncol(DVeg3)))
names(V1999) = names(DVeg3)
for(i in 1:length(SurvYears)){
Tmp.df = V1999 %>%
mutate(nLong = v.pnts$Long,
nLat = v.pnts$Lat,
DEM = v.pnts$DEM,
tpi = v.pnts$tpi,
tri = v.pnts$tri,
trr = v.pnts$trr,
slope = v.pnts$slope,
Fdir = v.pnts$Fdir,
aspect = v.pnts$aspect,
pRich = 0,
CDist = v.pnts$CDist,
Habitat = "Cosmo",
Year = SurvYears[i],
Time = paste(i))
if(i == 1){V.stack = Tmp.df}
else(V.stack = rbind(V.stack, Tmp.df))
}
V.stack$Plot = paste("Quad", 1:nrow(V.stack), sep="")
V.stack[,1:30] = 0
#Climate
V.stack$mPrec = with(LuTableP,
mPrec.tab[match(V.stack$Year,
SurvYears)])
V.stack$mTemp = with(LuTableT,
mT.tab[match(V.stack$Year,
SurvYears)])
V.stack$mPrec = as.numeric(as.character(V.stack$mPrec))
V.stack$mTemp = as.numeric(as.character(V.stack$mTemp))
#head(V.stack)Converting Projection
Coords4 = cbind(DVeg3$nLong, DVeg3$nLat)
DVeg4.pnts = SpatialPointsDataFrame(Coords4, DVeg3)
proj4string(DVeg4.pnts) = proj4string(SGI)
DVeg4.pnts = spTransform(DVeg4.pnts, proj4string(SGI_km))
DVeg4.pnts@data = DVeg4.pnts@data %>%
mutate(nLong = DVeg4.pnts@coords[,1],
nLat = DVeg4.pnts@coords[,2])Select Time Knots/Years
k = 5 #Number of time knots
DVeg4 = rbind(DVeg4.pnts@data) #, V.stack
DVeg4$Time = as.integer(DVeg4$Time)
DVeg4 = DVeg4 %>%
filter(Time <= k)Test for Multi-Linearity
suppressMessages(library(perturb))
Colin.df = DVeg4[,38:47]
CorCov = cor(Colin.df)
CI = colldiag(CorCov)
CI## Condition
## Index Variance Decomposition Proportions
## intercept DEM tpi tri trr slope Fdir aspect mPrec
## 1 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.002 0.000
## 2 1.925 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.021
## 3 2.328 0.003 0.000 0.064 0.000 0.000 0.000 0.001 0.007 0.000
## 4 2.921 0.016 0.000 0.023 0.000 0.000 0.000 0.016 0.000 0.005
## 5 3.885 0.004 0.000 0.001 0.000 0.000 0.000 0.004 0.085 0.012
## 6 4.795 0.000 0.001 0.000 0.002 0.000 0.001 0.039 0.000 0.031
## 7 5.201 0.001 0.000 0.009 0.001 0.000 0.002 0.007 0.186 0.040
## 8 21.728 0.158 0.545 0.150 0.094 0.001 0.008 0.196 0.141 0.152
## 9 32.643 0.542 0.001 0.035 0.896 0.040 0.045 0.503 0.384 0.495
## 10 64.667 0.276 0.452 0.719 0.007 0.958 0.945 0.233 0.196 0.244
## mTemp CDist
## 1 0.000 0.003
## 2 0.021 0.000
## 3 0.000 0.000
## 4 0.005 0.010
## 5 0.012 0.158
## 6 0.031 0.177
## 7 0.040 0.077
## 8 0.152 0.008
## 9 0.495 0.355
## 10 0.244 0.211Colin.df = Colin.df[,-c(2,7)]
CorCov = cor(Colin.df)
CI = colldiag(CorCov)
CI## Condition
## Index Variance Decomposition Proportions
## intercept DEM tri trr slope Fdir mPrec mTemp CDist
## 1 1.000 0.001 0.001 0.000 0.000 0.000 0.002 0.000 0.000 0.006
## 2 1.817 0.000 0.000 0.000 0.000 0.000 0.000 0.046 0.046 0.000
## 3 2.526 0.052 0.000 0.000 0.000 0.000 0.052 0.009 0.009 0.005
## 4 4.082 0.006 0.000 0.000 0.000 0.000 0.044 0.082 0.082 0.393
## 5 4.446 0.000 0.003 0.002 0.001 0.001 0.113 0.064 0.064 0.292
## 6 10.920 0.003 0.730 0.004 0.003 0.012 0.000 0.004 0.004 0.086
## 7 29.023 0.000 0.089 0.854 0.090 0.040 0.001 0.000 0.000 0.046
## 8 60.837 0.938 0.177 0.140 0.906 0.947 0.788 0.796 0.796 0.174Prior Specification for SPDE Model
locs = cbind(DVeg4$nLong, DVeg4$nLat)
A.mat = inla.spde.make.A(mesh=mesh2,
loc=locs,
group= DVeg4$Time)
spde1 = inla.spde2.pcmatern(mesh2, alpha = 2,
prior.range=c(0.01, 0.9),
prior.sigma=c(1, 0.5),
constr = TRUE)
field1 = inla.spde.make.index('field1', n.spde=spde1$n.spde, n.group=k)Error Precision
error.prec = list(hyper=list(theta=list(prior = "pc.prec", param=c(1, 0.01))))Construct Data Stack
DVeg4$Exp = ifelse(DVeg4$pRich >= 1, 0, 1)
FxSC.df = as.data.frame(apply(Colin.df, 2, scale))
apply(FxSC.df, 2, range)## DEM tri trr slope Fdir mPrec
## [1,] -0.5553666 -0.6722286 -0.7512135 -0.7229615 -1.297909 -1.751258
## [2,] 4.1097125 5.3730840 3.1926574 3.0611227 1.228619 1.241087
## mTemp CDist
## [1,] -1.573193 -1.630150
## [2,] 1.583402 1.743579FE.lst = list(c(field1,
list(intercept1 = 1)),
list(DEM = FxSC.df [,"DEM"],
CDist = FxSC.df[,"CDist"],
tri = FxSC.df [,"tri"],
trr = FxSC.df[,"trr"],
slope = FxSC.df [,"slope"],
Fdir = FxSC.df[,"Fdir"],
mPrec = FxSC.df[,"mPrec"],
mTemp = FxSC.df[,"mTemp"],
Year = DVeg4[,"Year"],
Plot = DVeg4[,"Plot"],
Habitat = DVeg4[,"Habitat"]))
Stack1 = inla.stack(data = list(SR = DVeg4$UNIPAN), # DVeg4$pRich
A = list(A.mat, 1),
effects = FE.lst,
tag = "rich.1")Construct and Execute Model
h.spec = list(theta=list(prior='pccor1', param=c(0, 0.9)))
HabPrior = list(theta=list(prior = "normal", param=c(0, 10)))
Frm0 = SR ~ -1 + intercept1 +
f(field1, model=spde,
group= field1.group,
control.group=list(model="ar1"),
hyper=h.spec)
#f(Habitat, model = "iid",
# hyper = HabPrior) +
#CDist + DEM + tri + slope + Fdir + mTemp + mPrec
#theta0 = Model0$internal.summary.hyperpar$mean
theta0 = c(-1.31, 1.48, 7.47) #previous run
Model0 = inla(Frm0,
data = inla.stack.data(Stack1, spde=spde1),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta0),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model0)##
## Call:
## c("inla(formula = Frm0, family = \"binomial\", data = inla.stack.data(Stack1, ", " spde = spde1), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.1, prec.intercept = 0.001), ", " control.mode = list(restart = TRUE, theta = theta0))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.3804 246.9816 1.0321 249.3941
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -0.7607 0.7683 -2.2697 -0.7605 0.7458 -0.76 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field1 0.0460 0.0143 0.0229 0.0445 0.0784 0.0414
## Stdev for field1 4.2549 0.6927 3.0810 4.1885 5.7994 4.0510
## GroupRho for field1 0.9966 0.0028 0.9894 0.9974 0.9995 0.9987
##
## Expected number of effective parameters(std dev): 63.61(0.00)
## Number of equivalent replicates : 23.11
##
## Deviance Information Criterion (DIC) ...: 1217.28
## Effective number of parameters .........: 60.95
##
## Watanabe-Akaike information criterion (WAIC) ...: 1213.75
## Effective number of parameters .................: 53.62
##
## Marginal log-Likelihood: -677.14
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedQuick Check
idat = inla.stack.index(Stack1, "rich.1")$data
cor(DVeg4$pRich, Model0$summary.linear.predictor$mean[idat])## [1] 0.2225647Joint Model
Spatial Indicies
#field1 = inla.spde.make.index('field1', n.spde=spde1$n.spde, n.group=k)
s1 = s2 = s3 = s12 = s13 = s23 = rep(1:spde1$n.spde, times=k)
g1 = g2 = g3 = g12 = g13 = g23 = rep(1:k, each=spde1$n.spde)
rho1p = list(theta=list(prior='pccor1', param=c(0, 0.9)))
ctr.g = list(model='ar1', hyper=rho1p)Model Formula
hc3 = hc2 = hc1 = list(theta=list(prior='normal', param=c(0,10))) #prior
J.Form = Y ~ -1 + intercept1 +
intercept2 +
intercept3 +
f(s1, model=spde,
ngroup=k, group=g1,
control.group=ctr.g) +
f(s2, model=spde,
ngroup=k, group=g2,
control.group=ctr.g) +
f(s3, model=spde,
ngroup=k, group=g3,
control.group=ctr.g) +
f(s12, copy="s1",
group=g12,
fixed=FALSE) +
f(s13, copy="s1",
group=g13,
fixed=FALSE) +
f(s23, copy="s2",
group=g23,
fixed=FALSE)Projector Matrix
(same for all)
locs = cbind(DVeg4$nLong, DVeg4$nLat)
A.mat = inla.spde.make.A(mesh=mesh2,
loc=locs,
n.group=k,
group= DVeg4$Time)Organizing Data for Moidel Fitting
#HETSUB
stack1 = inla.stack(
data=list(Y=cbind(DVeg4$HETSUB, NA, NA)),
A=list(A.mat),
effects=list(list(intercept1=1, s1=s1, g1=g1)))
#ERALUG
stack2 = inla.stack(
data=list(Y=cbind(NA, DVeg4$ERALUG, NA)), A=list(A.mat),
effects=list(list(intercept2=1, s2=s2, g2=g2,
s12=s12, g12=g12)))
#UNIPAN (focal species)
stack3 = inla.stack(
data=list(Y=cbind(NA, NA, DVeg4$UNIPAN)), A=list(A.mat),
effects=list(list(intercept3=1, s3=s3, g3=g3,
s13=s13, g13=g13, s23=s23, g23=g23)))
Joint.stk = inla.stack(stack1, stack2, stack3)Execute Joint Multi-Species Model
Model.J = inla(J.Form,
data = inla.stack.data(Joint.stk, spde=spde1),
family = c("binomial", "binomial","binomial"),
verbose = TRUE,
control.fixed = list(prec = 0.1, prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) Model Comparison
Checking Parsimony
#Model0 (spatial effect only)
Model0$dic$dic## [1] 1217.283Model0$waic$waic## [1] 1213.746#Model.J (joint model)
tapply(Model.J$dic$local.dic, Model.J$dic$family, sum)[3]## 3
## 1217.154tapply(Model.J$waic$local.waic, Model.J$dic$family, sum)[3]## 3
## 1214.234Prediction Locations
Points for visualization
Tmp.grd = makegrid(SGI, cellsize = 3)
Grd.pnts = SpatialPointsDataFrame(Tmp.grd, Tmp.grd)
proj4string(Grd.pnts) = proj4string(SGI)
Grd.pnts$ID = !is.na(over(Grd.pnts, SGI))[,1]
Grd.pnts = subset(Grd.pnts, ID == TRUE )
#Terra
Terra.df2 = as.data.frame(extract(Terra.stk,
spTransform(Grd.pnts, proj4string(Terra.stk)), method="simple"))
#Dist
P.pp = as.ppp(Grd.pnts)
Terra.df2$CDist = as.numeric(nncross(P.pp, C.psp)[,"dist"])
Terra.df2[is.na(Terra.df2)] = 0
Terra.df2 = apply(Terra.df2, 2, scale)
Grd.pnts@data = cbind(Grd.pnts@data, Terra.df2)Pred Results for Plotting
F.s1 = cbind(Model.J$summary.random$s1$mean, field1$field1.group)
F.s2 = cbind(Model.J$summary.random$s2$mean, field1$field1.group)
F.s3 = cbind(Model.J$summary.random$s3$mean, field1$field1.group)
F.s12 = cbind(Model.J$summary.random$s12$mean, field1$field1.group)
F.s13 = cbind(Model.J$summary.random$s13$mean, field1$field1.group)
F.s23 = cbind(Model.J$summary.random$s23$mean, field1$field1.group)
s1xmean = list()
s1xmean = split(F.s1[,1], F.s1[,2])
s2xmean = list()
s2xmean = split(F.s2[,1], F.s2[,2])
s3xmean = list()
s3xmean = split(F.s3[,1], F.s3[,2])
s12xmean = list()
s12xmean = split(F.s12[,1], F.s12[,2])
s13xmean = list()
s13xmean = split(F.s13[,1], F.s13[,2])
s23xmean = list()
s23xmean = split(F.s23[,1], F.s23[,2])
Grd.pnts2 = spTransform(Grd.pnts, proj4string(SGI_km))
GLoc = cbind(Grd.pnts2@coords[,1], Grd.pnts2@coords[,2])
Ap = inla.spde.make.A(mesh2, loc=GLoc)
blank.r = raster(nrow = 400, ncol = 500, ext = extent(SGI_km))
Domain.r = rasterize(SGI_km, blank.r,
field = 0,
background = NA)
for(i in 1:k){
Grd.pnts2$m1RE = drop(Ap %*% xmean[[i]])
tmp.r = rasterize(Grd.pnts2,
Domain.r,
"m1RE",
background = NA)
if(i == 1){M0.stk = tmp.r}
else{M0.stk = stack(M0.stk, tmp.r)}
}
names(M0.stk) = SurvYears[1:k]
for(i in 1:k){
Grd.pnts2$m1RE = plogis(drop(Ap %*% s1xmean[[i]]) +
drop(Ap %*% s2xmean[[i]]) +
drop(Ap %*% s3xmean[[i]]) +
drop(Ap %*% s12xmean[[i]]) +
drop(Ap %*% s13xmean[[i]]) +
drop(Ap %*% s23xmean[[i]]) +
Model.J$summary.fixed[1,1] +
Model.J$summary.fixed[2,1] +
Model.J$summary.fixed[3,1])
tmp.r = rasterize(Grd.pnts2,
Domain.r,
"m1RE",
background = NA)
if(i == 1){M0p.stk = tmp.r}
else{M0p.stk = stack(M0p.stk, tmp.r)}
}
names(M0p.stk) = SurvYears[1:k]Random Field Density
rng = seq(-8.1, 8.1, 0.01)
mCols = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4","#313695")
#brewer.pal(11, "RdBu")
cr = colorRampPalette(c(rev(mCols)),
bias = 1, space = "rgb")
RF0osi = levelplot(M0.stk,
layout=c(2, 3),
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(cex=1.5),
space = "right"),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SGI, col = "black", lwd = 1))
RF0osi
Probabilty of UNIPAN occurence
rng = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrRd")[2:9]
cr0 = rev(colorRampPalette(rev(mCols))(n = 299))
cr = colorRampPalette(c("tan", cr0),
bias = 2, space = "rgb")
M1Jop = levelplot(M0p.stk,
margin = FALSE,
xlab = " ", ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(labels=list(at=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0%", "25%", "50%", "75%", "100%"), cex=1.5),
labels=list(cex=12),
space = "right"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5))
M1Jop
Session Info
date()## [1] "Thu Nov 16 13:50:03 2017"sessionInfo()## R version 3.4.2 (2017-09-28)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] perturb_2.05 ggtree_1.8.2 treeio_1.0.2
## [4] phangorn_2.3.1 phytools_0.6-44 ape_5.0
## [7] spatstat_1.53-2 rpart_4.1-11 nlme_3.1-131
## [10] spatstat.data_1.1-1 bindrcpp_0.2 dplyr_0.7.4
## [13] dismo_1.1-4 PresenceAbsence_1.1.9 xtable_1.8-2
## [16] gridExtra_2.3 GISTools_0.7-4 MASS_7.3-47
## [19] ggthemes_3.4.0 ggmap_2.6.1 spdep_0.6-15
## [22] mapproj_1.2-5 maptools_0.9-2 Rmisc_1.5
## [25] plyr_1.8.4 rgeos_0.3-26 rasterVis_0.41
## [28] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-35
## [31] corrplot_0.84 factoextra_1.0.5 ggplot2_2.2.1
## [34] adegenet_2.1.0 ade4_1.7-8 FactoMineR_1.39
## [37] fields_9.0 maps_3.2.0 spam_2.1-1
## [40] dotCall64_0.9-04 raster_2.5-8 reshape2_1.4.2
## [43] rgdal_1.2-15 INLA_17.08.19 Matrix_1.2-11
## [46] sp_1.2-5 rgl_0.98.1 knitr_1.17
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.3-2 seqinr_3.4-5
## [3] rjson_0.2.15 deldir_0.1-14
## [5] rprojroot_1.2 hexbin_1.27.1
## [7] ggrepel_0.7.0 mvtnorm_1.0-6
## [9] splines_3.4.2 leaps_3.0
## [11] mnormt_1.5-5 polyclip_1.6-1
## [13] jsonlite_1.5 cluster_2.0.6
## [15] png_0.1-7 shiny_1.0.5
## [17] compiler_3.4.2 rvcheck_0.0.9
## [19] backports_1.1.1 assertthat_0.2.0
## [21] lazyeval_0.2.1 htmltools_0.3.6
## [23] tools_3.4.2 igraph_1.1.2
## [25] coda_0.19-1 gtable_0.2.0
## [27] glue_1.2.0 clusterGeneration_1.3.4
## [29] gmodels_2.16.2 fastmatch_1.1-0
## [31] Rcpp_0.12.13 msm_1.6.4
## [33] gdata_2.18.0 stringr_1.2.0
## [35] proto_1.0.0 mime_0.5
## [37] gtools_3.5.0 goftest_1.1-1
## [39] LearnBayes_2.15 zoo_1.8-0
## [41] scales_0.5.0 spatstat.utils_1.7-1
## [43] parallel_3.4.2 expm_0.999-2
## [45] animation_2.5 geosphere_1.5-7
## [47] stringi_1.1.5 plotrix_3.6-6
## [49] permute_0.9-4 boot_1.3-20
## [51] RgoogleMaps_1.4.1 rlang_0.1.4
## [53] pkgconfig_2.0.1 evaluate_0.10.1
## [55] purrr_0.2.4 tensor_1.5
## [57] bindr_0.1 htmlwidgets_0.9
## [59] labeling_0.3 magrittr_1.5
## [61] R6_2.2.2 combinat_0.0-8
## [63] foreign_0.8-69 mgcv_1.8-20
## [65] survival_2.41-3 scatterplot3d_0.3-40
## [67] abind_1.4-5 tibble_1.3.4
## [69] rmarkdown_1.7 jpeg_0.1-8
## [71] vegan_2.4-4 digest_0.6.12
## [73] flashClust_1.01-2 numDeriv_2016.8-1
## [75] tidyr_0.7.2 httpuv_1.3.5
## [77] munsell_0.4.3 viridisLite_0.2.0
## [79] quadprog_1.5-5