#################################################################### # # KittiwakeWinnerTraits: Calculates P(survival "quality" | R) for the # Kittiwake model. Note that survival is varying simultaneously for # stages 2--5. To do: can we also vary the probability of breeding in # a positively covarying way, per Cam et al. 2002? The breeding # states are 4, 5, 9, 10, 14, 15, 19, 20, 24, 25. # ####################################################################### rm(list=ls(all=TRUE)); graphics.off(); require(Matrix); require(statmod); source ("KittiwakeMatrices.R") source("Standard Graphical Pars.R") source("Utilities.R") source("MegamatrixFunctions.R"); #source("MatrixImage.R") maxKids = 80 mT = maxKids + 1 numTraitVals = 61 debugging = FALSE getProbRCondTrait = function (P) { B = mk_B() ## Construct A, the transition matrix for a (number of kids) x stage x ## age model. out = make_AxT (B, P, mT) A = out$A; ## 2D iteration matrix ## Create a survival probability vector surv = apply (A, 2, sum); die = 1-surv; ## And the fundamental matrix N <- solve(diag(ncol(A))-A) ## Make omega omega = matrix(0,nrow(N),ncol(N)); for(j in 1:ncol(N)) {omega[,j] = die*N[,j]} distAtBirth = matrix(0,25,mT) distAtBirth[1,1] = 1 # All chicks start off at age 1, stage 1, no kids distAtDeath = omega%*%vec(distAtBirth) # state vector distAtDeath = unvec(distAtDeath,nrow=25) # state matrix distKidsAtDeath = apply(distAtDeath,2,sum); return (distKidsAtDeath) } ################################################## # Variation in survival probabilities ################################################## ## rough estimates for kittiwake survival based on Cam et al. (2002) results q75 = 0.9 #75th percentile q25 = 0.7 #25th percentile IQR=q75-q25; sigma=IQR/1.35; mu = (q75+q25)/2; CVest = sigma/mu; # CV of actual survival 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; CVest = sigma/mu; 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) ## why 18*sigma? We want breeding probability to vary over an ## interval of 3*sigma=0.0888, which corresponds to varying the ## breeding probability factor over an interval of 1.6. If we want ## the breeding probability factor interval to correspond to 3 times ## the standard deviation of the breeding probability factor ## distribution, then it seems like the sd of the breeding probability ## factor should be 1.6/0.0888 * sigma. breedingProbFactorProb = dnorm (breedingProbFactorDist, mean=1, sd=18*sigma) breedingProbFactorProb = breedingProbFactorProb/sum(breedingProbFactorProb) if (debugging) { ## Verify that the breeding probability varies between 0.6455 - ## 1.5*sigma = 0.601 and 0.6455 + 1.5*sigma = 0.6899. Recall that ## Cam et al. consider failed breeding attempts as "breeding" and ## that breeding probability is conditional upon survival, so we ## should use the M matrices, not P, to determine flux into breeding ## states. (We should, however, use P to determine the stable age x ## stage structure.) foo = c(minBreedingProbFactor, maxBreedingProbFactor) bpBounds = rep(0,2) for (jj in 1:2) { breedingProbFactor = foo[jj] ## Modify reproductive stage transition matrices newM3 = M3; newM4 = M4; newM5 = M5; newM5[4:5,] = newM5[4:5,] * breedingProbFactor for (j in 1:5) { B = sum(newM5[4:5,j]); A = sum(newM5[1:3,j]); bfac = (1-B)/A; newM5[1:3,j]=newM5[1:3,j]*bfac; } P = makeP(pars) A = P + F w = domEig(A, verbose=FALSE)$w ## Flux into age 4 fluxes3 = newM3 %*% w[11:15] ## Flux into age 5 fluxes4 = newM4 %*% w[16:20] ## Flux into age 6 and older fluxes5 = newM5 %*% w[21:25] #bpBounds[jj] = (sum(fluxes3[3:5]) + sum(fluxes4[3:5]) + sum(fluxes5[3:5])) / # (sum(fluxes3) + sum(fluxes4) + sum(fluxes5)) bpBounds[jj] = (sum(fluxes4[3:5]) + sum(fluxes5[3:5]))/(sum(fluxes4) + sum(fluxes5)) } cat (bpBounds, "\n") } ## dist. of Prob(R | surv) for each value of surv probRCondSurv = matrix (0, mT, numTraitVals) ## dist. of Prob(surv | R) for each value of surv probSurvCondR = matrix (0, mT, numTraitVals) for (i in 1:numTraitVals) { ## Set survival qualities, breeding probs. newpars = survQualityDist[i,] breedingProbFactor = breedingProbFactorDist[i] ## Modify reproductive stage transition matrices newM3 = M3; newM4 = M4; newM5 = M5; newM5[4:5,] = newM5[4:5,] * breedingProbFactor for (j in 1:5) { B = sum(newM5[4:5,j]); A = sum(newM5[1:3,j]); bfac = (1-B)/A; newM5[1:3,j]=newM5[1:3,j]*bfac; } 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) ## Fill in probRCondSurv probRCondSurv[,i] = getProbRCondTrait(P) cat (i, "\n") } numerator = matrix (0, mT, numTraitVals) denominator = numeric (numTraitVals) ## numerator = P(R | surv, breedingProb) P(surv, breedingProb). for (i in 1:numTraitVals) ## if variation in survival and breeding probability is perfectly ## correlated, then P(surv, breedingProb) = P(surv). numerator[,i] = probRCondSurv[,i] * survQualityProb[i] ## denominator = P(R) probR = apply (numerator, 1, sum) ## What are the 25th and 75th percentiles for P(R)? R25 = min(which(cumsum(probR) >= 0.25)) R50 = min(which(cumsum(probR) >= 0.5)) R75 = min(which(cumsum(probR) >= 0.75)) R90 = min(which(cumsum(probR) >= 0.9)) R95 = min(which(cumsum(probR) >= 0.95)) R99 = min(which(cumsum(probR) >= 0.99)) plotR = seq(1, R95 + 1, by=1) ## What are the quantiles for survival quality? Q01 = min(which(cumsum(survQualityProb) >= 0.01)) Q05 = min(which(cumsum(survQualityProb) >= 0.05)) Q95 = min(which(cumsum(survQualityProb) >= 0.95)) Q99 = min(which(cumsum(survQualityProb) >= 0.99)) dp = abs(cumsum(survQualityProb)-0.5) Q50 = which(dp==min(dp)); plotQ=seq(Q01,Q99,by=1); probSurvCondR = matrix (0, mT, numTraitVals) for (i in 1:numTraitVals) probSurvCondR[,i] = numerator[,i]/probR layout (mat=rbind(matrix(1, 7, 6), rep(2, 6))) par(mgp=c(2,1,0), cex.axis=1.3, cex.lab=1.4); par (mar=c(2,6,3,1)) image (plotR, plotQ, probSurvCondR[plotR,plotQ], axes=FALSE, xlab="", ylab="") mtext ("x=Inverse logit survival", side=2, line=3.5, cex=1.3) axis (side=2, at=c(Q05,Q50,Q95), label=c(" 5th %tile", "Mean","95th %tile ")) axis (side=1, seq(1, 30, by=5), at=seq(1, 30, by=5)) #contour (plotR, 1:numTraitVals, probSurvCondR[plotR,], add=T) contour (plotR, plotQ, probSurvCondR[plotR,plotQ], add=T) abline (h=Q50, lty=2) ## add a skinny probability density for P(R) below the plot par (mar=c(2,5,0,1)) plot (plotR, probR[1:length(plotR)], axes=FALSE, type="l", lwd=2, xlab="", ylab="") axis (side=1, at=c(R25, R75, R90, R95), label=c("25th %tile", "75th %tile", "90th %tile", "95th %tile")) mtext ("Lifetime reproductive output", side=3, line=-2, cex=1.3) #dev.copy2eps(file="KittwakeWinnerTraitsFecSurv.eps"); #dev.copy2eps(file="../../Manuscripts/Rpartition/KittiwakeWinnerTraitsFecSurv.eps");