{VERSION 5 0 "SUN SPARC SOLARIS" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 2 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 1 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Headi ng 1" -1 3 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 1 1 2 2 2 2 1 1 1 1 }1 1 0 0 6 6 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 0" -1 256 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 2" -1 257 1 {CSTYLE "" -1 -1 " Courier" 1 14 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 2" -1 258 1 {CSTYLE "" -1 -1 "Courier" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2 " -1 259 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 1 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 " " }}}{EXCHG {PARA 0 "" 0 "" {TEXT 256 10 "Leslie.mws" }{TEXT -1 57 " (revised Sept. 8, \+ 2001) Maple worksheet to accompany " }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 28 "Brault and Caswell, 1993, " }{TEXT 262 40 "Pod-specific demography of killer-whales" }{TEXT -1 10 ", Ecol ogy " }{TEXT 261 2 "74" }{TEXT -1 12 ": 1444-1454 " }}{PARA 0 "" 0 "" {TEXT -1 4 "and " }}{PARA 0 "" 0 "" {TEXT -1 21 "Rodrigo Bernal,1998. \+ " }{TEXT 263 40 "Demography of the vegetable ivory palm " }{TEXT 264 19 "Pytelephas seemanii" }{TEXT 265 48 " in Colombia, and the impact o f seed harvesting." }{TEXT -1 28 " Journal of Applied Ecology " } {TEXT 260 2 "35" }{TEXT -1 7 ":64-74 " }}{PARA 0 "" 0 "" {TEXT -1 34 " and other Leslie-Lefkowitch models" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 329 "Stage-based models are a modification of the Leslie model. In this case, the assumption is that the stage of \+ an organism (size class) is a better indicator of reproductive perform ance and survival than age is. This assumption works well for many pl ants, and any organism in which size is correlated with survival and f ecundity." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 304 "The general form of the model treats the stages as elements of a \+ vector. The population projection matrix A (sometimes called a Lefkov ich matrix, or, among mathematicians, transition matrix) is used to ca lculate the number of individuals in each stage class at time t+1, bas ed on the numbers at time t. " }}{PARA 0 "" 0 "" {TEXT -1 21 " \+ " }}{PARA 0 "" 0 "" {TEXT -1 10 " " }{TEXT 257 2 "n(" }{TEXT -1 3 "t+1" }{TEXT 258 9 ") = A n(" }{TEXT -1 1 "t" } {TEXT 259 2 ") " }{TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 6 " \+ " }}{PARA 0 "" 0 "" {TEXT -1 44 "The (non-zero) elements of the matri x A are:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 68 "Pi the probabilities of surviving and remaining in the same stag e." }}{PARA 0 "" 0 "" {TEXT -1 63 "Gi the probabilities of surviving a nd moving to the next stage." }}{PARA 0 "" 0 "" {TEXT -1 90 "Fi the f ertility, or number of female offspring at time t+1, per adult female at time t." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "restart: wi th(linalg): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "A:= matrix (4,4, [ [0,F2,F3,F4], [G1,P2,0,0], [0,G2,P3,0], [0,0,G3,P4] ]);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 102 "Here are the respective matrices: Caswell goes with the orca (whale) model and Bernal with palm trees. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "Caswell:= matrix(4, 4, [ [0,0.0043,0.1132,0], [0.9775,0.9111,0,0],\n[0,0.0736,0.9534,0], [0,0 ,0.0452,0.9804] ] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 185 "Be rnal:= matrix(6,6,[[0.7052,0,0,14.8173,14.9770,14.7732],\n[0.1548,0.63 39,0,0,0,0],\n[0,0.0216,0.9100,0,0,0],\n[0,0,0.0330,0.9614,0,0],\n[0,0 ,0,0.0386,0.9220,0],\n[0,0,0,0,0.0188,0.9535]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 30 "Bern al stable age distribution" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "A:= evalm( Bernal ); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "n:= rowdim(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 257 " The elements on the top row are the reproductiv e outputs. The subdiagonal elements represent the probability of tran sition to the next stage. The diagonal elements represent the probabil ity of remaining in a stage. here are some questions to think about. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 "*In \+ the palm tree model what fraction of seeds remain in the 'seed bank' e ach year, that is, do not germinate? " }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 182 "**Of 1000 stage 1 individuals, how m any survive to become stage 2? If there are currently 1000 stage 2 an d 1000 stage 3, how many stage 1 individuals will be present after 1 y ear? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 143 "***Suppose an individual is in the highest stage. What is the pr obablility that this individual will still be alive 14 time steps from now?***" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 176 " Let's track what happens to a single burst of reproduct ion, say 1000 seeds or newborns. Each time step is computed by taking the matrix A times the current population vector." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "pop||0:= vector( n \+ , [1000, 0, 0, 0, 0, 0]); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "distrib= vector( n , [ 1, 0, 0, 0, 0, 0] ); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digits:=6; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "points:= [0, 1000]; lpoints:= [0, 3];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 142 "Next we' re going to track the population as it changes. We'll go for 50 years , but take a snapshot every 5 years. Feel free to run it longer." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for i from 0 by 5 to 50 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "pop||(i+5):= evalm ( A^5 &* pop||i ): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "total||(i+5):= sum( (pop||(i+5))['j'] , 'j' = 1 .. n) :" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "points:= points, [i+5, total|| (i+5)]: lpoints:= lpoints, [i+5, log10(total||(i+5))]:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 54 "distrib:= scalarmul( pop||(i+5), 1/ total||(i+ 5)): od;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 266 "What app ears to be happening? Is the population growing? What's happening in the individual stages? to the overall distribution? Try it over again with more than 50 years. If you get too much output change od; to od : and use the following line for the final step." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 150 "pop||i:= evalm(pop||i); total||i:= total||i ; distrib:= evalm(distrib); five_year_rate:= total||i/total||(i-5); a nnual_rate:= (five_year_rate)^(1/5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "plot([points], title = \"Total population as a functi on of time\" );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "plot([lpoints], title = \"Lo g10 of total pop. vs. time\");" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "Is there any pattern h ere? What about in the long term?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 66 "Now let's compare this with the eigenvalue - eigen vector analysis." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys: =eigenvects(A): # (output printing suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:= seq( eigsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors (output printing suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys[i] [2], i = 1 .. nops([eigsys])): # multiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "magnitudes:=op( map( abs, [EVal] )); # magnitudes of \+ eigenvalues" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 195 "Out of the comple te list of eigenvalues we find the maximum or \"dominant\" eigenvalue \+ (lambda = population growth rate), and the index to the corresponding \+ eigenvector (stable age distribution). " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "lambda:=max(magnitudes);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "member( lambda, [magnitudes], 'position'): indx:= po sition: # isolate the dominant eigenvalue" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "V:=op( EVec[indx] ); # find an eigenvector that is as sociated with this eigenvalue" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 66 "How does this compare with the annual rate that you found earlier?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "annual_rate;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 507 "If the population grows for many \+ generations with constant transition probabilities given by matrix A, \+ the population age or stage vector will approach an eigenvector associ ated with the dominant eigenvector of A. This is the stable age or s tage distribution vector. Typically the stable stage vector w is scal ed so that the sum of the elements is 1.0. To scale the stable stage \+ vector, we have to do some fiddling. Add up the elements in the vecto r; then divide all elements in the vector by this sum:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 45 "total:= sum( V['j'] , 'j' = 1 .. rowdim(A) ); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "We now have the stable age or stage distribution vector:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "w:=scalarmul(V, 1/total);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "How does this compare with the distribution vector that you fou nd iteratively?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "evalm(di strib);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 34 "Let's check; is A w = (lambda) w? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 337 " Remember this is our de finition of an eigenvalue and eigenvector - that the eigenvector multi plied by the matrix is the same as the eigenvector multiplied by the e igenvalue. If we got the correct eigenvalues and eigenvectors, we sho uld obtain the same answer by multiplying the stable age distribution \+ by matrix A or by the eigenvalue:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "evalm( A &* w ); scalarmul( w , lambda);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 31 "Caswel l stable age distribution" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "restart: with(linalg): with(plots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 110 "Caswell:=matrix(4, 4,[ [0,0.0043,0.1132,0], [0.9775, 0.9111,0,0],\n[0,0.0736,0.9534,0], [0,0,0.0452,0.9804] ] );" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "A:= evalm( Caswell ); " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "n:= rowdim(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 257 " The ele ments on the top row are the reproductive outputs. The subdiagonal el ements represent the probability of transition to the next stage. The \+ diagonal elements represent the probability of remaining in a stage. \+ here are some questions to think about." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 183 "*In the whale model what tells you \+ that juveniles do not reproduce? that older females do not reproduce? \+ Does this remind you of any other creature with a similar life stage p attern?*" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 182 "**Of 1000 stage 1 individuals, how many survive to become stage 2 ? If there are currently 1000 stage 2 and 1000 stage 3, how many stag e 1 individuals will be present after 1 year? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 143 "***Suppose an individual is in the highest stage. What is the probablility that this individu al will still be alive 14 time steps from now?***" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 176 " Let's track wha t happens to a single burst of reproduction, say 1000 seeds or newborn s. Each time step is computed by taking the matrix A times the curren t population vector." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "pop||0:= vector( n , [1000, 0, 0, 0 ]); " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "distrib= vector( n , [ 1, 0, 0, 0] ); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Digits:=6; \+ " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "points:= [0, 1000]; lpo ints:= [0, 3];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for i fro m 0 by 10 to 100 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "pop||(i+10): = evalm ( A^10 &* pop||i ):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "tota l||(i+10):= sum( pop||(i+10)['j'] , 'j' = 1 .. n):" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 54 "distrib:= scalarmul( pop||(i+10), 1/ total||(i+10) \+ ): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "points:= points, [(i+10), to tal||(i+10)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "lpoints:= lpoints, [(i+10), log10(total||(i+10))]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " od;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 383 "What appears to be happening? Is the population growing? Wha t's happening in the individual stages? to the overall distribution? T ry it over again with more than 20 steps. Suggestion: change the ; to \+ a : after the od to avoid screenfulls of output, and don't be too gree dy. Let's just print the last computation and get an idea of the grow th rate by comparing the last two totals." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 146 "pop||i:= evalm(pop||i); total||i:= total||i; dis trib:= evalm(distrib); ten_yr_rate:= total||i/total||(i-10); annual_ra te:= (ten_yr_rate)^(1/10);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "plot( [points] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "plot( [lpoints] );" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "Now let's compare this wi th the eigenvalue - eigenvector analysis." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output printing suppress ed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:= seq( eigsys[i][3], i \+ = 1 .. nops([eigsys])): # eigenvectors (output printing suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys[i][2], i = 1 .. nops([eigsys])): # m ultiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "magnitudes:=op( ma p( abs, [EVal] )); # magnitudes of eigenvalues" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 261 "Out of the complete list of eigenvalues we find the maximum or \"dominant\" eigenvalue (lambda = population growth rate), and the index to the corresponding eigenvector (stable age distributi on). How are these eigenvalues different from the ones we got with th e " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "lambda:=max(magnitude s);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "member( lambda, [magnitudes] , 'position'): indx:= position:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 508 "If the population grows for many generat ions with constant transition probabilities given by matrix A, the pop ulation age or stage vector will approach the eigenvector associated w ith the dominant eigenvector of A. This is the stable age or stage d istribution vector. Typically the stable stage vector w is scaled so \+ that the sum of the elements is 1.0. To scale the stable stage vector , we have to do some fiddling. Add up the elements in the vector; the n divide all elements in the vector by this sum:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "total:= sum( V['j'] , 'j' = 1 .. n );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "We now ha ve the stable age or stage distribution vector:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "w:=scalarmul(V, 1/total);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 89 "Compare wit h the annual growth rate and the long term distribution that we found \+ earlier." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "annual_rate; e val(distrib);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 34 "Let's check; is A w = (lambda) w? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 337 " Remember this is our de finition of an eigenvalue and eigenvector - that the eigenvector multi plied by the matrix is the same as the eigenvector multiplied by the e igenvalue. If we got the correct eigenvalues and eigenvectors, we sho uld obtain the same answer by multiplying the stable age distribution \+ by matrix A or by the eigenvalue:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "evalm( A &* w ); scalarmul( w , lambda);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 259 "" 0 "" {TEXT -1 46 "\"Fi sh model\" and a look at reproductive value " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "restart: with(linalg): with(plots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "Fish:= matrix( 3, 3, [ [0, 0.8, 2.2 ], [0.5, 0, 0], [0, 0.6, 0] ]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "A:= evalm( Fish ); n:= rowdim(A);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 35 "pop||0:= vector( n , [100, 0, 0]); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "distrib= vector( n , [ 1, 0, 0] ); \+ " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "total||0:= 100; " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "points:= [0, 100]; lpoints:= [0, 2];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for i from 0 by 1 to 30 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 35 "pop||(i+1):= evalm ( A &* pop||i ):" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 52 "total||(i+1):= sum( pop||(i+1)['j'] , 'j' = 1 \+ .. n):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "distrib:= scalarmul( pop| |(i+1), 1/ total||(i+1) ): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "poin ts:= points, [(i+1), total||(i+1)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "lpoints:= lpoints, [(i+1), log10(total||(i+1))]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "pop||i:= evalm(pop||i) ; total||i:= total||i; distrib:= evalm(distrib); annual_rate:= (tota l||i)/(total||(i-1));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "plot( [points] );" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 18 "plot( [lpoints] );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "Here's the stable age distribution anaysis." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output p rinting suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:= seq( \+ eigsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors (output printi ng suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsy s[i][1], i = 1 .. nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys[i][2], i = 1 .. nops([eigsys])): # m ultiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "magnitudes:=op( ma p( abs, [EVal] )); # magnitudes of eigenvalues" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 24 "lambda:=max(magnitudes);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "member( lambda, [magnitudes], 'position'): indx:= po sition:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "total:= sum( V['j'] , 'j' = 1 .. n \+ ):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "We now have the stable age or stage distribution vector:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "w:=scalarmul(V, 1/total);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "Compare with our 30 year iteration:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "annual_rate; evalm(distrib);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 155 "Here we transpose the A matrix in order to set up calculation of the reproductive value vector, whic h is the dominant eigenvector of the transposed matrix." }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 30 "AT:=transpose(A); Digits:= 10:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing \+ suppressed)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 120 "Find maximum eigenvalue and index to corresponding eig envector; this eigenvector is the reproductive value distribution." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "EVecT:= seq( eigsysT[i][3], i = 1 . . nops([eigsysT]) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) \+ ); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EMT:= seq( eigs ysT[i][2], i = 1 ..nops([eigsysT]) ): # multiplicities" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 137 "magsT:= op( map( abs , [EValT] )); # magnitud es of eigenvalues (Note: [ ] turns a sequence into a list; op turns a \+ list into a sequence.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "m u:= max( magsT); # get dominant eigenvalue" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "member(mu, [magsT], 'position'): indexT:= position: " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 "Convert eigenvector lists to actual vectors for use in sensitivity analysis." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "VT:= op( EVecT[indexT ] ); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 309 "This is the reproductive value distribution. Note that \+ vector v is scaled so the first element is 1.0. We do a scalar multi plication of all elements in the vector by the reciprocal of the first element. (Recall that a scalar multiple of an eigenvector is again a n eigenvector with the same eigenvalue.) " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "v:= scalarmul(VT, 1/VT[1] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 180 "What in formation does this vector give us? Which age or stage group is contri buting the most to reproduction, and why is this plausible in terms of both fecundity and survivorship?" }}{PARA 0 "" 0 "" {TEXT -1 1 "." }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "Let's confirm that we correctly \+ identified the eigenvalue/eigenvector pair by checking that AT v = (mu ) v. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "evalm( AT &* v ); scala rmul( v , mu);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}} {SECT 1 {PARA 3 "" 0 "" {TEXT -1 43 "Sensitivity analysis for Bernal \+ (optional)" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 831 "Now we will calculate the sensitivity of the system to c hanges in fecundity. This is the analysis proposed in Caswell. We ar e analytically determining the partial derivative of lambda (populatio n growth rate) with respect to each element in the projection matrix. \+ Thus we are asking: given a change in a particular element in the pro jection matrix, what is the corresponding change in population growth \+ rate? This allows us to determine which components of the life histor y have the biggest effects on population growth, and which may be most sensitive to natural selection. The first parameter we need is the i nner product of the two eigenvectors (stable age distribution and repr oductive value distribution). This is symbolized as in the Ca swell handout. This value is used as the denominator of the calculati on. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 185 "Bernal:= matrix(6,6,[[0.7052,0,0,1 4.8173,14.9770,14.7732],\n[0.1548,0.6339,0,0,0,0],\n[0,0.0216,0.9100,0 ,0,0],\n[0,0,0.0330,0.9614,0,0],\n[0,0,0,0.0386,0.9220,0],\n[0,0,0,0,0 .0188,0.9535]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A:=eval m(Bernal):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigen vects(A): # (output printing suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:= seq( eigsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors (output printing suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys])) ; # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys [i][2], i = 1 .. nops([eigsys])): # multiplicities" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 64 "magnitudes:=op( map( abs, [EVal] )); # magnitudes o f eigenvalues" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "lambda:=max(magnitudes);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "m ember( lambda, [magnitudes], 'position'): indx:= position:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 45 "total:= sum( V['j'] , 'j' = 1 .. rowdim(A) ) ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Calculate the scaled stable \+ age distribution:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "w:=scalarmul(V ,1/total);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 145 "Transpose the matr ix in order to set up calculation of the reproductive value vector, wh ich is the dominant eigenvector of the transposed matrix." }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing suppres sed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 120 "Find maximum eigenvalue and index to corresponding eigen vector; this eigenvector is the reproductive value distribution." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "EVecT:= seq( eigsysT[i][3], i = 1 . . nops([eigsysT]) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) \+ ); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EMT:= seq( eigs ysT[i][2], i = 1 ..nops([eigsysT]) ): # multiplicities" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 137 "magsT:= op( map( abs , [EValT] )); # magnitud es of eigenvalues (Note: [ ] turns a sequence into a list; op turns a \+ list into a sequence.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "m u:= max( magsT); # get dominant eigenvalue" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "member(mu, [magsT], 'position'): indexT:= position :" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 114 "Convert eigenvector lists to actual vectors for use in sensitivit y analysis. VT is the reproductive value vector." }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 26 "VT:= op( EVecT[indexT] ); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "v is the scaled reproductive value vector (scaled so t he first element is 1.0)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "v:=scal armul(VT,1/VT[1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 114 "Sensitivit y to fecundity is related both to the stable age distribution and to t he reproductive value distribution" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "IP:= innerprod(VT , V); # dot, or scalar, product \+ " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "sensitivity_to_fecundity:= [seq ( V[i]*VT[1] / IP , i = 1 .. rowdim(A) )]; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "Observe that it is immaterial whether we use the normalized vectors or not, as lon g as we do this consistently." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "[seq( w[i]*v[1] /innerprod( w, v) , i = 1 .. rowdim(A) )]; " }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 215 "In terms of the paper this corres ponds to the top row of matrix S, equation 7, in Caswell except that t he first and last elements have been suppressed because the original m atrix had zero entries in these locations." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "We do a similar calcuation to o btain the sensitivity to survival: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "sensitivity_to_survival:= [seq( V[i]*VT[i+1] / IP , i = 1 .. row dim(A)-1 )];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 401 "This vector corr esponds to the subdiagonal of matrix S. The entire sensitivity matri x can be calcuated by generalizing the above operations. Note there a re some entries that can be taken to be zero, or otherwise ignored, be cause the A matrix had zero entries (but this isn't really necessary, \+ since as we shall see in a moment we end up multiplying each entry of \+ S by the corresponding entry of A)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 106 "S:=matrix(rowdim(A),rowdim(A), [ seq( seq( V[i] * VT[j] / IP, i = 1.. rowdim(A)), j = 1.. rowdim(A) ) ]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 418 "The elasticity values are calculated by scaling the sens itivities by aij / lambda (Caswell handout). The resulting matrix is \+ reported as Table 4 on p 70. The elasticities are a bit easier to int erpret than the sensitivities, because they represent the proportional contributions of matrix elements to lambda. Recall the we can't use \+ the label E in Maple because this means the base of natural logarithms (2.718...)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "e:=matrix(rowdim(A ),rowdim(A),[ seq( seq( S[j,i] * A[j,i] / lambda, i = 1 .. rowdim(A)), j = 1.. rowdim(A)\\)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 401 "***Try seeing if you believe that the se nsitivities or elasticities are correct, following the pattern illustr ated below.*** We will change one of the entries in the original \"A \" matrix by 10% and see what percent change there is in lambda. Comp are this change with the corresponding value in the elasticity or sens itivity matrix. Let's first recall the original matrix, eigenvalues, eigenvectors:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "A:= evalm(A);" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "lambda[old]:= lambda; mu[ old]:= mu;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "RV:= evalm(v) ; SAD:= evalm(w);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 428 "Let's chan ge the P2 value by10% (we'll do this by multiplying P2 by 1.1) and see what percent change we get in lambda. This is the so-called [2,2] en try of the matrix A, that is, the number that appears in the first row (from top to bottom) and second column (from left to right). If the author's description of elasticity is correct, then we should get a c hange in lambda that is proportional to the elasticity value for P2." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "A[2,2]:=1.1*A[2,2]; # make the \+ change in the affected entry" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "A:=evalm(A); # confirm the matrix is OK" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "Now let's calculate \+ the eigenvalues and eigenvectors for this revised matrix." }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output printing s uppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "EVec:= seq( e igsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. \+ nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "m ags:= op( map( abs, [EVal] )); lambda[new]:=max(mags);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 63 "Compare t his new value for lambda with the old value of lambda:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "lambda[old];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 163 "This is the ratio of the new value to the old; the \+ fractional change in lambda is the portion after the decimal. Did you expect lthe new lambda to be larger? Why?" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 57 "ratio:= lambda[new]/lambda[old]; frac_change:= rati o - 1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 "Here is the elasticity \+ value. Compare it to the fractional change in lambda." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "elas:= e[2,2];" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 410 "Note that we made a 10 percent change in one elem ent of the matrix (a fractional change of 0.1). If you multiply this \+ elasticity value by the fractional change made in the matrix element, \+ you should get a number close to the fractional change in lambda. If you made a 100% change (a fractional change of 1.0) in the matrix ele ment, you should see a fractional change in lambda equal to the elasti city value." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "(elas * 0.1) - (frac_change); # close to zero?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "member( lambda[new], [mags], 'position'): indx:= pos ition:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "evalm(SAD); # here's the \+ original stable age distribution" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "\nCompare with the new stable age distribution. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 120 "V := op( EVec[indx] ); total:= sum( V['k '] ,' k' = 1 .. rowdim(A) ): w:= scalarmul(V, 1/total); # compare with SAD above" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 147 "Transpose the A ma trix in order to set up calculation of the reproductive value vector, \+ which is the dominant eigenvector of the transposed matrix." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing suppres sed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "EVecT:= seq( eigsys T[i][3], i = 1 .. rowdim(A) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "EValT:= seq( eigsysT[i][1], i = 1 .. rowdim(A) ) ; # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "magsT:= op( map( abs , [EValT] )); # magnitudes of eigenvalues " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "mu[new]:= max( magsT); member(mu[new], [magsT], 'po sition'): indexT:= position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "VT:= op( EVecT[indexT] ): v:= scalarmul(VT, 1/VT[1] ); # compare with RV above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "evalm(RV) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 365 "***Now, how did v and w change, given that 10% chan ge in a single entry of the matrix A? Considering the change that was made can you give an explanation for the result? Now change the surv ival probability G3 by 15% (multiply it by 1.5). (Don't forget to res et A[1,2] back to its original value first) What effect do you expect to see? What happens in fact?***" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 252 "The author states on p71 that \"an 86% r eduction in fecundity of all classes is required to reduce lambda from 1.0589 to 1.0\". Try making reductions in fecundity and see if your \+ observations confirm this: First we will get an original copy of matr ix A" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A:=evalm(Bernal);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "Then we define our fractional reduction in fecundity." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "fecundity_reduction:=0.86;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 271 "We will now reduce the seed pr oduction entries by this amount. We will not change A[1,1], because t hat element of the matrix defines the rate of retention of seeds in th e seed bank, not a fecundity. Therefore we change only elements 2 to \+ 6 in the top row of the matrix:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for i from 2 to rowdim(A) do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 " A[1,i]:=A[1,i]*(1-fecundity_reduction); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "evalm(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 235 "Now move back to the top of this section of the workshee t, and rerun from just below the original definition of matrix A. See if the values for lambda, stable stage distribution, and reproductive value distribution change as expected." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 44 "Sensitivity analysis for Caswell (optional)" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 99 "FIrst we need to get the eigenvalues and eigenvectors of both the \+ original matrix and its transpose" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "restart: with(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "Caswell:= matrix(4, 4,[ [0,0.0043,0.1132,0], [0.9775 ,0.9111,0,0],\n[0,0.0736,0.9534,0], [0,0,0.0452,0.9804] ] );" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "A:= evalm(Caswell):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (o utput printing suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec: = seq( eigsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors (output printing suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EV al:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); # eigenvalues" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys[i][2], i = 1 .. no ps([eigsys])): # multiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 " magnitudes:=op( map( abs, [EVal] )); # magnitudes of eigenvalues" } {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "lambda:=max (magnitudes);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "member( lambda, [m agnitudes], 'position'): indx:= position:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "total:= sum( V['j'] , 'j' = 1 .. rowdim(A) );" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Calculate the scaled stable age di stribution:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "w:=scalarmul(V,1/tot al);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 145 "Transpose the matrix in \+ order to set up calculation of the reproductive value vector, which is the dominant eigenvector of the transposed matrix." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing suppressed)" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 120 "Find maximum eigenvalue and index to corresponding eigenvector ; this eigenvector is the reproductive value distribution." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "EVecT:= seq( eigsysT[i][3], i = 1 .. nops([ eigsysT]) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) ); # eigenval ues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EMT:= seq( eigsysT[i][2], i \+ = 1 ..nops([eigsysT]) ): # multiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 137 "magsT:= op( map( abs , [EValT] )); # magnitudes of e igenvalues (Note: [ ] turns a sequence into a list; op turns a list in to a sequence.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "mu:= max ( magsT); # get dominant eigenvalue" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "member(mu, [magsT], 'position'): indexT:= position:" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 114 "C onvert eigenvector lists to actual vectors for use in sensitivity anal ysis. VT is the reproductive value vector." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "VT:= op( EVecT[indexT] ); " }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 78 "v is the scaled reproductive value vector (scaled so th e first element is 1.0)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "v:=scala rmul(VT,1/VT[1]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 940 "N ow we will calculate the sensitivity of the system to changes in fecun dity. This is the analysis proposed on p. 1447 in equation 3. We are analytically determining the partial derivative of lambda (population growth rate) with respect to each element in the projection matrix. \+ Thus we are asking: given a change in a particular element in the proj ection matrix, what is the corresponding change in population growth r ate? This allows us to determine which components of the life history have the biggest effects on population growth, and which may be most \+ sensitive to natural selection. The first parameter we need is the in ner product of the two eigenvectors (stable age distribution and repro ductive value distribution). This is symbolized as in equatio n 3 on p. 1447. This value is used as the denominator on the right ha nd size of equation 3. Observe that it is immaterial whether we use t he normalized vecotrs or not." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "IP:= innerprod(VT , V); # dot, or scalar, product " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "sensitivity_to_fecundity:= [seq( V[i]*VT[ 1] / IP , i = 1 .. 4 )]; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 111 "Observe that it is immaterial whether we use the normalized vectors or not, as long as we do this consistently ." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "[seq( w[i]*v[1] /inner prod( w, v) , i = 1 .. 4 )]; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 215 "In terms of the paper this correspon ds to the top row of matrix S on p. 1448, equation 7, except that the \+ first and last elements have been suppressed because the original matr ix had zero entries in these locations." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "We do a similar calcuation to obtain the sensitivity to survival: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 " sensitivity_to_survival:= [seq( V[i]*VT[i+1] / IP , i = 1 .. 3 )];" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 401 " This vector corresponds to the subdiagonal of matrix S. The entire s ensitivity matrix can be calcuated by generalizing the above operation s. Note there are some entries that can be taken to be zero, or other wise ignored, because the A matrix had zero entries (but this isn't re ally necessary, since as we shall see in a moment we end up multiplyin g each entry of S by the corresponding entry of A)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "S:=matrix(4,4, [ seq( seq( V[i] * VT[j] / IP, i = \+ 1.. 4), j = 1.. 4 ) ]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 420 "The elasticity values are calculated by \+ scaling the sensitivities by aij / lambda (equation 4 on p. 1447). Th e resulting matrix is reported as equation 8 on p. 1448. The elastici ties are a bit easier to interpret than the sensitivities, because the y represent the proportional contributions of matrix elements to lambd a. Recall the we can't use E in Maple because this means the base of \+ natural logarithms (2.718...)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "e :=matrix(4,4,[ seq( seq( S[j,i] * A[j,i] / lambda, i = 1 .. 4), j = 1. . 4)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 401 "***Try seeing if you believe that the sensitivities or e lasticities are correct, following the pattern illustrated below.*** \+ We will change one of the entries in the original \"A\" matrix by 10% \+ and see what percent change there is in lambda. Compare this change w ith the corresponding value in the elasticity or sensitivity matrix. \+ Let's first recall the original matrix, eigenvalues, eigenvectors:" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "A:= evalm(A);" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 37 "lambda[old]:= lambda; mu[old]:= mu;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "RV:= evalm(v); SAD:= evalm( w);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 416 "Let's change the F2 value from 0.0043 to 0.00473 (a 10% change ) and see what percent change we get in lambda. This is the so-called [1,2] entry of the matrix A, that is, the number that appears in the \+ first row (from top to bottom) and second column (from left to right). If Caswell's description of elasticity is correct, then we should g et a change in lambda that is proportional to the elasticity value for F2." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "A[1,2]:=.00473; # make th e change in the affected entry" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "A:=evalm(A); # confirm the matrix is OK" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "Now let's calculat e the eigenvalues and eigenvectors for this revised matrix." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output printin g suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "EVec:= seq ( eigsys[i][3], i = 1 .. 4): # eigenvectors " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "EVal:= seq( eigsys[i][1], i = 1 .. 4); # eigenva lues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "mags:= op( map( abs, [EVal] )); lambda[new]:=max(mags);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 118 "Compare this new value for lambda wi th the \"estimate\" of lambda at bottom of p 1447. This is the old va lue of lambda:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "lambda[ol d];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 163 "This is the ratio of the \+ new value to the old; the fractional change in lambda is the portion a fter the decimal. Did you expect lthe new lambda to be larger? Why?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "ratio:= lambda[new]/lambd a[old]; frac_change:= ratio - 1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 "Here is the elasticity value. Compare it to the fractional change in lambda." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "elas:= e[1,2 ];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 390 "Note that we made a 10 per cent change in fecundity (a fractional change of 0.1). If you multipl y this elasticity value by the fractional change made in the fecundity , you should get a number close to the fractional change in lambda. \+ If you made a 100% change (a fractional change of 1.0) in the fecundit y value, you should see a fractional change in lambda equal to the ela sticity value." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "(elas * 0 .1) - (frac_change); # close to zero?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "member( lambda[new], [mags], 'position'): location:= position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "evalm(SAD); # here's the original stable age distribution" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 47 "\nCompare with the new stable age distribution. " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 116 "V := op( EVec[location] ); \+ total:= sum( V['k'] ,' k' = 1 .. 4 ):\nw:= scalarmul(V, 1/total); # co mpare with SAD above" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 147 "Transpos e the A matrix in order to set up calculation of the reproductive valu e vector, which is the dominant eigenvector of the transposed matrix. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "EVecT:= seq ( eigsysT[i][3], i = 1 .. 4 ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "EValT:= seq( eigsysT[i][1], i = 1 .. 4 ); # eige nvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "magsT:= op( map( abs , [ EValT] )); # magnitudes of eigenvalues " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "mu[new]:= max( magsT); member(mu[new], [magsT], 'position'): indexT:= position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "VT: = op( EVecT[indexT] ): v:= scalarmul(VT, 1/VT[1] ); # compare with RV \+ above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "evalm(RV);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 372 "***Now, ho w did v and w change, given that 10% change in a single entry of the m atrix A? Considering the change that was made can you give an explana tion for the result? Now change the survival probability G3 by 15% fr om 0.0452 to 0.05198. (Don't forget to reset A[1,2] back to its origi nal value of 0.0043!) What effect do you expect to see? What happens in fact?***" }}}}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }