A bunch of functions with Lotka-ish demography formulas

LotkaFunctions.R

# Author: Tim Riffe
###############################################################################


LotkaRanalytic <-
		function(fx,Lx,x){
	R0 <- Rmomentn(fx,Lx,x,0)
	R1 <- Rmomentn(fx,Lx,x,1)
	R2 <- Rmomentn(fx,Lx,x,2)
	# alpha would have been:	 R1/R0	
	# Beta would have been:		 alpha^2-(R2/R0)	
	# to save myself confusion, they are altered and 
	# assigned directly to a and b for the cuadratic formula
	b <- R1/R0
	a <- .5*((b^2)-R2/R0)
	c <- -log(R0)
	r <- rep(0,2) # in case R0 = 1...
	# cuadratic formula, 2 solutions:
	r[1] <- (-b+sqrt((b^2)-4*a*c))/(2*a)
	r[2] <- (-b-sqrt((b^2)-4*a*c))/(2*a)
	# check using Lotka equation:
	nr1 <- sum(exp(-r[1]*x)*fx*Lx)
	nr2 <- sum(exp(-r[2]*x)*fx*Lx)
	# it should equal 1 exactly, but isn't quite so precise
	resid1 <- (1-nr1)^2
	resid2 <- (1-nr2)^2
	r <- ifelse(resid1<resid2,r[1],r[2])
	return(r)
}

LotkaRCoale <-
		function(fx,Lx,x){
	# from Coale, Ansley J. (1957) A New Method for Calculating Lotka's r- the Intrinsic Rate of Growth in a Stable Population.
	# Population Studies, Vol. 11 no. 1, pp 92-94
	R0 <- Rmomentn(fx,Lx,x,0)
	# first assuming a mean generation time of 29
	ri <- log(R0)/29
	for (i in 1:15){ # 10 is more than enough!
		deltai <- sum(exp(-ri*x)*fx*Lx)-1
		# the mean generation time self-corrects 
		# according to the error produced by the Lotka equation
		ri <- ri+(deltai/(29-(deltai/ri)))
	}
	return(ri)	
}

LotkaRoptim <-
		function(fx,Lx,x){
	findR <- function(r){
		# this is the Lotka equation, with the 1 moved over to the right side!
		(1-sum(exp(-r*x)*fx*Lx))^2
	}
	optimize(f=findR,seq(-1,1,.0001),tol=.00000000001)$minimum
}

MeanGeneration <-
		function(fx,Lx,x){
	r <- LotkaRCoale(fx,Lx,x)
	R0 <- Rmomentn(fx,Lx,x,0)
	R1 <- Rmomentn(fx,Lx,x,1)
	R2 <- Rmomentn(fx,Lx,x,2)
	(R1/R0)+r*.5*(((R1/R0)^2)-(R2/R0))
}

plot.R0Decomp <-
		function(x,xlim=c(0,60),col=c("#19B355","#D8FD0C","#DA1F04"),main,ylim,dx=T,leg.pos="topright"){
	if (missing(main)){
		main <- paste("R0 difference decomposition; R01 = ",round(x$R01,digits=4),", R02 = ",round(x$R02,digits=4))
	}
	# age range to be plotted
	age <- xlim[1]:xlim[2]
	
	if (dx==T){
		Epsilon <- rbind(x$Fert[(age+1)],x$SRB[(age+1)],x$Mort[(age+1)])
		leg <- c(paste("Fertility ",round(sum(x$Fert),4)),paste("SRB     ",round(sum(x$SRB),4)),paste("dx          ",round(sum(x$Mort),4)))
	}
	else {
		Epsilon <- rbind(x$Fert[(age+1)],x$SRB[(age+1)],x$Surv[(age+1)])
		leg <- c(paste("Fertility ",round(sum(x$Fert),4)),paste("SRB     ",round(sum(x$SRB),4)),paste("Survival ",round(sum(x$Surv),4)))
	}
	
	EPSpos <- Epsilon * .5*(sign(Epsilon)+1) 	
	EPSneg <- Epsilon * .5*abs(sign(Epsilon)-1) 
	
	# ylim (for looping)
	if (missing(ylim)){
		ylim <- range(pretty(c(EPSpos,EPSneg)))
	}
	ymax <- ylim[2]
	ymin <- ylim[1]
	barplot(EPSpos,space=0,col=col,ylab="contrib to diff in R0",xlab="age",main=main,ylim=ylim)
	barplot(EPSneg,add=T,space=0,col=col)
	segments(age[which((age)%%5==0)]-xlim[1],ymin,age[which((age)%%5==0)]-xlim[1],0,lty=2,col="#00000030")
	text(age[which(age%%5==0)]-xlim[1],ymin,labels=age[which(age%%5==0)],pos=1,xpd=T)
	legend(leg.pos,fill=col,legend=leg)
}

plot.R0DecompKitagawa <-
		function(x,xlim=c(10,60),col=c("#19B355","#D8FD0C","#DA1F04"),main,ylim,leg.pos="topright"){
	if (missing(main)){
		main <- paste("R0 difference decomposition; R01 = ",round(x$R01,digits=4),", R02 = ",round(x$R02,digits=4))
	}
	# age range to be plotted
	age <- xlim[1]:xlim[2]
	Epsilon <- rbind(x$Fert[(age+1)],x$SRB[(age+1)],x$Surv[(age+1)])
	leg <- c(paste("Fertility",round(sum(x$Fert),4)),paste("SRB     ",round(sum(x$SRB),4)),paste("Survival ",round(sum(x$Surv),4)))
	
	EPSpos <- Epsilon * .5*(sign(Epsilon)+1) 	
	EPSneg <- Epsilon * .5*abs(sign(Epsilon)-1) 
	
	# ylim (for looping)
	if (missing(ylim)){
		ylim <- range(pretty(c(EPSpos,EPSneg)))
	}
	ymax <- ylim[2]
	ymin <- ylim[1]
	barplot(EPSpos,space=0,col=col,ylab="contrib to diff in R0",xlab="age",main=main,ylim=ylim)
	barplot(EPSneg,add=T,space=0,col=col)
	segments(age[which((age)%%5==0)]-xlim[1],ymin,age[which((age)%%5==0)]-xlim[1],0,lty=2,col="#00000030")
	text(age[which(age%%5==0)]-xlim[1],ymin,labels=age[which(age%%5==0)],pos=1,xpd=T)
	legend(leg.pos,fill=col,legend=leg)
}

R0Decomp <-
		function(Dx1,Nx1,Bx1,px1,Dx2,Nx2,Bx2,px2){
	require(LifeTable)
	N <- length(Nx1)
	LT1 <- LT(Nx1,Dx1,axmethod="schoen",axsmooth=T,mxsmooth=F,radix=1,type="single-age")
	LT2 <- LT(Nx2,Dx2,axmethod="schoen",axsmooth=T,mxsmooth=F,radix=1,type="single-age")
	dx1 <- LT1$dx ; dx2 <- LT2$dx ;lx1 <- LT1$lx ; lx2 <- LT2$lx
	Lx1 <- (lx1[1:(N-1)]+lx1[2:N])/2 ; Lx2 <- (lx2[1:(N-1)]+lx2[2:N])/2
	fx1 <- Bx1/Nx1	;	fx2 <- Bx2/Nx2
# # # # # # # # # # # # # # # # # # # # 	
	MortComponent3 <- MortComponent2 <- MortComponent <- SurvComponent <- FertComponent <- AgeComponent <- rep(0,(N-1))
# # # # # # # # # # # # # # # # # # # # 
	R01 <- sum(px1[-N]*fx1[-N]*Lx1)
	R02 <- sum(px2[-N]*fx2[-N]*Lx2)
	R0diff <- R01-R02
# # # # # # # # # # # # # # # # # # # # 
	AgeComponent 	<- (px1[-N]*fx1[-N]*Lx1)-(px2[-N]*fx2[-N]*Lx2)
	FertComponent <- (px1[-N]+px2[-N])*(fx1[-N]-fx2[-N])*(Lx1+Lx2)*.25
	SexRatComponent <- (px1[-N]-px2[-N])*(fx1[-N]+fx2[-N])*(Lx1+Lx2)*.25
	SurvComponent 	<- (Lx1-Lx2)*(px1[-N]*fx1[-N]+px2[-N]*fx2[-N])*.5
# # # # # # # # # # # # # # # # # # # # 
	for (i in 1:(N-1)){
		MortComponent[i]	<- (dx2[i]-dx1[i])*(sum(((px1*fx1)[(i+1):N]+(px2*fx2)[(i+1):N]))+.5*((px1*fx1)[i]+(px2*fx2)[i]))*.5					
	}
	output <- list(Epsilon=AgeComponent,Fert=FertComponent,SRB=SexRatComponent,Surv=SurvComponent,Mort=MortComponent,R01=R01,R02=R02,R0diff=R0diff)
	class(output) <- "R0Decomp"
	return(output)
}

R0DecompKitagawa <-
		function(fx1,px1,Lx1,fx2,px2,Lx2){
	SurvComponent <- FertComponent <- AgeComponent <- rep(0,length(fx1))
	N <- length(fx1)
	R01 <- sum(fx1*px1*Lx1)
	R02 <- sum(fx2*px2*Lx2)
	R0diff <- R01-R02
	AgeComponent 		<- (fx1*Lx1)-(fx2*Lx2)
	FertComponent 		<- (fx1-fx2)*(Lx1+Lx2)*(px1+px2)*.25
	SexRatComponent 	<- (fx1+fx2)*(Lx1+Lx2)*(px1-px2)*.25
	SurvComponent 		<- (Lx1-Lx2)*(fx1*px1+fx2*px2)*.5
	output <- list(Epsilon=AgeComponent,Fert=FertComponent,SRB=SexRatComponent,Surv=SurvComponent,R01=R01,R02=R02,R0diff=R0diff)
	class(output) <- "R0DecompKitagawa"
	return(output)
}

Rmomentn <-
		function(fx,Lx,x,n=0){
	sum((x^n)*fx*Lx)
}

summary.R0Decomp <-
		function(x){
	R01 <- round(x$R01,4)
	R02 <- round(x$R02,4)
	R0diff <- round(x$R0diff,4)
	Epsilon <- round(R01-R02,4)
	Mortality <- round(sum(x$Mort),4)
	Fertility <- round(sum(x$Fert),4)
	SexRatio <- round(sum(x$SRB),4)
	
	results <- matrix(nrow=3,ncol=2)
	results[1,] <- c(Fertility,round((Fertility/R0diff)*100,2))
	results[2,] <- c(SexRatio,round((SexRatio/R0diff)*100,2))
	results[3,] <- c(Mortality,round((Mortality/R0diff)*100,2))
	colnames(results) <- c("absolute","percent")
	rownames(results) <- c("Fertility","SexRatio","Mortality")
	line1 <- "\nDecomposition of R0 Difference"
	line2 <- paste("\nR01 =",R01," ; R02 =",R02,"")
	line3 <- paste("\nThe difference, Epsilon, =",Epsilon,"\n\n")
	cat("\n##############################################################")
	cat(line1,line2,line3)
	print(results)
	cat("\n##############################################################")
}

summary.R0DecompKitagawa <-
		function(x){
	R01 <- round(x$R01,4)
	R02 <- round(x$R02,4)
	R0diff <- round(x$R0diff,4)
	Epsilon <- round(R01-R02,4)
	Survival <- round(sum(x$Surv),4)
	Fertility <- round(sum(x$Fert),4)
	SexRatio <- round(sum(x$SRB),4)
	
	results <- matrix(nrow=3,ncol=2)
	results[1,] <- c(Fertility,round((Fertility/R0diff)*100,2))
	results[2,] <- c(SexRatio,round((SexRatio/R0diff)*100,2))
	results[3,] <- c(Survival,round((Survival/R0diff)*100,2))
	colnames(results) <- c("absolute","percent")
	rownames(results) <- c("Fertility","SexRatio","Survival")
	line1 <- "\n3 Factor nested Kitagawa Decomposition of R0 Difference"
	line2 <- paste("\nR01 =",R01," ; R02 =",R02,"")
	line3 <- paste("\nThe difference, Epsilon, =",Epsilon,"\n\n")
	cat("\n##############################################################")
	cat(line1,line2,line3)
	print(results)
	cat("\n##############################################################")
}