###Ambystoma CJS Models Goodness-of-fit Testing ## Arianne Messerman, Raymond Semlitsch, and Manuel Leal ## February 2019 # Load data for species-specific frequency amb <- read.table("Ambystoma_PenCode.txt",header=T,colClasses="character") head(amb) amb$spp<-as.factor(amb$spp) amb$tag<-as.factor(amb$tag) amb$block<-as.factor(amb$block) amb$pen.join<-as.factor(amb$pen.join) amb$pen.ID<-as.factor(amb$pen.ID) str(amb) ##Subset by species aman<-subset(amb, amb$spp=="AMAN") str(aman) amma<-subset(amb, amb$spp=="AMMA") str(amma) amte<-subset(amb, amb$spp=="AMTE") str(amte) # Create Capture History (CH) matrix from entire data set CH<-t(array(as.numeric(unlist(strsplit(amb$history,""))),dim=c(nchar(amb$history[1]),length(amb$tag)))) CH[1:11,] amb$history[10] # Check to make sure first non-zero is accurate... looks good # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) # Identifies first occasion of marking (f). Character ! stands for "is not" f <- apply(CH, 1, get.first) #------------------------------------------------------------------ # Create Capture History (CH) matrix from AMAN only data set CHaman<-t(array(as.numeric(unlist(strsplit(aman$history,""))),dim=c(nchar(aman$history[1]),length(aman$tag)))) CHaman[1:11,] aman$history[10] # Check to make sure first non-zero is accurate... looks good # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) # Identifies first occasion of marking (f). Character ! stands for "is not" faman <- apply(CHaman, 1, get.first) #------------------------------------------------------------------ # Create Capture History (CH) matrix from AMMA only data set CHamma<-t(array(as.numeric(unlist(strsplit(amma$history,""))),dim=c(nchar(amma$history[1]),length(amma$tag)))) CHamma[1:11,] amma$history[10] # Check to make sure first non-zero is accurate... looks good # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) # Identifies first occasion of marking (f). Character ! stands for "is not" famma <- apply(CHamma, 1, get.first) #------------------------------------------------------------------ # Create Capture History (CH) matrix from AMMA only data set CHamte<-t(array(as.numeric(unlist(strsplit(amte$history,""))),dim=c(nchar(amte$history[1]),length(amte$tag)))) CHamte[1:11,] amte$history[10] # Check to make sure first non-zero is accurate... looks good # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) # Identifies first occasion of marking (f). Character ! stands for "is not" famte <- apply(CHamte, 1, get.first) #------------------------------------------------------------------ ##Create vector of unequal study period intervals intervals<-read.csv("interval_primary.csv",header=T) str(intervals) int<-c(intervals$interval) length(int) is.numeric(int) library(R2ucare) #Frequency by species amb.freq <- read.table("Amb_freq.txt",header=T,colClasses="character") a.hist<-matrix(as.numeric(unlist(strsplit(as.character(amb.freq$history),""))),nrow=length(amb.freq$history),byrow=T) a.freq<-as.numeric(amb.freq$freq) a.sp<-amb.freq$spp sum(a.freq) mask = (a.sp == 'AMAN') aman.hist = a.hist[mask,] aman.freq = a.freq[mask] mask = (a.sp == 'AMMA') amma.hist = a.hist[mask,] amma.freq = a.freq[mask] mask = (a.sp == 'AMTE') amte.hist = a.hist[mask,] amte.freq = a.freq[mask] #Frequency over all species amb1 <- read.table("Amb_freq_nospp.txt" ,header=T,colClasses="character") head(amb1) a1.hist<-matrix(as.numeric(unlist(strsplit(as.character(amb1$history),""))),nrow=length(amb1$history),byrow=T) a1.freq<-as.numeric(amb1$freq) sum(a1.freq) overall_CJS(a1.hist, a1.freq)#p=0.001, X^2=88.552, df=51 overall_CJS(aman.hist, aman.freq) #p=0.122, X^2=40.312, df=31 overall_CJS(amma.hist, amma.freq) #p=0.171, X^2=38.319, df=31 overall_CJS(amte.hist, amte.freq) #p=0.471, X^2=32.917, df=33 ## None significant, implying adequate fit of model to data #Calculate C-hat, i.e. amount of overdispersion c.amb<-88.552/51 #1.736, #c-hat greater than 1, contribution to the QAICc value from the model likelihood will #decline, and thus the relative penalty for a given number of parameters K will increase c.aman<-40.312/31 #1.300 c.amma<-38.319/31 #1.236 c.amte<-32.917/33 #0.998 test3sr_aman = test3sr(aman.hist, aman.freq) #Not useful because no difference in first occasion marking test3sm_aman = test3sm(aman.hist, aman.freq) #Not useful because no difference in first occasion marking test2ct_aman = test2ct(aman.hist, aman.freq) #Tests for nonrandom recapture effects test2cl_aman = test2cl(aman.hist, aman.freq) #Tests for nonrandom recapture effects # display results: test3sr_aman test3sm_aman test2ct_aman test2cl_aman #None significant test2ct_amma = test2ct(amma.hist, amma.freq) test2cl_amma = test2cl(amma.hist, amma.freq) # display results: test2ct_amma #Short-term recapture effect test2cl_amma test2ct_amte = test2ct(amte.hist, amte.freq) test2cl_amte = test2cl(amte.hist, amte.freq) # display results: test2ct_amte test2cl_amte #None Significant test2ct_amb = test2ct(a1.hist, a1.freq) test2cl_amb = test2cl(a1.hist, a1.freq) # display results: test2ct_amb #Short-term recapture effect test2cl_amb