#####  function to build transition matrices

matrix.builder <- function(data){
attach(data)  
#### deterministic matrix
matD = c((1-gH),gH*(1-gTu)*(1-gPu),gH*gTu*(1-gPu),gH*(1-gTu)*gPu,gH*gTu*gPu);
matD = cbind(matD,c((eH),(1-eH)*(1-gTp)*(1-gPt),(1-eH)*gTp*(1-gPt),(1-eH)*(1-gTp)*gPt,(1-eH)*gTp*gPt));
matD = cbind(matD,c((eH),(1-eH)*eTp*(1-gPT),(1-eH)*(1-eTp)*(1-gPT),(1-eH)*eTp*gPT,(1-eH)*(1-eTp)*gPT));
matD = cbind(matD,c((eH),(1-eH)*(1-gTP)*ePt,(1-eH)*gTP*ePt,(1-eH)*(1-gTP)*(1-ePt),(1-eH)*gTP*(1-ePt)));
matD = cbind(matD,c((eH),(1-eH)*eTP*ePT,(1-eH)*(1-eTP)*ePT,(1-eH)*eTP*(1-ePT),(1-eH)*(1-eTP)*(1-ePT)))
#### Stochastic matrices
gHY = c(gH2004,gH2005,gH2006,gH2007,gH2008,gH2009);eHY = c(eH2004,eH2005,eH2006,eH2007,eH2008,eH2009);matR = array(0,c(5,5,6))
for (a in 1:6){gH = gHY[a];eH = eHY[a];
mat = c((1-gH),gH*(1-gTu)*(1-gPu),gH*gTu*(1-gPu),gH*(1-gTu)*gPu,gH*gTu*gPu);
mat = cbind(mat,c((eH),(1-eH)*(1-gTp)*(1-gPt),(1-eH)*gTp*(1-gPt),(1-eH)*(1-gTp)*gPt,(1-eH)*gTp*gPt));
mat = cbind(mat,c((eH),(1-eH)*eTp*(1-gPT),(1-eH)*(1-eTp)*(1-gPT),(1-eH)*eTp*gPT,(1-eH)*(1-eTp)*gPT));
mat = cbind(mat,c((eH),(1-eH)*(1-gTP)*ePt,(1-eH)*gTP*ePt,(1-eH)*(1-gTP)*(1-ePt),(1-eH)*gTP*(1-ePt)));
mat = cbind(mat,c((eH),(1-eH)*eTP*ePT,(1-eH)*(1-eTP)*ePT,(1-eH)*eTP*(1-ePT),(1-eH)*(1-eTP)*(1-ePT)))
               matR[,,a] = mat}
out <- list(matD = matD, matR = matR)
detach(data)
return(out)
}

#####  parameter estimates used in model are found in TransProbs.csv file

TP <- read.csv("TransProbs.csv")  ###   need to put data file in the working directory

#####  Build the transition matrices from underlying parameters

perD <- matrix.builder(TP[1,])$matD
perT <- matrix.builder(TP[1,])$matR
ephD <- matrix.builder(TP[2,])$matD
ephT <- matrix.builder(TP[2,])$matR

#### calculating the stationary state distribution and sensitivities in Figure 1

SSD.ephD <- sens(ephD)$SSD
SSD.perD <- sens(perD)$SSD

dedR.ephD <- sens(ephD)$sensitivities
dedR.perD <- sens(perD)$sensitivities

dedR.ephD.proportional <- markov.sens(ephD)$prop # fig 1
dedR.perD.proportional <- markov.sens(perD)$prop # fig 1


### derivatives of transition rates to each of the sub-parameters

dR.dTheta <- array(0,c(5,5,12))
cb <- c("gH","eH","gTu","gTp","gTP","eTp","eTP","gPu","gPt","gPT","ePt","ePT")
dR.dTheta[1,1,] <- attr(eval(deriv(~(1-gH),cb)),"gradient")
dR.dTheta[2,1,] <- attr(eval(deriv(~gH*(1-gTu)*(1-gPu),cb)),"gradient")
dR.dTheta[3,1,] <- attr(eval(deriv(~gH*gTu*(1-gPu),cb)),"gradient")
dR.dTheta[4,1,] <- attr(eval(deriv(~gH*(1-gTu)*gPu,cb)),"gradient")
dR.dTheta[5,1,] <- attr(eval(deriv(~gH*gTu*gPu,cb)),"gradient")

dR.dTheta[1,2,] <- attr(eval(deriv(~eH,cb)),"gradient")
dR.dTheta[2,2,] <- attr(eval(deriv(~(1-eH)*(1-gTp)*(1-gPt),cb)),"gradient")
dR.dTheta[3,2,] <- attr(eval(deriv(~(1-eH)*gTp*(1-gPt),cb)),"gradient")
dR.dTheta[4,2,] <- attr(eval(deriv(~(1-eH)*(1-gTp)*gPt,cb)),"gradient")
dR.dTheta[5,2,] <- attr(eval(deriv(~(1-eH)*gTp*gPt,cb)),"gradient")

dR.dTheta[1,3,] <- attr(eval(deriv(~eH,cb)),"gradient")
dR.dTheta[2,3,] <- attr(eval(deriv(~(1-eH)*eTp*(1-gPT),cb)),"gradient")
dR.dTheta[3,3,] <- attr(eval(deriv(~(1-eH)*(1-eTp)*(1-gPT),cb)),"gradient")
dR.dTheta[4,3,] <- attr(eval(deriv(~(1-eH)*eTp*gPT,cb)),"gradient")
dR.dTheta[5,3,] <- attr(eval(deriv(~(1-eH)*(1-eTp)*gPT,cb)),"gradient")

dR.dTheta[1,4,] <- attr(eval(deriv(~eH,cb)),"gradient")
dR.dTheta[2,4,] <- attr(eval(deriv(~(1-eH)*(1-gTP)*ePt,cb)),"gradient")
dR.dTheta[3,4,] <- attr(eval(deriv(~(1-eH)*gTP*ePt,cb)),"gradient")
dR.dTheta[4,4,] <- attr(eval(deriv(~(1-eH)*(1-gTP)*(1-ePt),cb)),"gradient")
dR.dTheta[5,4,] <- attr(eval(deriv(~(1-eH)*gTP*(1-ePt),cb)),"gradient")

dR.dTheta[1,5,] <- attr(eval(deriv(~eH,cb)),"gradient")
dR.dTheta[2,5,] <- attr(eval(deriv(~(1-eH)*eTP*ePT,cb)),"gradient")
dR.dTheta[3,5,] <- attr(eval(deriv(~(1-eH)*(1-eTP)*ePT,cb)),"gradient")
dR.dTheta[4,5,] <- attr(eval(deriv(~(1-eH)*eTP*(1-ePT),cb)),"gradient")
dR.dTheta[5,5,] <- attr(eval(deriv(~(1-eH)*(1-eTP)*(1-ePT),cb)),"gradient")

#####  sensitivities stationary state distribution to each of the lower-level parameters

####  Sensitivities for sub-parameters in Table 2

dedCB.ephD <- matrix(0,12,5)
for (a in 1:5){
  for (b in 1:12){
		dedCB.ephD[b,a] <- sum(dR.dTheta[,,b]*dedR.ephD[,,a])
}}
rownames(dedCB.ephD) = cb
colnames(dedCB.ephD) = c("e1","e2","e3","e4","e5")

dedCB.perD <- matrix(0,12,5)
for (a in 1:5){
  for (b in 1:12){
		dedCB.perD[b,a] <- sum(dR.dTheta[,,b]*dedR.perD[,,a])
}}
rownames(dedCB.perD) = cb
colnames(dedCB.perD) = c("e1","e2","e3","e4","e5")

                        
#####  Derived parameters with aggregated lower level parameters
                        
dedCB2.ephD <- rbind(dedCB.ephD[1:2,],colSums(dedCB.ephD[3:5,]),colSums(dedCB.ephD[6:7,]),colSums(dedCB.ephD[8:10,]),colSums(dedCB.ephD[11:12,]))
rownames(dedCB2.ephD) = c("gH","eH","gT","eT","gP","eP")      

dedCB2.perD <- rbind(dedCB.perD[1:2,],colSums(dedCB.perD[3:5,]),colSums(dedCB.perD[6:7,]),colSums(dedCB.perD[8:10,]),colSums(dedCB.perD[11:12,]))
rownames(dedCB2.perD) = c("gH","eH","gT","eT","gP","eP")      
 
dsoloToadsdCB2.perD <- dedCB2.perD[,3]                       
dsoloToadsdCB2.ephD <- dedCB2.ephD[,3]                    
                                                                                            
dcooccurdCB2.perD <- dedCB2.perD[,3]*-SSD.perD[5]/((sum(SSD.perD[c(3,5)]))^2) + dedCB2.perD[,5]*SSD.perD[3]/((sum(SSD.perD[c(3,5)]))^2)                       
dcooccurdCB2.ephD <- dedCB2.ephD[,3]*-SSD.ephD[5]/((sum(SSD.ephD[c(3,5)]))^2) + dedCB2.ephD[,5]*SSD.ephD[3]/((sum(SSD.ephD[c(3,5)]))^2)                       
                                                                                                                                                                                        
######  Temporal sensitivities for derived parameters to lower-level parameters

#####  generate new set of matrices with small perturbation to parameters                        
tp <- TP[1,]   ###  set as row 1 for perennial                       
q <- 0.0001
mat.gH <- matrix.builder(tp + c(rep(q,6),rep(0,18)))$matR   ####  This builds a new site of transition matrices with q added to gH
mat.eH <- matrix.builder(tp + c(rep(0,6),rep(q,6),rep(0,12)))$matR
mat.gT <- matrix.builder(tp + c(rep(0,14),rep(q,3),rep(0,7)))$matR   
mat.eT <- matrix.builder(tp + c(rep(0,17),rep(q,2),rep(0,5)))$matR
mat.gP <- matrix.builder(tp + c(rep(0,19),rep(q,3),rep(0,3)))$matR                        
mat.eP <- matrix.builder(tp + c(rep(0,22),rep(q,2)))$matR        

#####  numerical aproximation of the sensitivities of the stationary state distribution to a lower-level parameter                         
#####  the function lower.Tsens works by inputting the original and perturbed matrices (R)
                         
deT.dgH <- lower.Tsens(perT,mat.gH,q)
deT.deH <- lower.Tsens(perT,mat.eH,q)
deT.dgT <- lower.Tsens(perT,mat.gT,q)
deT.deT <- lower.Tsens(perT,mat.eT,q)
deT.dgP <- lower.Tsens(perT,mat.gP,q)
deT.deP <- lower.Tsens(perT,mat.eP,q)

per.deT.dlower <- data.frame(deT.dgH$sens,deT.deH$sens,deT.dgT$sens,deT.deT$sens,deT.dgP$sens,deT.deP$sens)

####  in the case of a derived parameters we can extract the sequence of state distributions at each time step and calculate the derived parameter for each step                        

dcooccurT.dgH <- (deT.dgH$p.T[5]/(deT.dgH$p.T[3]+deT.dgH$p.T[5])-deT.dgH$e.T[5]/(deT.dgH$e.T[3]+deT.dgH$e.T[5]))/q
dcooccurT.deH <- (deT.deH$p.T[5]/(deT.deH$p.T[3]+deT.deH$p.T[5])-deT.deH$e.T[5]/(deT.deH$e.T[3]+deT.deH$e.T[5]))/q
dcooccurT.dgT <- (deT.dgT$p.T[5]/(deT.dgT$p.T[3]+deT.dgT$p.T[5])-deT.dgT$e.T[5]/(deT.dgT$e.T[3]+deT.dgT$e.T[5]))/q
dcooccurT.deT <- (deT.deT$p.T[5]/(deT.deT$p.T[3]+deT.deT$p.T[5])-deT.deT$e.T[5]/(deT.deT$e.T[3]+deT.deT$e.T[5]))/q
dcooccurT.dgP <- (deT.dgP$p.T[5]/(deT.dgP$p.T[3]+deT.dgP$p.T[5])-deT.dgP$e.T[5]/(deT.dgP$e.T[3]+deT.dgP$e.T[5]))/q
dcooccurT.deP <- (deT.deP$p.T[5]/(deT.deP$p.T[3]+deT.deP$p.T[5])-deT.deP$e.T[5]/(deT.deP$e.T[3]+deT.deP$e.T[5]))/q
per.dcooccurT.dlower <- data.frame (dcooccurT.dgH,dcooccurT.deH,dcooccurT.dgT,dcooccurT.deT,dcooccurT.dgP,dcooccurT.deP)  
###  Because the derived parameter total toads is simple addition it can be calculated by
###  adding the sensitivities for 3rd and 5th states rather than going back to original simulated states
per.dsoloToadsT.dlower <- per.deT.dlower[3,]

####  SAME FOR EPHEMERAL                      
tp <- TP[2,]   ###  set as row 2 for ephemeral                       
q <- q.eph/10000
mat.gH <- matrix.builder(tp + c(rep(q,6),rep(0,18)))$matR
mat.eH <- matrix.builder(tp + c(rep(0,6),rep(q,6),rep(0,12)))$matR
mat.gT <- matrix.builder(tp + c(rep(0,14),rep(q,3),rep(0,7)))$matR   
mat.eT <- matrix.builder(tp + c(rep(0,17),rep(q,2),rep(0,5)))$matR
mat.gP <- matrix.builder(tp + c(rep(0,19),rep(q,3),rep(0,3)))$matR                        
mat.eP <- matrix.builder(tp + c(rep(0,22),rep(q,2)))$matR        

deT.dgH <- lower.Tsens(ephT,mat.gH,q)
deT.deH <- lower.Tsens(ephT,mat.eH,q)
deT.dgT <- lower.Tsens(ephT,mat.gT,q)
deT.deT <- lower.Tsens(ephT,mat.eT,q)
deT.dgP <- lower.Tsens(ephT,mat.gP,q)
deT.deP <- lower.Tsens(ephT,mat.eP,q)
eph.deT.dlower <- data.frame(deT.dgH$sens,deT.deH$sens,deT.dgT$sens,deT.deT$sens,deT.dgP$sens,deT.deP$sens)

dcooccurT.dgH <- (deT.dgH$p.T[5]/(deT.dgH$p.T[3]+deT.dgH$p.T[5])-deT.dgH$e.T[5]/(deT.dgH$e.T[3]+deT.dgH$e.T[5]))/q
dcooccurT.deH <- (deT.deH$p.T[5]/(deT.deH$p.T[3]+deT.deH$p.T[5])-deT.deH$e.T[5]/(deT.deH$e.T[3]+deT.deH$e.T[5]))/q
dcooccurT.dgT <- (deT.dgT$p.T[5]/(deT.dgT$p.T[3]+deT.dgT$p.T[5])-deT.dgT$e.T[5]/(deT.dgT$e.T[3]+deT.dgT$e.T[5]))/q
dcooccurT.deT <- (deT.deT$p.T[5]/(deT.deT$p.T[3]+deT.deT$p.T[5])-deT.deT$e.T[5]/(deT.deT$e.T[3]+deT.deT$e.T[5]))/q
dcooccurT.dgP <- (deT.dgP$p.T[5]/(deT.dgP$p.T[3]+deT.dgP$p.T[5])-deT.dgP$e.T[5]/(deT.dgP$e.T[3]+deT.dgP$e.T[5]))/q
dcooccurT.deP <- (deT.deP$p.T[5]/(deT.deP$p.T[3]+deT.deP$p.T[5])-deT.deP$e.T[5]/(deT.deP$e.T[3]+deT.deP$e.T[5]))/q
eph.dcooccurT.dlower <- data.frame (dcooccurT.dgH,dcooccurT.deH,dcooccurT.dgT,dcooccurT.deT,dcooccurT.dgP,dcooccurT.deP)  
eph.dsoloToadsT.dlower <- eph.deT.dlower[3,]

####  data used in figure 2

fig2.data <- matrix(c(dcooccurdCB2.perD,per.dcooccurT.dlower,dsoloToadsdCB2.perD,per.dsoloToadsT.dlower,dcooccurdCB2.ephD,eph.dcooccurT.dlower,dsoloToadsdCB2.ephD,eph.dsoloToadsT.dlower),6,8)