Hello RMark Users,
I am attempting to run my first survival analysis using RMark, but am faced with odd results that are quite confusing. My dataset includes a capture history of 14 encounters for each individual as well as 4 covariates (season, site, sex, and age) that I wish to use to inform survival (S) and recapture probability (p) of 238 individuals I've tracked over the past 3 years. To do this, I ran the following code:
Encounter <- read.csv("EncounterDataNewSeason.csv",
colClasses=c("character", "factor", "factor","factor",
"factor","factor","factor"),
header = TRUE, sep = ",")
#Process Data and filter groups
Surv <- process.data(Encounter, model = "Burnham", groups = c("Sex", "Site", "Season", "Age"))
#Create Design Data
Surv.ddl <- make.design.data(Surv)
#Create function to run models specified within function
run.surv <- function() {
S.sex <- list(formula=~Sex)
S.site <- list(formula=~Site)
S.season <- list(formula=~Season)
S.age <- list(formula=~Age)
s.dot <- list(formula=~1)
p.sex <- list(formula=~Sex)
p.site <- list(formula=~Site)
p.season <- list(formula=~Season)
p.age <- list(formula=~Age)
p.dot <- list(formula=~1)
F.dot <- list(formula=~1)
r.dot <- list(formula=~1)
# Create model list
model.list = create.model.list("Burnham")
Surv.results = mark.wrapper(model.list, data = Surv,
ddl = Surv.ddl, output = FALSE, invisible = TRUE, threads = 2)
# Return model table and list of models
return(Surv.results)
}
#View results
Surv.results = run.surv()
Surv.results
This code seems to work, and I am able to compare AIC values for different models representing each of the combinations of covariates (ex. S(sex), p(site); S(sex), p(age); etc.) The issue arises once I run just my top model (which I specify as follows):
top <- Surv.results$S.season.p.site.r.dot.F.dot
top$results
This provides output, but the $beta estimates are very large (not 0 to 1 as I would expect for survival estimates), $beta confidence intervals are matching (lcl = ucl on many occasions), and $real estimates appear to include all covariates in the model even though the model should only include 2 of the covariates (ex. model of S(season), p(site) has $real estimates for models with sex, site, season, and age all incorporated).
As this is my first time working with RMark, I am unsure if output like this is normal or if the results are nonsensical like I suspect. Any guidance in correctly interpreting the output and/or pointing out issues in the code would be greatly appreciated! Thank you so much in advance for your time and thoughts on this!