#################################################################### # # ArtemisiaRPartition: partitions the variance of lifetime # reproductive output, lifespan, and size at death into a contribution # from trait variation and a contribution from luck. # # X = trait, Z = state in matrix model, age x repr. state classification # The traits here are logit(survival) for the 4 adult breeding states # Juvenile survival is 1 (by definition; onluy survivors to adulthood are counted) ####################################################################### rm(list=ls(all=TRUE)); graphics.off(); require(Matrix); require(statmod); path ="/home/robin/ipm/drivers/Rpartition/empirical" out=try(setwd(path),silent=TRUE); if(class(out)=="try-error") { path=ifelse(.Platform$OS.type=="windows","c:/repos/drivers/Rpartition","~/repos/drivers/Rpartition"); setwd(path); } source ("KittiwakeMatrices.R") source("Standard Graphical Pars.R") source("Utilities.R") source("MegamatrixFunctions.R"); numTraitVals = 61 mz = 25; # size of the age/stage transition matrix ## Prob. of breeding pb = c(0, 0, 0, 1, 1) # prob. of breeding ## Variance in per capita # of recruits. ## Assumes those with 2 or 3 recruits have equal prob of 2 or 3. sigbsq = c(0, 0, 0, 0, 0.25) # rough estimates for kittiwake survival based on Cam et al. (2002) results q75 = log(0.9/0.1); #75th percentile on inverse-logit scale ("quality") q25 = log(0.7/0.3) #25th percentile on inverse-logit scale ("quality") IQR=q75-q25; sigma=IQR/1.35; # mu = (q75+q25)/2; # estimated mean, sd, and CV of quality CVest = sigma/mu; # ## what is the corresponding CV in survival (instead of "quality")? survs = rnorm(100000,mean=mu,sd=sigma); survs = exp(survs)/(1+exp(survs)); CVest.surv=sd(survs)/mean(survs) sdSurvQuality = CVest * pars minSurvQuality = pars - 2.5*sdSurvQuality maxSurvQuality = pars + 2.5*sdSurvQuality ## matrix of survival qualities survQualityDist = matrix (0, numTraitVals, 4) for (i in 1:4) { survQualityDist[,i] = seq(minSurvQuality[i], maxSurvQuality[i], length.out=numTraitVals) } ## vector of associated probabilities survQualityProb = dnorm (survQualityDist[,1], mean=pars[1],sd=sdSurvQuality[1]) survQualityProb = survQualityProb/sum(survQualityProb) ################################################## # Variation in breeding probabilities ################################################## ## estimate CV of breeding prob. from Cam et al. 2002 q25 = 0.94 q75 = 0.98 IQR=q75-q25; sigma=IQR/1.35; mu = (q75+q25)/2; CVBreedingProb = sigma/mu; ## This is the range of breedingProbFactor that we need to get a +/- ## 1.5 std dev. variation in breeding probability. minBreedingProbFactor = 1 - 25.92*CVBreedingProb maxBreedingProbFactor = 1 + 25.92*CVBreedingProb ## Breeding prob. adjustment factors breedingProbFactorDist = seq(minBreedingProbFactor, maxBreedingProbFactor, length.out=numTraitVals) expRCondXZ = varRCondXZ = expLCondXZ = varLCondXZ = matrix (0, mz, numTraitVals) ## Loop over different trait values; the trait is survival and breeding ## probability "quality" for (ii in 1:numTraitVals) { ## Set survival qualities, breeding probs. newpars = survQualityDist[ii,] breedingProbFactor = breedingProbFactorDist[ii] ## Modify reproductive stage transition matrices. (Only M5 will ## actually change, the way we have things now.) newM3 = M3; newM4 = M4; newM5 = M5; newM5[4:5,] = newM5[4:5,] * breedingProbFactor for (j in 1:5) { B = sum(newM5[4:5,j]) ## flux into breeding stages A = sum(newM5[1:3,j]) ## flux into non-breeding stages bfac = (1-B)/A; ## alter flux into non-breeding stages so that column sums of M5 ## remain equal to 1. newM5[1:3,j]=newM5[1:3,j]*bfac; } ## Make the growth/survival matrix P s = c(1, exp(newpars)/(1 + exp(newpars))) P = matrix (0, 25, 25) P[5*1 + 1:5, 5*0 + 1:5] = M1 * matrix (rep(s1, 5), 5, 5, byrow=TRUE) P[5*2 + 1:5, 5*1 + 1:5] = M2 * matrix (rep(s2, 5), 5, 5, byrow=TRUE) P[5*3 + 1:5, 5*2 + 1:5] = newM3 * matrix (rep(s3, 5), 5, 5, byrow=TRUE) P[5*4 + 1:5, 5*3 + 1:5] = newM4 * matrix (rep(s, 5), 5, 5, byrow=TRUE) P[5*4 + 1:5, 5*4 + 1:5] = newM5 * matrix (rep(s, 5), 5, 5, byrow=TRUE) N = solve (diag(mz) - P) ## expected lifetime reproduction conditional on state expRCondXZ[,ii] = apply (F %*% N, 2, sum) ## Expected lifetime conditional on state (for the current trait value) expLCondXZ[,ii] = apply (N, 2, sum) ## variance of lifetime reproduction conditional on state (for the current trait value) rbarPib = expRCondXZ[,ii] %*% P # \bar{r} \pi_b r2 = (sigbsq + (pb*b)^2 + 2*pb*b*rbarPib) %*% N varRCondXZ[,ii] = r2 - expRCondXZ[,ii]^2 ## Variance in lifespan conditional on trait value varLCondXZ[,ii] = apply (2*N %*% N - N, 2, sum) - expLCondXZ[,ii]^2 } ## Set up birth size distribution - all are immature, age 1. initSizeDist = c(1, rep(0, mz-1)); ## eq. 2 in Kittiwake-CalculatingVarR.tex expZVarRCondXZ = initSizeDist %*% varRCondXZ varZExpRCondXZ = initSizeDist %*%(expRCondXZ^2) - (initSizeDist %*% expRCondXZ)^2 varRCondX = expZVarRCondXZ + varZExpRCondXZ expZVarLCondXZ = initSizeDist %*% varLCondXZ varZExpLCondXZ = initSizeDist %*% expLCondXZ^2 - (initSizeDist %*% expLCondXZ)^2 varLCondX = expZVarLCondXZ + varZExpLCondXZ ## eq. 3 in Kittiwake-CalculatingVarR.tex expRCondX = initSizeDist %*% expRCondXZ expLCondX = initSizeDist %*% expLCondXZ ## eq. 1 in Kittiwake-CalculatingVarR.tex expXVarRCondX = sum (survQualityProb * varRCondX) varXExpRCondX = sum (survQualityProb * expRCondX^2) - sum(survQualityProb * expRCondX)^2 expXVarLCondX = sum (survQualityProb * varLCondX) varXExpLCondX = sum (survQualityProb * expLCondX^2) - sum (survQualityProb * expLCondX)^2 luckContrib = expXVarRCondX traitContrib = varXExpRCondX luckContribL = expXVarLCondX traitContribL = varXExpLCondX traitPercentageR = traitContrib / (traitContrib + luckContrib) traitPercentageL = traitContribL / (traitContribL + luckContribL) cat ("Proportion of V(R) from trait var. = ", traitPercentageR, "\n") cat ("Proportion of V(L) from trait var. = ", traitPercentageL, "\n")