zc_fit.R

### Load packages ###
# Must have 'marked' v >= 1.1.8
library(marked)
# The splines package is only necessary for fitting b-spline curves used in the paper
# It is not required for the multivate models in the marked package
library(splines)


### Import and process data ###

chdata=read.table("mvms.dat",header=T,sep="\t",stringsAsFactors=FALSE)
# create weight anomaly by subtracting sex-specific weight
chdata$weight[chdata$sex=="F"]=chdata$weight[chdata$sex=="F"]-mean(chdata$weight[chdata$sex=="F"])
chdata$weight[chdata$sex=="M"]=chdata$weight[chdata$sex=="M"]-mean(chdata$weight[chdata$sex=="M"])
# Process data for multivariate models in marked
dp=process.data(chdata,model="mvmscjs",strata.labels=list(area=c("A","S"),ltag=c("+","-","u"),rtag=c("+","-","u")))


### Make design data

ddl=make.design.data(dp)
# Create pup variable for Phi
ddl$Phi$pup=ifelse(ddl$Phi$Age==0,1,0)
ddl$Phi$sex=factor(ddl$Phi$sex)


### Fix parameters for references cells and structural zeros ###

# Detection model
# Set final year (2014)  p=0 (no resight data) for ANI
ddl$p$fix = ifelse(ddl$p$Time==17 & ddl$p$area=="A", 0, ddl$p$fix)

# Delta model
# create indicator variables for 'unknown' tag observations
ddl$delta$obs.ltag.u = ifelse(ddl$delta$obs.ltag=="u", 1, 0)
ddl$delta$obs.rtag.u = ifelse(ddl$delta$obs.rtag=="u", 1, 0)

# Psi model
# Set Psi to 0 for cases which are not possible - missing tag to having tag
ddl$Psi$fix[as.character(ddl$Psi$ltag)=="-"&as.character(ddl$Psi$toltag)=="+"]=0
ddl$Psi$fix[as.character(ddl$Psi$rtag)=="-"&as.character(ddl$Psi$tortag)=="+"]=0
# Create indicator variables for transitioning between states
ddl$Psi$AtoS=ifelse(ddl$Psi$area=="A"&ddl$Psi$toarea=="S",1,0)  # ANI to SMI movement
ddl$Psi$StoA=ifelse(ddl$Psi$area=="S"&ddl$Psi$toarea=="A",1,0)  # SMI to ANI movement
ddl$Psi$lpm=ifelse(ddl$Psi$ltag=="+"&ddl$Psi$toltag=="-",1,0)   # Losing left tag
ddl$Psi$rpm=ifelse(ddl$Psi$rtag=="+"&ddl$Psi$tortag=="-",1,0)   # Losing right tag
ddl$Psi$sex=factor(ddl$Psi$sex)

# Fit model 
mod=crm(dp,ddl,model.parameters=list(Psi=list(formula=~-1+ AtoS:sex + AtoS:sex:bs(Age) + StoA:sex + StoA:sex:bs(Age) + I(lpm+rpm) +I(lpm+rpm):Age + lpm:rpm),
                                     p=list(formula=~time*area),
                                     delta=list(formula= ~ -1 + obs.ltag.u + obs.rtag.u + obs.ltag.u:obs.rtag.u),
                                     Phi=list(formula=~sex*bs(Age)+pup:weight+area)), hessian=TRUE)





save(list=ls(), file="mvms_zc_results_annotated.RData")

zc_fit_summary.R

library(marked)
library(splines)
library(mvtnorm)
library(dplyr)
library(ggplot2)
theme_set(theme_bw(base_size = 12))
load("mvms_zc_results.RData")

# Obtain variance-covariance matrix and coefficient estimates
vcv = mod$results$beta.vcv
vcv_names = rownames(mod$results$beta.vcv)

## Delta
mod$model.parameters$delta$formula
delta = data.frame(expand.grid(obs.ltag.u=c(0,1), obs.rtag.u=c(0,1)))
delta_dm = model.matrix(mod$model.parameters$delta$formula,data=delta)
delta_coef = mod$results$beta$delta
delta_vcv = vcv[grep("delta", vcv_names),grep("delta", vcv_names)]
delta_rep = rmvnorm(10000, delta_coef, delta_vcv)
delta$estimate = exp(delta_dm%*%delta_coef)/sum(exp(delta_dm%*%delta_coef))
pred_rep=apply(delta_rep, 1, FUN=function(b,D){exp(D%*%b)/sum(exp(D%*%b))}, D=delta_dm)
delta = cbind(delta, t(apply(pred_rep, 1, quantile, prob=c(0.025, 0.975))))
names(delta)[4:5] = c("lower","upper")
delta$obs = with(delta, paste(ifelse(obs.ltag.u==1,"Unk","Obs"),ifelse(obs.rtag.u==1,"Unk","Obs"),sep='/'))
p_delta = ggplot(delta) + geom_bar(aes(x=obs, y=estimate), stat="identity", fill=gray(0.5)) + 
  geom_errorbar(aes(x=obs, ymin=lower, ymax=upper), width=0.25) +
  ylab(expression(delta)) + xlab("Observed tag status (Left/Right)") 
ggsave("p_delta.pdf", p_delta, width=5, height = 3)
# odds of missing second tag given first was not determined
delta_dep = exp(delta_coef[3])
delta_dep_ci = exp(delta_coef[3] +c(-1,1)*1.96*sqrt(delta_vcv[3,3]))


## Phi
mod$model.parameters$Phi$formula
#pup
Phi_pup = rbind(
  data.frame(
    sex="F",
    weight=seq(min(ddl$Psi$weight[ddl$Psi$sex=="F"]), max(ddl$Psi$weight[ddl$Psi$sex=="F"]),0.1)
  ),
  data.frame(
    sex="M",
    weight=seq(min(ddl$Psi$weight), max(ddl$Psi$weight),0.1)
  )
)
Phi_pup_dm = model.matrix(~sex  + weight, data=Phi_pup)
Phi_coef = mod$results$beta$Phi
Phi_vcv = vcv[grep("Phi", vcv_names),grep("Phi", vcv_names)]
Phi_rep = rmvnorm(10000, Phi_coef, Phi_vcv)
pred_rep=apply(Phi_rep[,c(1,2,10)], 1, FUN=function(b,D){plogis(D%*%b)}, D=Phi_pup_dm)
xxx = t(apply(pred_rep, 1, quantile, prob=c(0.025, 0.975)))
colnames(xxx) = c("lower","upper")
Phi_pup = cbind(Phi_pup,xxx)
Phi_pup$estimate = plogis(Phi_pup_dm%*%Phi_coef[c(1,2,10)])
p_pup_surv = ggplot(data=Phi_pup) + geom_line(aes(x=weight, y=estimate, color=sex)) +
  geom_ribbon(aes(ymin=lower, ymax=upper, x=weight, fill=sex), alpha=0.33) + 
  ylab(expression(phi[0])) + xlab("Mass anomaly (kg)")
ggsave("p_pup_surv.pdf", p_pup_surv, width=5, height=3)

#adult
Phi_ad = data.frame(
  expand.grid(
    sex=c("F","M"),
    Age=c(1:17),
    area=c("A","S"),
    pup=0,
    weight=0
  )
)
Phi_ad_dm = model.matrix(mod$model.parameters$Phi$formula, Phi_ad)  
Phi_ad$estimate = plogis(Phi_ad_dm%*%Phi_coef)
xxx = Phi_rep %>% apply(., 1, FUN=function(b,D){plogis(D%*%b)}, D=Phi_ad_dm) %>% apply(., 1, quantile, prob=c(0.025, 0.975)) %>% t(.)
colnames(xxx)=c("lower","upper")
Phi_ad = cbind(Phi_ad,xxx)
p_ad_surv = ggplot(data=Phi_ad) + geom_line(aes(x=Age, y=estimate, color=area)) +
  geom_ribbon(aes(ymin=lower, ymax=upper, x=Age, fill=area), alpha=0.33) + 
  facet_wrap(~sex) +
  ylab(expression(phi[Age])) + xlab("Age (yr)")
ggsave("p_ad_surv.pdf", p_ad_surv, width=5, height=3)


## p
mod$model.parameters$p$formula
p=predict(mod,ddl=ddl,parameter="p",unique=TRUE,select=c("area","Time"))
p$year = p$Time+1997
p_dm = model.matrix(mod$model.parameters$p$formula,data=p)
p_coef = mod$results$beta$p
#p_dm = p_dm[,colnames(p_dm) %in% names(p_coef)]
p_vcv = vcv[grep("p.", vcv_names,fixed=TRUE),grep("p.", vcv_names,fixed=TRUE)]
p_rep = rmvnorm(10000, p_coef, p_vcv, method="svd")
pred_rep=apply(p_rep, 1, FUN=function(b,D){plogis(D%*%b)}, D=p_dm)
xxx = t(apply(pred_rep, 1, quantile, prob=c(0.025, 0.975)))
colnames(xxx) = c("lower","upper")
p = cbind(p,xxx)
p=p[p$area!="A" | p$year!=2014,]
p_det = ggplot(data=p) + geom_line(aes(x=year, y=estimate, color=area)) +
  geom_ribbon(aes(ymin=lower, ymax=upper, x=year, fill=area), alpha=0.33) + 
  ylab(expression(p[Area])) + xlab("Year")
ggsave("p_det.pdf", p_det, width=5, height=3)


## Psi
mod$model.parameters$Psi$formula
Psi = data.frame(
  expand.grid(
    AtoS = c(0,1),
    StoA = c(0,1),
    Age=c(0:17),
    sex=c("F","M"),
    lpm=0,
    rpm=0
  )
)
Psi = with(Psi,Psi[StoA + AtoS==1,])
Psi$travel = ifelse(Psi$StoA==1, "S to A", "A to S")
Psi_dm = model.matrix(mod$model.parameters$Psi$formula,data=Psi)
Psi_coef = mod$results$beta$Psi
Psi_vcv = vcv[grep("Psi", vcv_names),grep("Psi", vcv_names)]
Psi_rep = rmvnorm(10000, Psi_coef, Psi_vcv, method="svd")
Psi$estimate = plogis(Psi_dm%*%Psi_coef)
xxx = apply(Psi_rep, 1, FUN=function(b,D){plogis(D%*%b)}, D=Psi_dm) %>% apply(.,1,quantile, prob=c(0.025,0.975)) %>% t(.)
colnames(xxx)=c("lower","upper")
Psi = cbind(Psi,xxx)
# movement between areas
p_move=ggplot(data=Psi) + geom_line(aes(x=Age, y=estimate, color=sex)) + geom_ribbon(aes(x=Age, ymin=lower, ymax=upper, fill=sex), alpha=0.25) + 
  facet_wrap(~travel, scales="free_y", ncol=1) + xlab("Age (yrs)") + ylab("Prob. of movement between areas\n")
ggsave("p_move.pdf",p_move, width=5, height=5)
# tag loss
tl_idx = c(grep("lpm",names(Psi_coef)))
tl_coef = Psi_coef[tl_idx]
tl_vcv = Psi_vcv[tl_idx,tl_idx]
tl_rep = rmvnorm(10000, tl_coef, tl_vcv)

make_trans = function(par, dep=TRUE){
  out=array(0, dim=c(18,4,4))
  out[,1,1] = 1
  out[,1,2] <- out[,1,3] <- exp(par[1]+par[2]*c(0:17))
  out[,1,4] = exp(2*par[1]+2*par[2]*c(0:17)+dep*par[3])
  out[,2,2] = 1
  out[,2,4] = exp(par[1]+par[2]*c(0:17))
  out[,3,3] = 1
  out[,3,4] = exp(par[1]+par[2]*c(0:17))
  out[,4,4] = 1
  for(i in 1:18){out[i,,] = sweep(out[i,,], 1, rowSums(out[i,,]), "/")}
  return(out)
}

single_loss = function(par, status="double"){
  if(status=="double"){
    return(apply(make_trans(par), 1, function(x)sum(x[1,c(2,4)]))) 
  } else if(status=="single"){
    return(apply(make_trans(par), 1, function(x)x[2,4]))
  } else stop("unknown status")
}
double_loss = function(par,dep,pdf=FALSE, state=4){
  psi=make_trans(par,dep)
  out=rep(NA,dim(psi)[1])
  x=t(matrix(c(1,0,0,0)))
  out[1]=x[state]
  for(i in 2:(dim(psi)[1])){
    x=x%*%psi[i-1,,]
    out[i]=x[state]
  }
  if(!pdf) {
    return(out)
  } else return(c(0,diff(out)))
}
tag_data = data.frame(Age=c(0:17), single_loss=single_loss(tl_coef), status="double", lower=NA, upper=NA)
tag_data = rbind(tag_data, data.frame(Age=c(0:17), single_loss=single_loss(tl_coef,"single"), status="single", lower=NA, upper=NA))
xxx = apply(tl_rep, 1, single_loss)
tag_data$lower[tag_data$status=="double"] = apply(xxx,1,quantile, prob=0.025)
tag_data$upper[tag_data$status=="double"] = apply(xxx,1,quantile, prob=0.975)
xxx = apply(tl_rep, 1, single_loss, status="single")
tag_data$lower[tag_data$status=="single"] = apply(xxx,1,quantile, prob=0.025)
tag_data$upper[tag_data$status=="single"] = apply(xxx,1,quantile, prob=0.975)

p_tag_loss = ggplot(data=tag_data) + geom_line(aes(x=Age, y=single_loss, color=status)) + geom_ribbon(aes(x=Age, ymin=lower, ymax=upper, fill=status), alpha=0.33) +
  ylab("Probability of tag loss\n") 
ggsave("p_tag_loss.pdf", p_tag_loss, width=5, height=3)

double_loss_data = data.frame(Age=c(0:17))
double_loss_data$dl = double_loss(tl_coef,T)
xxx = apply(tl_rep, 1, double_loss, dep=TRUE)
double_loss_data$lower = apply(xxx,1,quantile, prob=0.025)
double_loss_data$upper = apply(xxx,1,quantile, prob=0.975)
double_loss_data$type = "CDF"

xxx = apply(tl_rep, 1, double_loss, dep=TRUE, pdf=TRUE)

double_loss_data = rbind(double_loss_data, 
                         data.frame(
                           Age=c(0:17),
                           dl = double_loss(tl_coef, T, T),
                           lower = apply(xxx,1, quantile, prob=0.025),
                           upper = apply(xxx,1, quantile, prob=0.975),
                           type = "PDF"
                         ))

p_double_loss = ggplot(data=double_loss_data)+geom_line(aes(x=Age,y=dl))+
  facet_wrap(~type, nrow = 2, scales = "free_y") + 
  geom_ribbon(aes(x=Age, ymin=lower, ymax=upper), alpha=0.33) + ylab("Probability of double loss\n")
ggsave("p_double_loss.pdf",p_double_loss, width=5, height = 5)

# marginal loss
double_loss(tl_coef, T, F, 2) + double_loss(tl_coef, T, F, 4)
xxx = apply(tl_rep, 1, double_loss, dep=TRUE, pdf=FALSE, 2)+ apply(tl_rep, 1, double_loss, dep=TRUE, pdf=FALSE, 4)
t(apply(xxx, 1, quantile, prob=c(0.025, 0.975)))
dl_mode = apply(xxx, 2, function(x)which(x==max(x))-1)

# odds of losing second tag given other was lost
tl_dep = exp(tl_coef[3])
tl_dep_ci = exp(tl_coef[3] +c(-1,1)*1.96*sqrt(tl_vcv[3,3]))


save(list=ls(), file="MVMS_zc_results_summary.RData")