Cougar GAM
# Read in data
setwd('C:/Users/Maggie/Desktop/USU_COUGAR_FOLDERS')
photo <- read.csv('cougar_photos.csv')
# Convert times to minute of the day (range: 0:1439)
time2 <- photo$hour*60 + photo$minute
photo <- cbind(photo,time2)
# Pull out minutes when PUCO detected and save as df
dets <- photo %>%
group_by(time2) %>%
summarize(dets=sum(PUCO))
dets <- as.data.frame(dets)
# fill in detections as one, instead of number of photos, per minute
dets$dets[1:16] <- rep(1,16)
# Create df including every minute with 0's for absence (no observations)
dets2 <- as.data.frame(matrix(NA,nrow=1441,ncol=2))
names(dets2) <- c("minute","detections")
for (i in 1:1441){
dets2$minute[i] = i-1
dets2$detections[i] = 0
}
# Change minutes with detections to 1
for (i in dets$time2){
dets2$detections[i+1] = 1
}
# Run basic GAM and plot results, plus lines for data points
gam <- gam(detections ~ s(minute,bs="cc"),data=dets2)
plot(gam)
abline(h=0,col="green")
for (i in dets$time2) abline(v=i,col="blue")
# Run predictions for full diel cyle and plot (code from Michel)
newdat <- data.frame(minute=seq(0,1440,1))
out <- predict.gam(gam, newdata=newdat,type="response",se.fit=TRUE,bs="cc")
out2 <- data.frame(minute=seq(0,1440,1),
fit=out$fit,
se=out$se.fit)
# generate C.I.s
out2$upCI <- out2$fit + (1.96* out2$se)
out2$lowCI <- out2$fit - (1.96* out2$se)
ggplot(out2) +
geom_line(aes(y=fit, x=minute)) + # predictions
geom_ribbon(aes(x=minute, ymin=lowCI, ymax=upCI), alpha=0.15) + # confidence band
geom_rug(sides="bl",data=dets,mapping=aes(x=time2)) + # detections
xlab("Minute of the diel cycle") +
ylab("Predicted GAM response") +
ggtitle("GAM of Predicted Cougar Activity throughout the Diel Cycle")
##
## Family: gaussian
## Link function: identity
##
## Formula:
## detections ~ s(minute, bs = "cc")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.011103 0.002753 4.033 5.81e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(minute) 4.453 8 1.263 0.0371 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.00573 Deviance explained = 0.88%
## GCV = 0.010966 Scale est. = 0.010925 n = 1441
Mule Deer GAM - presence/absence
# Pull out minutes when ODHE detected and save as df
ODHE <- read.csv("photo_data_21JAN18.csv")
time_min <- ODHE$Hour*60 + ODHE$minute
ODHE <- cbind(ODHE, time_min)
ODHE1 <- ODHE %>%
filter(Species1=="ODHE")
## Warning: package 'bindrcpp' was built under R version 3.4.3
ODHE2 <- rep(1, length(ODHE1$CreationDatetime))
ODHE1 <- cbind(ODHE1, ODHE2)
deer <- ODHE1 %>%
filter(!Season %in% c("SPR16","SUM16")) %>%
group_by(time_min) %>%
summarize(dets=sum(ODHE2))
deer <- as.data.frame(deer)
# Create df including every minute with 0's for absence (no observations)
deer2 <- as.data.frame(matrix(NA,nrow=1441,ncol=2))
names(dets2) <- c("minute","detections")
for (i in 1:1441){
deer2$minute[i] = i-1
deer2$detections[i] = 0
}
# Change minutes with detections to 1
for (i in deer$time_min){
deer2$detections[i+1] = 1
}
# Run basic GAM and plot results, plus lines for data points
gam2 <- gam(detections ~ s(minute,bs="cc"),data=deer2)
plot(gam2)
abline(h=0,col="green")
for (i in deer$time_min) abline(v=i,col="blue")
# Run predictions for full diel cyle and plot (code from Michel)
newdat <- data.frame(minute=seq(0,1440,1))
out <- predict.gam(gam2, newdata=newdat,type="response",se.fit=TRUE,bs="cc")
out3 <- data.frame(minute=seq(0,1440,1),
fit=out$fit,
se=out$se.fit)
# generate C.I.s
out3$upCI <- out3$fit + (1.96* out3$se)
out3$lowCI <- out3$fit - (1.96* out3$se)
ggplot(out3) +
geom_line(aes(y=fit, x=minute)) + # predictions
geom_ribbon(aes(x=minute, ymin=lowCI, ymax=upCI), alpha=0.15) + # confidence band
geom_rug(sides="bl",data=deer,mapping=aes(x=time_min)) + # detections
xlab("Minute of the diel cycle") +
ylab("Predicted GAM response") +
ggtitle("GAM of Predicted Mule Deer Activity - presence/absence")
##
## Family: gaussian
## Link function: identity
##
## Formula:
## detections ~ s(minute, bs = "cc")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.085357 0.006992 12.21 <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(minute) 7.378 8 19.97 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0982 Deviance explained = 10.3%
## GCV = 0.070869 Scale est. = 0.070457 n = 1441
Mule Deer GAM - detections
# Change minutes with detections to number detections
deer2[deer2$minute %in% deer$time_min, ]$detections <- deer$dets
# Run basic GAM and plot results, plus lines for data points
gam3 <- gam(detections ~ s(minute,bs="cc"),data=deer2)
plot(gam3)
abline(h=0,col="green")
for (i in deer$time_min) abline(v=i,col="blue")
# Run predictions for full diel cyle and plot (code from Michel)
newdat <- data.frame(minute=seq(0,1440,1))
out <- predict.gam(gam3, newdata=newdat,type="response",se.fit=TRUE,bs="cc")
out4 <- data.frame(minute=seq(0,1440,1),
fit=out$fit,
se=out$se.fit)
# generate C.I.s
out4$upCI <- out4$fit + (1.96* out4$se)
out4$lowCI <- out4$fit - (1.96* out4$se)
ggplot(out4) +
geom_line(aes(y=fit, x=minute)) + # predictions
geom_ribbon(aes(x=minute, ymin=lowCI, ymax=upCI), alpha=0.15) + # confidence band
geom_rug(sides="bl",data=deer,mapping=aes(x=time_min)) + # detections
xlab("Minute of the diel cycle") +
ylab("Predicted GAM response") +
ggtitle("GAM of Predicted Mule Deer Activity - detections")
##
## Family: gaussian
## Link function: identity
##
## Formula:
## detections ~ s(minute, bs = "cc")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.55586 0.06331 8.781 <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(minute) 7.133 8 12.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0641 Deviance explained = 6.88%
## GCV = 5.8077 Scale est. = 5.7749 n = 1441