#####  Sensitivity analysis of probability of a patch being occupied to changes in the patch area
#####  Example comes from Hanski 1994 incidence function model where rescue effects are not included
#####  Sensitivity is for a single patch of mean size

rel.A = 1  	###  This parameter can be used to examine the sensitivity of patch occupancy 
           		###  to increases in patch size for patches of different relative 
           		###  size compared to the mean patch size.

####  Species 1 – M. cinxia

e = 0.044		### estimated parameter from Hanski 1994
x = 0.501		### estimated parameter	
A = 0.41*rel.A		### mean patch size
E = e/(A^x)		### predicted extinction rate for mean size patch 

D1 = c()		### empty vector to hold output
for (a in 1:100){
C = a*0.01				### colonization rate varied from 0.01 to 1
R <- matrix(c((1-C),C,E,(1-E)),2,2)	### transition matrix (eq. 16)
R.sens <- sens(R)$sensitivities		### sensitivities of matrix elements
dT.dE <- matrix(c(0,0,1,-1),2,2)		### derivative of transition rate elements to extinction
dT.dC <- matrix(c(-1,1,0,0),2,2)		### derivative of transition rate elements to colonization
dE.dA <- -e*x*(A^(-x-1))			### derivative of extinction to patch size (eq. 17)
dOCC.dA <- sum(dE.dA*dT.dE*R.sens[,,2])	### full derivative of occupancy to area (eq. 18)
D1[a] = dOCC.dA 			
}

####  Species 2 – H. comma

e = 0.010
x = 1.009
A = 0.94*rel.A
E = e/(A^x)

D2 = c()
for (a in 1:100){
C = a*0.01
R <- matrix(c((1-C),C,E,(1-E)),2,2)
R.sens <- sens(R)$sensitivities
dT.dE <- matrix(c(0,0,1,-1),2,2)
dT.dC <- matrix(c(-1,1,0,0),2,2)
dE.dA <- -e*x*(A^(-x-1))
dOCC.dA <- sum(dE.dA*dT.dE*R.sens[,,2])
D2[a] = dOCC.dA 
}

####  Species 3 – S. orion

e = 0.036
x = 0.790
A = 0.10*rel.A


E = e/(A^x)

D3 = c()
for (a in 1:100){
C = a*0.01
R <- matrix(c((1-C),C,E,(1-E)),2,2)
R.sens <- sens(R)$sensitivities
dT.dE <- matrix(c(0,0,1,-1),2,2)
dT.dC <- matrix(c(-1,1,0,0),2,2)
dE.dA <- -e*x*(A^(-x-1))
dOCC.dA <- sum(dE.dA*dT.dE*R.sens[,,2])
D3[a] = dOCC.dA 
}

#####   Recreates plots from figure 1 in the text

par(mfrow = c(3,1))
plot(c(1:100)/100,D1, type = "l",xlab = "M. cinxia")
plot(c(1:100)/100,D2, type = "l",xlab = "H. comma")
plot(c(1:100)/100,D3, type = "l",xlab = "S. orion")