Instantly share code, notes, and snippets.

# noamross/DRUG_code_LY.Rmd Created Dec 3, 2012

What would you like to do?
Matrix population models by Lauren Yamane

# A discrete time, age-structured model of a salmon population (semelparous) that can live to age 5, with fishing and environmental stochasticity (by Lauren Yamane, presented to UC Davis R Users' Group)

Parameter values

``````a=60 #alpha for Beverton-Holt stock-recruitment curve
b=0.00017 #beta for Beverton-Holt
tf=2000 #number of time steps
N0=c(100,0,0,0,0) #initial population vector for 5 age classes
s=0.28 #survival rate with fishing
e=0.1056 #fraction of age 3 fish that spawn early (age 4 is mean spawner age)
l=0.1056 #fraction of age 4 fish that spawn late as age 5 fish
sx=c(s,s,(s*(1-e)),(s*(l))) #survival vector for all ages with fishing, spawners die after spawning
t<-1 #start off at time=1
``````

Make a function for the age-structured matrix

``````#Make a function for the age-structured matrix with fishing
AgeStructMatrix_F = function(sx,a,b,tf,N0) {
sig_r=0.3
ncls=length(N0) #Number of age classes
Nt_F=matrix(0,tf,ncls) #Initialize output matrix with time steps as rows, age classes as columns
Nt_F[1,]=N0 #put initial values into first row of output matrix
for(t in 1:(tf-1)) { #for time step t in 1:1999
Pt= (e*Nt_F[t,3])+((1-l)*Nt_F[t,4])+Nt_F[t,5] #number of spawners
Nt_F[t+1,1] = ((a*Pt)/(1+(b*Pt)))*(exp(sig_r*rnorm(1,mean=0, sd=1))) #number of recruits with environmental stochasticity
Nt_F[t+1,2:ncls] = sx*Nt_F[t,1:(ncls-1)] #number of age classes 2-5
}
return(Nt_F)
}
``````

Run model by calling function

``````Nt_F=AgeStructMatrix_F(sx,a,b,tf,N0)
``````

Plot of time series with all 5 age classes

``````matplot(1:tf,Nt_F,type="l",xlab="Time",ylab="Population size", main="Age-structured model with fishing")
``````

Export text file for plot with age classes

``````colnames(Nt_F)<-c("Recruits", "Age2", "Age3", "Age4", "Age5") #Rename columns
write.table(Nt_F, "~/Desktop/Nt_F.txt", col.names=TRUE, row.names=FALSE, quote=FALSE, sep=",") #comma separated
``````

Or plot of time series for total population

``````Nt_F_totals=rowSums(Nt_F)
plot(c(1:tf),Nt_F_totals, type="l",xlab="Time",ylab="Population size", col="blue")
``````

Export text file for plot with total population

``````write.table(Nt_F_totals, "~/Desktop/Nt_F_totals.txt", col.names=FALSE, row.names=FALSE, quote=FALSE, sep=",")
``````

# Stability analysis of deterministic population model

This assumes that the partial differential equations that make up the Jacobian matrix have been determined analytically. Below are the parameters that were associated with my age-structured salmon model with fishing. I can compare the stability of the population with fishing to that with no fishing

``````e<-0.1056
l<-0.1056
s<-0.28 #survival with fishing
s_noF<-0.85 #survival with no fishing
c<-1/((e*(s^2)+(1-e)*(1-l)*(s^3)+(1-e)*l*(s^4)))
c_noF<-1/((e*(s_noF^2)+(1-e)*(1-l)*(s_noF^3)+(1-e)*l*(s_noF^4)))
alpha<-60
beta<-0.00017
``````

Fishing Jacobian

``````one_F<-c(0,0,((alpha*e)/(1+(1/(2*c))*(-c+alpha+sqrt((c-alpha)^2)))^2),((alpha*(1-l))/(1+(1/(2*c))*(-c+alpha+sqrt((c-alpha)^2)))^2),(alpha/(1+(1/(2*c))*(-c+alpha+sqrt((c-alpha)^2)))^2))
two_F<-c(s,0,0,0,0)
three_F<-c(0,s,0,0,0)
four_F<-c(0,0,((1-e)*s),0,0)
five_F<-c(0,0,0,(l*s),0)
jacobian_F<-rbind(one_F,two_F,three_F,four_F,five_F)
jacobian_F
``````

Fishing Eigenvalues

``````ev_F<-eigen(jacobian_F)
ev1_F=c(Re(ev_F\$values[1]),Im(ev_F\$values[1]))
ev2_F=c(Re(ev_F\$values[2]),Im(ev_F\$values[2]))
ev3_F=c(Re(ev_F\$values[3]),Im(ev_F\$values[3]))
ev4_F=c(Re(ev_F\$values[4]),Im(ev_F\$values[4]))
ev5_F=c(Re(ev_F\$values[5]),Im(ev_F\$values[5]))
``````

No Fishing Jacobian

``````one_noF<-c(0,0,((alpha*e)/(1+(1/(2*c_noF))*(-c_noF+alpha+sqrt((c_noF-alpha)^2)))^2),((alpha*(1-l))/(1+(1/(2*c_noF))*(-c_noF+alpha+sqrt((c_noF-alpha)^2)))^2),(alpha/(1+(1/(2*c_noF))*(-c_noF+alpha+sqrt((c_noF-alpha)^2)))^2))
two_noF<-c(s_noF,0,0,0,0)
three_noF<-c(0,s_noF,0,0,0)
four_noF<-c(0,0,((1-e)*s_noF),0,0)
five_noF<-c(0,0,0,(l*s_noF),0)
jacobian_noF<-rbind(one_noF,two_noF,three_noF,four_noF,five_noF)
jacobian_noF
``````

No Fishing Eigenvalues

``````ev_noF<-eigen(jacobian_noF)
ev1_noF=c(Re(ev_noF\$values[1]),Im(ev_noF\$values[1]))
ev2_noF=c(Re(ev_noF\$values[2]),Im(ev_noF\$values[2]))
ev3_noF=c(Re(ev_noF\$values[3]),Im(ev_noF\$values[3]))
ev4_noF=c(Re(ev_noF\$values[4]),Im(ev_noF\$values[4]))
ev5_noF=c(Re(ev_noF\$values[5]),Im(ev_noF\$values[5]))

# Setup for plot, just remember to expand the plot interface first
``````

Plot the eigenvalues with legend to compare stability of both deterministic models

``````layout(matrix(c(1,0,2,0,3,0), 2,3, byrow=TRUE)) #Layout with 2 rows, 3 columns,
#the matrix tells you the number of the figure in each location.
#Since the layout is byrow, it fills the first row first and then the 2nd row
#no fishing plot
plot(ev1_noF[1],ev1_noF[2], type="p", pch=0, cex=2, xlab="Real", ylab="Imaginary", xlim=c(-1.25,1.25),ylim=c(-1.25,1.25))
title("No fishing", adj=0) #plot dominant eigenvalue, with x-coordinate as real part and y-coordinate as imaginary part
draw.circle(0,0,radius=1, border="black", lty=1, lwd=1) #draw unit circle
points(ev2_noF[1],ev2_noF[2], pch=1, cex=2) #2nd eigenvalue
points(ev3_noF[1],ev3_noF[2], pch=2, cex=2) #3rd eigenvalue
points(ev4_noF[1],ev4_noF[2], pch=6, cex=2) #4th eigenvalue
points(ev5_noF[1],ev5_noF[2], pch=5, cex=2) #5th eigenvalue
abline(a=NULL, b=NULL, h=0, lty=3)
abline(a=NULL, b=NULL, v=0, lty=3)

#fishing plot
plot(ev1_F[1],ev1_F[2], type="p", pch=0, cex=2, xlab="Real", ylab="Imaginary", xlim=c(-1.25,1.25),ylim=c(-1.25,1.25))
points(ev2_F[1],ev2_F[2], pch=1, cex=2)
points(ev3_F[1],ev3_F[2], pch=2, cex=2)
points(ev4_F[1],ev4_F[2], pch=6, cex=2)
points(ev5_F[1],ev5_F[2], pch=5, cex=2)
abline(a=NULL, b=NULL, h=0, lty=3)
abline(a=NULL, b=NULL, v=0, lty=3)

#Legend
plot(1:10, type="n", axes=FALSE, xlab="", ylab="", frame.plot=TRUE)