Plotting an additive Unmarked occupancy model with ggplot

[Caleb Durbin]

Hello all, 

I am trying to figure out the best way to plot my additive model and I am not familiar ggplot and I would like to get a more interpretable plot than the one I have now. The water covariate is the distance to water and the woody covariate is the percentage of woody cover. 

Any help would be greatly appreciated! Thank you!

Best, 

Caleb

Here is the code along with the (ugly) plot: 

#  bookkeeping
rm(list=ls())

library(unmarked)
library(AICcmodavg)


setwd("D:/Hadley, SER, Kirwin, Lonsinger Data")
# ------------------------- Read data ------------------------------


# Import data and check structure
kancoon.data <- read.csv("ALL2023Kansasdata.csv", row.names=1)
str(kancoon.data)
#kancoon.yd2 <- kancoon.data[kancoon.data$detection == "detected",] #subsetted data of just detections
#kancoon.ynd <- kancoon.data[kancoon.data$detection == "not detected",] #subsetted data of just nondetections
# Pull out count matrix and covert to binary
kancoon.y <- kancoon.data[,4:39]
#coon.y <- cbind(deployment_id=coon.data[,2],coon.y)
kancoon.y1 <- kancoon.y # Make a copy


# Create a covariate called sitelccovs
sitecovs <- read.csv("KansasSitecovs.csv")

# Transect-specific covariates
month <- c(sitecovs$Month)
site <- c(sitecovs$Site)
water <- c(sitecovs$distance_to_allwatersources)
woody <- c(sitecovs$WC_2022_RAP_nlcdfix)

water.mean <- mean(water)
water.sd <- sd(water)
water.z <- (water-water.mean)/water.sd

woody.mean <- mean(woody)
woody.sd <- sd(woody)
woody.z <- (woody-woody.mean)/woody.sd



(covs <- data.frame(month = month, site = site, water = water, woody = woody))



# Create unmarkedFrame
library(unmarked)
kanraccoon.umf <- unmarkedFrameOccu(y=kancoon.y1,
                                     siteCovs=covs)
summary(kanraccoon.umf)



# the intercept model
#
(intercept <- occu(~1 ~1, kanraccoon.umf))

summary(intercept)

backTransform(intercept, type="state")
backTransform(intercept, type="det")

#predict(intercept, type='det')

Cand.modDet<- list()

####################DETECTION MODELS###################################################################################        
###### occu(~detection.covariate ~1, amwo.umf)####################################

Cand.modDet[[1]] <-occu(~1 ~1, kanraccoon.umf)
Cand.modDet[[2]] <-occu(~month ~1, kanraccoon.umf)
Cand.modDet[[3]] <-occu(~site ~1, kanraccoon.umf)
Cand.modDet[[4]] <-occu(~water ~1, kanraccoon.umf)
Cand.modDet[[5]] <-occu(~woody ~1, kanraccoon.umf)



ModnamesDet <- c("intercept1","month2","site3","water4","woody5")
               

##model selection table based on AICc
aictab(cand.set = Cand.modDet, modnames = ModnamesDet)
# top model was water but decided to keep the detection constant
##set up candidate models          
Cand.modPSI<- list()

#Abundance models
Cand.modPSI[[1]] <-occu(~1 ~1, kanraccoon.umf)
Cand.modPSI[[2]] <-occu(~1~water, kanraccoon.umf)
Cand.modPSI[[3]] <-occu(~1~woody, kanraccoon.umf)
Cand.modPSI[[4]] <-occu(~1~water+woody, kanraccoon.umf)
Cand.modPSI[[5]] <-occu(~1~water*woody, kanraccoon.umf)
Cand.modPSI[[6]] <-occu(~1~water+ I(water^2), kanraccoon.umf)



##assign names to each model
ModnamesPSI <- c("null_1","water_2","woody_3","water_woody_add_4","water_woody_inter_5","water_quad_6")
##model selection table based on AICc
aictab(cand.set = Cand.modPSI, modnames = ModnamesPSI)
####water+woody was the top-ranked model. ####
Cand.modPSI[[4]]




#occupancy with water +woody model 
(fm25 <- occu(~1~water+woody, kanraccoon.umf))


predwater= c(rep(50, times = 201))
predwater= append(predwater,rep(250, times = 201))
predwater= (water=append(predwater,rep(450, times = 201)))                

newData2<- data.frame(woody =rep(seq(min(woody),100, by=.5),times = 3))
newData2 <- cbind(newData2,water = predwater)
W.PSI <- predict(fm25, type="state", newdata=newData2, appendData=TRUE)
head(W.PSI)



# Plot it
par(family="serif")
plot(Predicted[1:201] ~ woody[1:201], W.PSI, type="l",
     xlab="Percentage of Woody Cover",
     ylab="Occupancy")

lines(lower[1:201] ~ woody[1:201], W.PSI, type="l", col=gray(0.5))
lines(upper[1:201] ~ woody[1:201], W.PSI, type="l", col=gray(0.5))

lines(Predicted[202:402] ~ woody[202:402], W.PSI, type="l", col="blue")
lines(upper[202:402] ~ woody[202:402], W.PSI, type="l", col="red")
lines(lower[202:402] ~ woody[202:402], W.PSI, type="l", col="red")


lines(Predicted[403:603] ~ woody[403:603], W.PSI, type="l", col="purple")
lines(upper[403:603] ~ woody[403:603], W.PSI, type="l", col="brown")
lines(lower[403:603] ~ woody[403:603], W.PSI, type="l", col="brown")

[Jim Baldwin]

If there are two predictors, then a contour plot with contours of equal occupancy probabilities would seem to be a good choice.  (Or a spinning 3D plot which would be a bit more exciting.)

Jim

--
*** Three hierarchical modeling email lists ***
(1) unmarked (this list): for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology: for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/d02ffc24-edc7-4220-8a3f-23476828f2c9n%40googlegroups.com.

[Caleb Durbin]

Thanks Jim!

I am interested in using the contour plot. I’m just not sure how to go about creating this plot and how to interpret the plot after it is created. I’m relatively new to R coding. If there is some code and/or articles to help me code this plot I would be very grateful to have that sent my way. 

I conducted single season occupancy models. If there is a way to make the confidence intervals a lighter color fill of the line color and the line color is bolder and thicker, then i think the graph I had might work but I do want the best graph I can get to display my results. Thank you again for the help! 

Best, 

Caleb

Caleb Durbin
AWB®

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/CAJCKiwzNHmO2vuRkmokM-2VntGXZ3iVAybzwRRbk76OkbJXfhw%40mail.gmail.com.

[Jim Baldwin]

Here is some code using the contour function (but not using ggplot2) and data from the unmarked documentation.

# Import data

library(unmarked)

wt <- read.csv(system.file("csv","widewt.csv", package="unmarked"))

y <- wt[,2:4]

siteCovs <-  wt[,c("elev", "forest", "length")]

wt <- unmarkedFrameOccu(y = y, siteCovs = siteCovs)

# Run model

fm <- occu(~ 1 ~ elev*length, wt)

# Generate list of increasing values for elev and length

ngrid <- 50   # Number of grid points

xvals <- min(siteCovs$elev, na.rm=TRUE) +

  (max(siteCovs$elev, na.rm=TRUE) - min(siteCovs$elev, na.rm=TRUE))*c(0:ngrid)/ngrid

yvals <- min(siteCovs$length, na.rm=TRUE) +

  (max(siteCovs$length, na.rm=TRUE) - min(siteCovs$length, na.rm=TRUE))*c(0:ngrid)/ngrid

# Generate matrix of occupancy probabilities

x.elev <- rep(xvals, ngrid+1)

y.length <- rep(yvals, each=ngrid+1)

zvals <- matrix(predict(fm, type="state", new=data.frame(elev=x.elev, length=y.length))$Predict, nrow=ngrid+1)

# Produce contour plot

contour(xvals, yvals, zvals, nlevels=10, las=1,

  xlab="Elevation", ylab="Length", main="Predicted occupancy probabilities")

Jim

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/CAMrqDrrY_yOmFi0vjdZUpzq3j3jmPhmoNtvzp2w3cOkPhKXEGA%40mail.gmail.com.