setwd("/Users/gdunshea/Documents/PaperAnalysis/ObserverEvaluation_V3/Rscripts&data") library(lattice) library(plyr) library(reshape) library(reshape2) gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting only marine mammal sightings in the transect levels(gpbin$transect) levels(gpbin$transect_position) levels(gpbin$realdt) recapFCcom1 <- subset(gpbin, realdt %in% c("Dolphin", "Dugong") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom1 <- droplevels(recapFCcom1) levels(recapFCcom1$transect) levels(recapFCcom1$transect_position) str(recapFCcom1) ## Got all dolphin and dugong sightings, now getting complete cases for analysis completeFun <- function(data, desiredCols) { completeVec <- complete.cases(data[, desiredCols]) return(data[completeVec, ]) } count(recapFCcom1, c("realdt", "disagno")) recapFCcom <- completeFun(recapFCcom1, c("realgroup", "realdt", "team1", "transect_position", "beaufort", "turbidity", "glare")) str(recapFCcom) count(recapFCcom, c("team1", "realdt", )) table(recapFCcom$team1, recapFCcom$turbidity) table(recapFCcom$team1, recapFCcom$realdt) count(recapFCcom, c( "glare","disagno" )) ### only one Turbidity = 1 sighting, initially dropping this to examine turbidity effects ### Also West STBD had no dugong sightings with a group size disagreement, so can't include team*datatype interaction ## Dropping single T=1 sightings recapFCcom <- subset(recapFCcom, turbidity !=1) str(recapFCcom) library(MuMIn) m31 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday, data=recapFCcom, family="binomial") m32 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+realgroup:team1, data=recapFCcom, family="binomial") m33 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(glare)+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m34 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(glare)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m35 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m36 <- glm(disagno~realgroup+team1+realdt+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m37 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m38 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m31,m32,m33,m34,m35,m36,m37, rank="AICc") m351 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(turbidity)+as.factor(beaufort)+realgroup:team1, data=recapFCcom, family="binomial") m352 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m353 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m354 <- glm(disagno~realgroup+team1+realdt+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m355 <- glm(disagno~realgroup+team1+transect_position+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m31,m32,m33,m34,m35,m36,m37,m351,m352,m353,m354,m355, rank="AICc") m3521 <- glm(disagno~realgroup+team1+realdt+transect_position+as.factor(turbidity)+realgroup:team1, data=recapFCcom, family="binomial") m3522 <- glm(disagno~realgroup+team1+realdt+transect_position+cumday+realgroup:team1, data=recapFCcom, family="binomial") m3523 <- glm(disagno~realgroup+team1+realdt+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m3524 <- glm(disagno~realgroup+team1+transect_position+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m31,m32,m33,m34,m35,m36,m37,m351,m352,m353,m354,m355,m3521,m3522,m3523,m3524, rank="AICc") m35231 <- glm(disagno~realgroup+team1+realdt+as.factor(turbidity)+realgroup:team1, data=recapFCcom, family="binomial") m35232 <- glm(disagno~realgroup+team1+realdt+cumday+realgroup:team1, data=recapFCcom, family="binomial") m35233 <- glm(disagno~realgroup+team1+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m31,m32,m33,m34,m35,m36,m37,m351,m352,m353,m354,m355,m3521,m3522,m3523,m3524,m35231,m35232,m3523, rank="AICc") m352321 <- glm(disagno~realgroup+team1+realdt+realgroup:team1, data=recapFCcom, family="binomial") m352322 <- glm(disagno~realgroup+team1+cumday+realgroup:team1, data=recapFCcom, family="binomial") m352323 <- glm(disagno~realgroup+team1+realgroup:team1, data=recapFCcom, family="binomial") m352324 <- glm(disagno~realgroup+team1+cumday, data=recapFCcom, family="binomial") m352325 <- glm(disagno~realgroup+team1+realdt, data=recapFCcom, family="binomial") m352326 <- glm(disagno~realgroup+realdt+cumday, data=recapFCcom, family="binomial") m352327 <- glm(disagno~team1+realdt+cumday, data=recapFCcom, family="binomial") m35232null <- glm(disagno~1, data=recapFCcom, family="binomial") model.sel(m31,m32,m33,m34,m35,m36,m37,m351,m352,m353,m354, m355,m3521,m3522,m3523,m3524,m35231,m35232,m3523, m352321,m352322,m352323,m352324,m352325,m352326,m352327,m35232null, rank="AICc") ################################### ####### Refitting with full cases gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting only marine mammal sightings in the transect recapFCcom1 <- subset(gpbin, realdt %in% c("Dolphin", "Dugong") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom <- completeFun(recapFCcom1, c("realgroup", "team1")) m352323 <- glm(disagno~realgroup+team1+realgroup:team1, data=recapFCcom, family="binomial") summary(m352323) drop1(m352323, test="Chi") anova(m352323, test="Chi") ## getting predictions for final model newdat <- expand.grid(team1=levels(recapFCcom$team1), realgroup=seq(2,25,0.5)) preds <- predict(m352323, newdata=newdat, se.fit=T, interval="confidence", level=0.95, type="response") newdat$fit <- preds$fit newdat$sefit<- preds$se.fit newdat levels(newdat$team) ndep <- subset(newdat, team1=="east_Port") ndes <- subset(newdat, team1=="east_Stbd") ndwp <- subset(newdat, team1=="west_Port") ndws <- subset(newdat, team1=="west_Stbd") ### Making postscript plot pdf("group-size.pdf", height = 3, width = 7.48) par(mfrow=c(1,3), mar=c(4,4,2,1)) plot(ndep$fit~ndep$realgroup, pch=NA, xlim=c(2,25), ylim=c(0,1), xaxt='n', frame.plot=F, xlab="Group Size", ylab="Probability of Group Size Mismatch") axis(1, at=c(2,5,10,15,20,25)) lines(smooth.spline(ndep$realgroup, ndep$fit, df=50), col="red") polygon(c(ndep$realgroup, rev(ndep$realgroup)), c(ndep$fit-ndep$sefit , rev(ndep$fit+ndep$sefit)),col = rgb(1,0,0, 0.25), border = NA) lines(smooth.spline(ndes$realgroup, ndes$fit, df=50), col="coral") ndes$sefit2 <- ifelse(ndes$fit+ndes$sefit>1,1,ndes$fit+ndes$sefit) polygon(c(ndes$realgroup, rev(ndes$realgroup)), c(ndes$fit-ndes$sefit , rev(ndes$sefit2)),col = rgb(1,127/255,80/255, 0.25), border = NA) lines(smooth.spline(ndwp$realgroup, ndwp$fit, df=50), col="blue") polygon(c(ndwp$realgroup, rev(ndwp$realgroup)), c(ndwp$fit-ndwp$sefit , rev(ndwp$fit+ndwp$sefit)),col = rgb(0,0,1, 0.25), border = NA) lines(smooth.spline(ndws$realgroup, ndws$fit, df=50), col="cornflowerblue") polygon(c(ndws$realgroup, rev(ndws$realgroup)), c(ndws$fit-ndws$sefit , rev(ndws$fit+ndws$sefit)),col = rgb(100/255,149/255,237/255, 0.25), border = NA) legend(17.5,0.45, lty=1, col=c("red", "coral", "blue", "cornflowerblue"), c("East Port", "East Stbd", "West Port","West Stbd"), ncol=1, bty="n", cex=0.85, xpd=NA) points(jitter(recapFCcom$disagno, factor=0.2)~recapFCcom$realgroup, col=rgb(0,0,0)) text(-4,1.1,labels=c("A."), xpd=NA, cex=2) ############################################# ######################### #### Dolphins Seperately gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting only marine mammal sightings in the transect levels(gpbin$transect) levels(gpbin$transect_position) recapFCcom1 <- subset(gpbin, realdt %in% c("Dolphin") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom1 <- droplevels(recapFCcom1) levels(recapFCcom1$transect) levels(recapFCcom1$transect_position) str(recapFCcom1) ## Got all dolphin and dugong sightings, now getting complete cases for analysis completeFun <- function(data, desiredCols) { completeVec <- complete.cases(data[, desiredCols]) return(data[completeVec, ]) } recapFCcom <- completeFun(recapFCcom1, c("realgroup", "team1", "transect_position", "beaufort", "turbidity", "glare")) str(recapFCcom) count(recapFCcom, c("team1")) table(recapFCcom$team1, recapFCcom$turbidity) count(recapFCcom, c("team1","turbidity", "disagno" )) ### only one Turbidity = 1 sighting, initially dropping this to examine turbidity effects ## Dropping single T=1 sightings recapFCcom <- subset(recapFCcom, turbidity !=1) str(recapFCcom) m1 <- glm(disagno~realgroup*team1+transect_position*as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday, data=recapFCcom, family="binomial") summary(m1) m2 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m3 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") library(MuMIn) model.sel(m1,m2,m3, rank="AICc") m31 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday, data=recapFCcom, family="binomial") m32 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+realgroup:team1, data=recapFCcom, family="binomial") m33 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m34 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m35 <- glm(disagno~realgroup+team1+transect_position+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m36 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m37 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m38 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m1,m2,m3, m31,m32,m33,m34,m35,m36,m37, rank="AICc") m361 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+realgroup:team1, data=recapFCcom, family="binomial") m362 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m363 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m364 <- glm(disagno~realgroup+team1+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m365 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m1,m2,m3, m31,m32,m33,m34,m35,m36,m37,m361,m362,m363,m364,m365, rank="AICc") m3621 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(turbidity)+realgroup:team1, data=recapFCcom, family="binomial") m3622 <- glm(disagno~realgroup+team1+as.factor(glare)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m3623 <- glm(disagno~realgroup+team1+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") m3624 <- glm(disagno~realgroup+team1+as.factor(glare)+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m1,m2,m3, m31,m32,m33,m34,m35,m36,m37,m361,m362,m363,m364,m365,m3621,m3622,m3623,m3624, rank="AICc") m36231 <- glm(disagno~realgroup+team1+as.factor(turbidity)+realgroup:team1, data=recapFCcom, family="binomial") m36232 <- glm(disagno~realgroup+team1+cumday+realgroup:team1, data=recapFCcom, family="binomial") m36233 <- glm(disagno~realgroup+team1+as.factor(turbidity)+cumday+realgroup:team1, data=recapFCcom, family="binomial") model.sel(m1,m2,m3, m31,m32,m33,m34,m35,m36,m37, m361,m362,m363,m364,m365,m3621,m3622,m3623,m3624, m36231,m36232,m36233, rank="AICc") m362321 <- glm(disagno~realgroup+team1+cumday, data=recapFCcom, family="binomial") m362322 <- glm(disagno~realgroup+team1+realgroup:team1, data=recapFCcom, family="binomial") m362323 <- glm(disagno~realgroup+team1+cumday+realgroup:team1, data=recapFCcom, family="binomial") m362324 <- glm(disagno~realgroup+team1+realgroup:team1, data=recapFCcom, family="binomial") m362325 <- glm(disagno~realgroup+team1, data=recapFCcom, family="binomial") m362326 <- glm(disagno~realgroup, data=recapFCcom, family="binomial") m362327 <- glm(disagno~team1, data=recapFCcom, family="binomial") m36232null <- glm(disagno~1, data=recapFCcom, family="binomial") model.sel(m1,m2,m3, m31,m32,m33,m34,m35,m36,m37,m361,m362,m363,m364,m365,m3621,m3622,m3623, m36231,m36232 ,m362321,m362322,m362325,m362326,m362327,m36232null, rank="AICc") ################################# ## Getting complete cases gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting only marine mammal sightings in the transect recapFCcom1 <- subset(gpbin, realdt %in% c("Dolphin") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom <- completeFun(recapFCcom1, c("realgroup", "team1")) m362322 <- glm(disagno~realgroup+team1+realgroup:team1, data=recapFCcom, family="binomial") drop1(m362322, test="Chi") summary(m362322) ## getting predictions for final model newdat <- expand.grid(team1=levels(recapFCcom$team1), realgroup=seq(2,25,0.5)) preds <- predict(m36413, newdata=newdat, se.fit=T, interval="confidence", level=0.95, type="response") newdat$fit <- preds$fit newdat$sefit<- preds$se.fit levels(newdat$team) ndep <- subset(newdat, team1=="east_Port") ndes <- subset(newdat, team1=="east_Stbd") ndwp <- subset(newdat, team1=="west_Port") ndws <- subset(newdat, team1=="west_Stbd") ### Making postscript plot #pdf("group-size.pdf", height = 3.7, width = 3.54) par(mar=c(4,4,2,1)) plot(ndep$fit~ndep$realgroup, pch=NA, xlim=c(2,25), ylim=c(0,1), xaxt='n', frame.plot=F, xlab="Group Size", ylab=NA) axis(1, at=c(2,5,10,15,20,25)) lines(smooth.spline(ndep$realgroup, ndep$fit, df=50), col="red") polygon(c(ndep$realgroup, rev(ndep$realgroup)), c(ndep$fit-ndep$sefit , rev(ndep$fit+ndep$sefit)),col = rgb(1,0,0, 0.25), border = NA) lines(smooth.spline(ndes$realgroup, ndes$fit, df=50), col="coral") ndes$sefit2 <- ifelse(ndes$fit+ndes$sefit>1,1,ndes$fit+ndes$sefit) polygon(c(ndes$realgroup, rev(ndes$realgroup)), c(ndes$fit-ndes$sefit , rev(ndes$sefit2)),col = rgb(1,127/255,80/255, 0.25), border = NA) lines(smooth.spline(ndwp$realgroup, ndwp$fit, df=50), col="blue") polygon(c(ndwp$realgroup, rev(ndwp$realgroup)), c(ndwp$fit-ndwp$sefit , rev(ndwp$fit+ndwp$sefit)),col = rgb(0,0,1, 0.25), border = NA) lines(smooth.spline(ndws$realgroup, ndws$fit, df=50), col="cornflowerblue") polygon(c(ndws$realgroup, rev(ndws$realgroup)), c(ndws$fit-ndws$sefit , rev(ndws$fit+ndws$sefit)),col = rgb(100/255,149/255,237/255, 0.25), border = NA) legend(17.5,0.45, lty=1, col=c("red", "coral", "blue", "cornflowerblue"), c("East Port", "East Stbd", "West Port","West Stbd"), ncol=1, bty="n", cex=0.85, xpd=NA) points(jitter(recapFCcom$disagno, factor=0.2)~recapFCcom$realgroup, col=rgb(0,0,0)) text(-4,1.1,labels=c("B."), xpd=NA, cex=2) ################################################### ################################ ## Examining dugongs seperately. gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting only marine mammal sightings in the transect levels(gpbin$transect) levels(gpbin$transect_position) levels(gpbin$realdt) recapFCcom1 <- subset(gpbin, realdt %in% c("Dugong") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom1 <- droplevels(recapFCcom1) levels(recapFCcom1$transect) levels(recapFCcom1$transect_position) str(recapFCcom1) count(recapFCcom1, c( "glare","disagno" )) count(recapFCcom1, c( "beaufort","disagno" )) ## Got all dolphin and dugong sightings, now getting complete cases for analysis completeFun <- function(data, desiredCols) { completeVec <- complete.cases(data[, desiredCols]) return(data[completeVec, ]) } ## There is one case each at turbidity 4 and at glare 3 - dropping both these to examine if ## there are any glare or turbidity effects recapFCcom <- completeFun(recapFCcom1, c("realgroup", "beaufort", "glare")) recapFCcom <- subset(recapFCcom, glare !=3) recapFCcom <- subset(recapFCcom, beaufort !=4) str(recapFCcom) count(recapFCcom, c( "beaufort","disagno" )) count(recapFCcom, c( "realgroup")) count(recapFCcom, c("transect_position", "glare" )) m1 <- glm(disagno~realgroup+as.factor(beaufort)+as.factor(glare), data=recapFCcom, family="binomial") summary(m1) m2 <- glm(disagno~realgroup+as.factor(beaufort), data=recapFCcom, family="binomial") m3 <- glm(disagno~realgroup+as.factor(glare), data=recapFCcom, family="binomial") m4 <- glm(disagno~as.factor(beaufort)+as.factor(glare), data=recapFCcom, family="binomial") m5 <- glm(disagno~realgroup, data=recapFCcom, family="binomial") m6 <- glm(disagno~as.factor(beaufort), data=recapFCcom, family="binomial") m7 <- glm(disagno~as.factor(glare), data=recapFCcom, family="binomial") library(MuMIn) model.sel(m1,m2,m3,m4,m5,m6,m7, rank="AICc") drop1(m5, test="Chi") summary(m5) gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting o complete cases for final model gpbin <- read.csv("amGRPbinaryNEW.csv") recapFCcom1 <- subset(gpbin, realdt %in% c("Dugong") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom <- completeFun(recapFCcom1, c("realgroup")) str(recapFCcom) m5 <- glm(disagno~realgroup, data=recapFCcom, family="binomial") drop1(m5, test="Chi") summary(m5) ## getting predictions for final model summary(recapFCcom$realgroup) newdat1 <- expand.grid(realgroup=seq(2,7,1)) preds <- predict(m5, newdata=newdat1, se.fit=T, interval="confidence", level=0.95, type="response") newdat1$fit <- preds$fit newdat1$sefit<- preds$se.fit newdat1 plot(newdat1$fit~newdat1$realgroup, pch=NA, xlim=c(2,7), ylim=c(0,1), xaxt='n', frame.plot=F, xlab="Group Size", ylab=NA) axis(1, at=c(2:7)) lines(newdat1$realgroup, newdat1$fit, col="black") polygon(c(newdat1$realgroup, rev(newdat1$realgroup)), c(newdat1$fit-newdat1$sefit , rev(newdat1$fit+newdat1$sefit)),col = rgb(0,0,0, 0.25), border = NA) points(jitter(recapFCcom$disagno, factor=0.2)~recapFCcom$realgroup, col=rgb(0,0,0)) text(0.6,1.1,labels=c("C."), xpd=NA, cex=2) dev.off() ## No relation in dugongs between team or group size and probability of group disagreement, but data are too sparse (74 recaptures in transect, 8 grp misspecifications) ## there are 10% of sightings where group size is disagreed, 6/8 of which were a disagreement between whether there were 2 or 1 animal, ## of these 3/6 were where one observer missed a calf and called the group size 1. ######################### #### Turtles Seperately gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) ## Getting only marine mammal sightings in the transect levels(gpbin$transect) levels(gpbin$transect_position) levels(gpbin$) recapFCcom1 <- subset(gpbin, realdt %in% c("Turtle") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom1 <- droplevels(recapFCcom1) levels(recapFCcom1$transect) levels(recapFCcom1$transect_position) str(recapFCcom1) ## Got all dolphin and dugong sightings, now getting complete cases for analysis completeFun <- function(data, desiredCols) { completeVec <- complete.cases(data[, desiredCols]) return(data[completeVec, ]) } recapFCcom <- completeFun(recapFCcom1, c("realgroup", "team1", "transect_position", "beaufort", "turbidity", "glare")) str(recapFCcom) table(recapFCcom$team1, recapFCcom$turbidity) table(recapFCcom$team1, recapFCcom$disagno) library(reshape) library(reshape2) count(recapFCcom, c("team1", "realgroup" )) ### only one Turbidity = 1 sighting, initially dropping this to examine turbidity effects ### Also West STBD had no dugong sightings with a group size disagreement, so can't include team*datatype interaction ## Because all sightings of turtle groups >4 are from one team - East Starboard ## - except two sightings from other teams (one from EP (grpsize 1), the other from WP (grpsize=7)), sightings ## >4 were excluded from this analysis because strong effect of ES with group size. recapFCcom <- subset(recapFCcom, realgroup <5) count(recapFCcom, c("realgroup")) m1 <- glm(disagno~realgroup+team1+transect_position*as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday, data=recapFCcom, family="binomial") summary(m1) step(m1) ## reducing from the step results m2 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m3 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+realgroup:team1, data=recapFCcom, family="binomial") library(MuMIn) model.sel(m1,m2,m3, rank="AICc") m21 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday , data=recapFCcom, family="binomial") m22 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m23 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m24 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m25 <- glm(disagno~realgroup+team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m26 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m27 <- glm(disagno~team1+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") model.sel(m2,m3, m21,m22,m23,m24,m25,m26,m27, rank="AICc") m261 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday , data=recapFCcom, family="binomial") m262 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m263 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(turbidity)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m264 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m265 <- glm(disagno~transect_position+as.factor(glare)+as.factor(turbidity)+as.factor(beaufort)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") model.sel(m2,m3, m21,m22,m23,m24,m25,m26,m27, m261,m262,m263,m264,m265, rank="AICc") m2631 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(turbidity)+cumday , data=recapFCcom, family="binomial") m2632 <- glm(disagno~realgroup+transect_position+as.factor(glare)+as.factor(turbidity)+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m2633 <- glm(disagno~realgroup+transect_position+as.factor(glare)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m2634 <- glm(disagno~transect_position+as.factor(glare)+as.factor(turbidity)+cumday+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") model.sel(m2,m3, m21,m22,m23,m24,m25,m26,m27, m261,m262,m263,m264,m265,m2631,m2632,m2633,m2634, rank="AICc") m26321 <- glm(disagno~realgroup+transect_position+as.factor(glare)+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m26322 <- glm(disagno~transect_position+as.factor(glare)+as.factor(turbidity)+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m263211 <- glm(disagno~transect_position+as.factor(glare)+transect_position:as.factor(glare) , data=recapFCcom, family="binomial") m263212 <- glm(disagno~realgroup+transect_position , data=recapFCcom, family="binomial") m263213 <- glm(disagno~realgroup+as.factor(glare) , data=recapFCcom, family="binomial") m263214 <- glm(disagno~transect_position+as.factor(glare) , data=recapFCcom, family="binomial") m263215 <- glm(disagno~realgroup , data=recapFCcom, family="binomial") m263216 <- glm(disagno~transect_position , data=recapFCcom, family="binomial") m263217 <- glm(disagno~as.factor(glare) , data=recapFCcom, family="binomial") mnull <- glm(disagno~as.factor(glare) , data=recapFCcom, family="binomial") model.sel(m2,m3, m21,m22,m23,m24,m25,m26,m27, m261,m262,m263,m264,m265,m2631,m2632,m2633,m2634,m26321,m26322, m263211, m263212, m263213, m263214, m263215, m263216, m263217, mnull, rank="AICc") drop1(m26321, test="Chi") ################################# ## Getting complete cases gpbin <- read.csv("amGRPbinaryNEW.csv") str(gpbin) recapFCcom1 <- subset(gpbin, realdt %in% c("Turtle") & transect_position %in% c("low","medium","high","veryhigh") & transect !="off transect" & transect !="Off transect"& transect !="Off Transect") recapFCcom1 <- droplevels(recapFCcom1) recapFCcom <- completeFun(recapFCcom1, c("realgroup")) str(recapFCcom) ### only one Turbidity = 1 sighting, initially dropping this to examine turbidity effects ### Also West STBD had no dugong sightings with a group size disagreement, so can't include team*datatype interaction ## Because all sightings of turtle groups >4 are from one team - East Starboard ## - except two sightings from other teams (one from EP (grpsize 1), the other from WP (grpsize=7)), sightings ## >4 were excluded from this analysis because strong effect of ES with group size. recapFCcom <- subset(recapFCcom, realgroup <5) count(recapFCcom, c("realgroup", "disagno" )) m263215 <- glm(disagno~realgroup , data=recapFCcom, family="binomial") mnull <- glm(disagno~realgroup , data=recapFCcom, family="binomial") drop1(m263215, test="Chi") summary(m263215) ## getting predictions for final model newdat <- expand.grid(realgroup=seq(1,4,1)) preds <- predict(m263215, newdata=newdat, se.fit=T, interval="confidence", level=0.95, type="response") newdat$fit <- preds$fit newdat$sefit<- preds$se.fit newdat levels(newdat$team) ndep <- subset(newdat, team1=="east_Port") ndes <- subset(newdat, team1=="east_Stbd") ndwp <- subset(newdat, team1=="west_Port") ndws <- subset(newdat, team1=="west_Stbd")