#### functions in Appendix A required to run the analysis #### Golden Eagle Data - Biological Conservation 2009 expit <- function(x) (exp(x)/(1+exp(x))) h.series = c(0.0703,0.0453,0.0251,0.0143,0.0085,0.0070,0.0083,0.0076,0.0354,0.0445,0.0812,0.0775,0.0540,0.0052,0.0087,0.0120,0.0625,0.2520,0.2620) #p0i = -0.72631 ### intercept logit function for psi-0 #p1i = 2.559808 ### intercept logit function for psi-1 #p2i = 3.090135 ### intercept logit function for psi-2 #p0d = -1.41284 ### slope for disturbance effect logit function for psi-0 #p1d = 0.556733 ### slope for disturbance effect logit function for psi-1 #p2d = 0.100006 ### slope for disturbance effect logit function for psi-2 #Ri = -0.775 ### intercept logit function for R-0 #Rh = 5.948795 ### slope for hare effect logit function for R-0 #dist = 25/93 ### Proportion of disturbed sites #hare = 0.057 ### Mean annual hare index eagle <- data.frame(p0i = -0.72631, p1i = 2.559808, p2i = 3.090135, p0d = -1.41284, p1d = 0.556733, p2d = 0.100006, Ri = -0.775, Rh = 5.94879, dist = 25/93, hare = 0.057) matrix.builder <- function(data,h.series){ attach(data) PSI0 = expit(p0i + p0d*dist) PSI1 = expit(p1i + p1d*dist) PSI2 = expit(p2i + p2d*dist) R0 = expit(Ri + Rh*hare) R1 = expit(Ri + Rh*hare) R2 = expit(Ri + Rh*hare) #### deterministic matrix matD = c((1-PSI0),(1-PSI1),(1-PSI2)) matD = rbind(matD,c(PSI0*(1-R0),PSI1*(1-R1),PSI2*(1-R2))) matD = rbind(matD,c(PSI0*R0,PSI1*R1,PSI2*R2)) #### Stochastic matrices matR = array(0,c(3,3,19)) for (a in 1:19){ R0 = expit(Ri + Rh*h.series[a]) R1 = expit(Ri + Rh*h.series[a]) R2 = expit(Ri + Rh*h.series[a]) mat = c((1-PSI0),(1-PSI1),(1-PSI2)) mat = rbind(mat,c(PSI0*(1-R0),PSI1*(1-R1),PSI2*(1-R2))) mat = rbind(mat,c(PSI0*R0,PSI1*R1,PSI2*R2)) matR[,,a] = mat} out <- list(matD = matD, matR = matR) detach(data) return(out) } t.mat <- matrix.builder(eagle,h.series)\$matD #### Sensitvities of Stable State Distribution to changes in Transition Rates and Stable State Distribution t.sens = sens(t.mat)\$sensitivities ### Equation 6 W = sens(t.mat)\$SSD ### Equation 4 #### Sensitivities of transition rates using Uniform compensation dW.dTU <- markov.sens(t.mat)\$unif ### Equation 7 #### Sensitivities of transition rates using proportional compensation dW.dTP <- markov.sens(t.mat)\$prop ### Equation 8 #### Derivative matrices for each of the parameters used in the conditional binomials - like Equation 12 dR.dTheta <- array(0,c(3,3,6)) CB <- c("R0","R1","R2","PSI0","PSI1","PSI2") dR.dTheta[1,1,] <- attr(eval(deriv(~1-PSI0,CB)),"gradient") dR.dTheta[1,2,] <- attr(eval(deriv(~1-PSI1,CB)),"gradient") dR.dTheta[1,3,] <- attr(eval(deriv(~1-PSI2,CB)),"gradient") dR.dTheta[2,1,] <- attr(eval(deriv(~PSI0*(1-R0),CB)),"gradient") dR.dTheta[2,2,] <- attr(eval(deriv(~PSI1*(1-R1),CB)),"gradient") dR.dTheta[2,3,] <- attr(eval(deriv(~PSI2*(1-R2),CB)),"gradient") dR.dTheta[3,1,] <- attr(eval(deriv(~PSI0*R0,CB)),"gradient") dR.dTheta[3,2,] <- attr(eval(deriv(~PSI1*R1,CB)),"gradient") dR.dTheta[3,3,] <- attr(eval(deriv(~PSI2*R2,CB)),"gradient") #### Sensitivity of Stable Stage Distribution to changes in conditional binomial parameters - Equation 9 dW.dR0 <- c(sum(dR.dTheta[,,1]*t.sens[,,1]),sum(dR.dTheta[,,1]*t.sens[,,2]),sum(dR.dTheta[,,1]*t.sens[,,3])) dW.dR1 <- c(sum(dR.dTheta[,,2]*t.sens[,,1]),sum(dR.dTheta[,,2]*t.sens[,,2]),sum(dR.dTheta[,,2]*t.sens[,,3])) dW.dR2 <- c(sum(dR.dTheta[,,3]*t.sens[,,1]),sum(dR.dTheta[,,3]*t.sens[,,2]),sum(dR.dTheta[,,3]*t.sens[,,3])) dW.dPSI0 <- c(sum(dR.dTheta[,,4]*t.sens[,,1]),sum(dR.dTheta[,,4]*t.sens[,,2]),sum(dR.dTheta[,,4]*t.sens[,,3])) dW.dPSI1 <- c(sum(dR.dTheta[,,5]*t.sens[,,1]),sum(dR.dTheta[,,5]*t.sens[,,2]),sum(dR.dTheta[,,5]*t.sens[,,3])) dW.dPSI2 <- c(sum(dR.dTheta[,,6]*t.sens[,,1]),sum(dR.dTheta[,,6]*t.sens[,,2]),sum(dR.dTheta[,,6]*t.sens[,,3])) #### Derivatives of each of the conditional binomials to changes in covariate #### uses the derivative equation for the logit link - Equation 14 dPSI0.dDIST <- (p0d*exp(p0d*dist+p0i))/((exp(2*p0d*dist+2*p0i))+(2*exp(p0d*dist+p0i))+1) dPSI1.dDIST <- (p1d*exp(p1d*dist+p1i))/((exp(2*p1d*dist+2*p1i))+(2*exp(p1d*dist+p1i))+1) dPSI2.dDIST <- (p2d*exp(p2d*dist+p2i))/((exp(2*p2d*dist+2*p2i))+(2*exp(p2d*dist+p2i))+1) dR.dHARE <- (Rh*exp(Rh*hare+Ri))/((exp(2*Rh*hare+2*Ri))+(2*exp(Rh*hare+Ri))+1) #### Sensitivity of Stable State Distribution to changes in covariates - Equation 15 dW1.dHARE <- sum(dR.dHARE*dW.dR0[1])+ sum(dR.dHARE*dW.dR1[1])+ sum(dR.dHARE*dW.dR2[1]) dW2.dHARE <- sum(dR.dHARE*dW.dR0[2])+ sum(dR.dHARE*dW.dR1[2])+ sum(dR.dHARE*dW.dR2[2]) dW3.dHARE <- sum(dR.dHARE*dW.dR0[3])+ sum(dR.dHARE*dW.dR1[3])+ sum(dR.dHARE*dW.dR2[3]) dW1.dDIST <- sum(dPSI0.dDIST*dW.dPSI0[1])+ sum(dPSI1.dDIST*dW.dPSI1[1])+ sum(dPSI2.dDIST*dW.dPSI2[1]) dW2.dDIST <- sum(dPSI0.dDIST*dW.dPSI0[2])+ sum(dPSI1.dDIST*dW.dPSI1[2])+ sum(dPSI2.dDIST*dW.dPSI2[2]) dW3.dDIST <- sum(dPSI0.dDIST*dW.dPSI0[3])+ sum(dPSI1.dDIST*dW.dPSI1[3])+ sum(dPSI2.dDIST*dW.dPSI2[3]) #### Sensitivity of super-parameters (proportion of occupied sites where reproduction was successful #### and proportion of sites occupied) to changes in covariates dRSUCC.dHARE = dW2.dHARE*-W[3]/((W[2]+W[3])^2) + dW3.dHARE*W[2]/((W[2]+W[3])^2) dOCC.dHARE = dW2.dHARE + dW3.dHARE dRSUCC.dDIST = dW2.dDIST*-W[3]/((W[2]+W[3])^2) + dW3.dDIST*W[2]/((W[2]+W[3])^2) dOCC.dDIST = dW2.dDIST + dW3.dDIST ##### Recreates values in Table B1 in the appendix Table1 <- cbind(as.vector(t.mat),as.vector(dW.dTU[,,1]),as.vector(dW.dTU[,,2]),as.vector(dW.dTU[,,3])) Table1 <- cbind(Table1,as.vector(dW.dTP[,,1]),as.vector(dW.dTP[,,2]),as.vector(dW.dTP[,,3])) ##### TEMPORALLY AND SPATIALLY VARIABLE SENSITIVITIES ##### Sensitivities are calculated for the stationary state distribution where R varies among years according to hare abundance ##### and among site by whether the site is disturbance prone or not ##### Examples give sensitivity of the stationary state distribution to changes in a transition probability (R^[i,j]), ##### a conditional binomial parameter (e.g., psi or B), and to the environmental covariates dist (disturbance) and hare (hare index) matR <- matrix.builder(eagle,h.series)\$matR ## matrices including temporal variation in hare abundance #### Perturbation to the transition probabilities assuming uniform and proportional compensation #### Compare these values to those in Table B1 to see how accounting for temporal variaion changes sensitivities deTdR.unif <- unif.Tsens(matR) deTdR.prop <- prop.Tsens(matR) #### Perturbation to the conditional binomial parameter PSI1 q = 0.0001 eagle2 <- eagle + c(0,q,rep(0,8)) matR2 <- matrix.builder(eagle2,h.series)\$matR deTdPSI1 <- lower.Tsens(matR,matR2,q)\$sens #### calculations for Table B2 with environmental variability eagle.d0 <- eagle; eagle.d0["dist"]= 0 ### data for sites with no disturbance eagle.d1 <- eagle; eagle.d1["dist"]= 1 ### data for sites with disturbance h.series2 <- h.series + q ### perturbed h.series mat.dist0 <- matrix.builder(eagle.d0,h.series)\$matR mat.dist1 <- matrix.builder(eagle.d1,h.series)\$matR mat.dist0p <- matrix.builder(eagle.d0,h.series2)\$matR mat.dist1p <- matrix.builder(eagle.d1,h.series2)\$matR ssd.d0 <- ssd.T(mat.dist0) ssd.d1 <- ssd.T(mat.dist1) deTSdpropdist <- ssd.d1-ssd.d0 ### sensitivity to changes in the proportion disturbed deTdhare.undisturbed <- lower.Tsens(mat.dist0,mat.dist0p,q)\$sens deTdhare.disturbed <- lower.Tsens(mat.dist1,mat.dist1p,q)\$sens deTSdhare <- 68/93*deTdhare.undisturbed+25/93*deTdhare.disturbed ### sensitivity to changes in hare index