“First Law” of Population Dynamics (like Newton’s First ...



Bio 292: Population Ecology Fall, 2011

Instructor: Bill Morris

MW 2:50-4:05pm, 130 BioSci

Tentative syllabus:

Week: Topic: Readings in

Morris & Doak book:

1-2 Density-independent models based on the total number Ch. 1-3

of individuals in a population under a randomly

varying environment;

Introduction to programming in R

3-4 More complex models based on total population size: Ch. 4

negative and positive density dependence, environmental

autocorrelation, and catastrophes/bonanzas

5-6 Deterministic projection matrix models for age- or Ch. 6, 7

size-structured populations

7 Sensitivity analysis for deterministic projection matrix Ch. 9 to p. 351

models

8-9 Stochastic projection matrix models; Ch. 8

Additional complications in matrix models

10 Meta-population and other spatial models Ch. 10-11

11-12 Models of species interactions (if there is time)

13 Students present results of research projects

Course requirements:

1) Do problem sets and worksheets (not graded; only credit/no credit)

2) Class project:

- create and analyze a population model using your own data or data from the literature

- report results to the class in the last week of the semester

Day 1: Intro and Understanding the effects of variability

Introduction and Goals of the Course:

Definition of a population: the set of individuals of a single species in a defined area

Goals of Population Ecology

To find answers to questions such as:

- Why are there so many (or few) of a given species in a particular place?

- Why do numbers change (or not) over time or space?

- Why do we see the observed ratio of different sized or aged individuals in a population?

- How does a population persist for an indefinite period of time despite the fact that every individual in the population will die relatively soon (that is, what keeps births and deaths in approximate balance)?

- What is the probability that a given population will go extinct by a given future time, AND HOW IS THAT PROBABILITY INFLUENCED BY HUMAN ACTIVITIES?

Why Population Ecol MUST be quantitative:

ALL of the preceding questions involve numbers. To answer them in anything but the most cursory way requires quantitative tools. In applications, it is often not enough to know what factors influence, e.g., extinction risk. We need to have RELATIVE answers: - which factors have the most influence on populations

- which of several populations is most likely to avoid extinction

- which management strategy will reduce extinction risk more?

To provide such answers requires…

Quantitative tools:

1) mathematics

deterministic

stochastic

2) computing languages – MATLAB, C, R

We will learn about both in this course

>>Equivalent of Newton’s Law of Constant Motion:

constant PER-CAPITA growth in absence of “other forces”;

constant per-capita growth means GEOMETRIC increase (or decline)

Simplest population model:

N(t) = no. of individuals in population in year t (literally, at census point t)

Assume births immediately follow the census, then individuals may (or may not) survive to the next census.

B = average no. of surviving newborns PER CAPITA before the next census

D = fraction of adults dying before the next census (also PER CAPITA)

N(t+1) = B N(t) + (1 – D) N(t)

N(t+1) = (B – D + 1) N(t)

N(t+1) = lambda N(t)

where lambda = (B – D + 1) measures excess of births over deaths

if B = D, lambda = 1, no change in population

if B < D, lambda < 1, population declines

if B > D, lambda > 1, population grows

Prediction: total population size grows or declines GEOMETRICALLY but at a constant PER-CAPITA rate, lambda

Simulate growth of populations with constant lambda

FIRST, BASIC INTRO TO R

Opening R

Changing directory to a folder for the current project (I always do this first thing)

In Console window:

lam=1.1

lam

N=10

N=lam*N

N

n (case sensitive)

Iteration:

N=lam*N; N

alternate up-arrow and Enter

To iterate, better to write a “script” or program:

opening a script

saving a script

running a script

Script will need to:

- repeat N=lam*N tmax times

- use 3 different lams: 1

- plot results vs. time

To do this, we will need to:

- define constants (to make the script generic)

- use matrices: one name for many items N=matrix(c(4,3,6,2,1),1,5)

accessing items: N[3] (=6) using square brackets

- allocate memory to store numbers at each time – e.g., N1=matrix(0,tmax+1,1)

- repeat an action using a “for loop” – for(t in 1:tmax)

what does 1:tmax do?

- use a plotting function – help(matplot) or ?matplot

THEN WRITE FIRST PROGRAM IN R – lambda.const.R

5 (or 6) basic components of every program (or script):

[ clear memory - rm(list=ls(all=TRUE)) ]

define constants

allocate memory to variables ( to speed program )

assign initial values to variables

perform actions (e.g., loops)

output results

# define constants

lam1=0.9

lam2=1

lam3=1.1

tmax=20

n0=10

# allocate memory

N1=N2=N3=matrix(0,tmax+1,1) # could also use N1=N2=N3=numeric(tmax+1)

# initialize variables

N1[1]=N2[1]=N3[1]=n0

# perform iterations

for (t in 1:tmax){

N1[t+1]=lam1*N1[t]

N2[t+1]=lam2*N2[t]

N3[t+1]=lam3*N3[t]

}

# first run above and look at N1, N2, and N3;

# then add the following to plot results and run again

matplot(0:tmax,cbind(N1,N2,N3),xlab="Year,t",ylab="N(t)",main="Lambda Constant",type="b",pch=19)

NEXT:

First deviation from the first law: population growth rate is not constant

Now

N(t+1) = lambda(t) N(t)

lambda(t)>1 in years when B>D

lambda(t)> generate npops uniform nos. between lmin and lmax

lambar=1

dl=.1

x=runif(5000,lambar-dl,lambar+dl)

hist(x, xlim=c(.5,1.5))

similar syntax for other prob. distributions; e.g. Normal:

x=rnorm(5000,.1,.1)

hist(x, xlim=c(-1,1))

2) Multiplying a row of a matrix by a row vector:

N=matrix(10,10,5)

L=runif(5,.9,1.1)

L

N[2,]=L*N[1,]

N

L=runif(5,.9,1.1)

N[3,]=L*N[2,]

N

Now work on writing your own GENERIC script:

SCRIPT (this is lambda.variable.R):

tmax=50

npops=20

lambar=1.01

dl=.05

n0=10

lmin=lambar-dl

lmax=lambar+dl

n=matrix(0,tmax+1,npops) # 2D array n will hold all npops trajectories

n[1,]=n0

for(t in 1:tmax) n[t+1,]=n[t,]*runif(npops,lmin,lmax) # EXPLAIN THIS LINE

matplot(0:tmax, n, type="l", xlab="Year", ylab="Population size", main="Lambda variable",ylim=c(0,1.05*max(n)))

Start with npop=1 and run a few times.

Then increase the number of populations to 50, then slowly increase dl, and see what happens

It would be convenient to look at n(t) on the log scale

Here, use stoc.lambda.R instead of modifying the program above

Clearly, lambar doesn’t describe population growth, because most populations can decline even if lambar>1

What is a better measure of the population growth rate?

Geometric mean as a measure of long-term growth rate

N(1) = L(0) N(0)

N(2) = L(1) N(1) = L(1) L(0) N(0)

N(t) = L(t-1) L(t-2) … L(1) L(0) N(0)

What L (call it Lg), when multiplied by itself t times, would equal L(t-1) L(t-2) … L(1) L(0)? This is a good measure of AVERAGE ANNUAL GROWTH

Solve Lg^t = [ L(t-1) L(t-2) … L(1) L(0) ] for Lg

Answer: Lg = [ L(t-1) L(t-2) … L(1) L(0) ]^(1/t)

But, this is NOT the arithmetic mean of the L’s, which is:

(1/t) [ L(t-1) + L(t-2) + …+ L(1) + L(0) ] = lambar (or La)

Instead, it is the GEOMETRIC MEAN. That’s why we called it Lg.

The geometric mean is more strongly depressed by an L x units below the arithmetic mean than it is increased by an L x units above the arith. mean (unlike the arith. mean itself)

Why? Math fact: rel. small values decrease a product more than rel. large values increase it.

For example, what if a single L(t) = 0? What is Lg? 0.

Another example showing deviations above and below La are not equal in terms of their effects on Lg:

La = 1.1; dL = 0.2

Case 1:

even years lambda = La

odd years lambda is La plus dL

Lg=sqrt(1.1 x 1.3 ) = 1.196, which is .096 above La

Case 2:

even years lambda = La

odd years lambda is La minus dL

Lg=sqrt(1.1 x 0.9 ) = 0.995, which is .105 below La

A third way to look at Lg:

A useful approximation for the geometric mean in terms of the arithmetic mean and the year-to-year variance of lambda, Var(lam) ( Var is average squared deviation from La )

Lg ≈ La exp{ - Var(lam)/ [ 2 La^2 ] }

If Var(lam) = 0, exp{ - Var(lam)/ [ 2 La^2 ] } = 1 and Lg = La. But whenever

Var(lam) > 0, exp{ - Var(lam)/ [ 2 La^2 ] } < 1 and Lg < La

On board: plot Lg vs Var(lam) according to approximation: exponential decline with increasing Var(lam)

Show that Lg does better at “splitting the difference” between possible trajectories than does La.

PROGRAM: median.vs.lamg.&.lama.R

This program uses another way to compute Lg:

if Lg = [ L(t-1) L(t-2) … L(1) L(0) ]^(1/t)

then

log Lg = log{ [ L(t-1) L(t-2) … L(1) L(0) ]^(1/t) }

log Lg = (1/t) * [ log{L(t-1)} + log{L(t-2)} + … log{L(1)} + log{L(0)} ]

= arithmetic mean of log(L)’s

So, Lg = exp{ arithmetic mean of log(L) }

Important messages:

1) Environmental stochasticity, by introducing variation in the lambdas, actually depresses long-term population growth relative to the prediction of lambar. We can even have a population that has an (arithmetic) average growth >1 but that is virtually certain to decline over the long run if variation is sufficiently high.

2) So when we relax the assumption of the first law that lambda is constant, not only do we get variation in population growth from year to year, we get a lower long-term rate of growth.

Even if Lg>1, some possible trajectories may reach low population sizes. For a sexually reproducing species, we might consider 1 to be effective extinction (“quasi-extinction”). Some trajectories might hit 1 even if Lg>1.

Next Big Topic

What is the probability that populations will go (quasi-)extinct by some future time? Effects of mean and variance in annual growth rate

{{ SKIP THE FOLLOWING IF PRESSED FOR TIME:

First:

uniform variation not realistic, and can’t increase variation indefinitely because lambda can’t be negative, but it often has a weaker constraint above than below

Alternative: normal variation, but then lambda could still be negative

Better alternative: let lambda follow a lognormal distribution

if X ~ normal(mean=mu,sd=sig)

then Y=exp(X) ~ lognormal

lognormal.demo.R

}}}

POSSIBLY START with YGB as motivation for wanting to quantify Prob(quasi-ext.)

Note that in running stoc.lambda.R, population size was LOGNORMALLY DISTRIBUTED

The Lognormal distribution:

mu=.1

sig=.5

x=rnorm(10000,mu,sig);

split.screen(c(2,1))

screen(1)

hist(x,breaks=50)

screen(2)

hist(exp(x),breaks=50)

This means that the LOG of N will follow a CHANGING Normal distribution

But this resembles the physical process of DIFFUSION (e.g. of molecules of a gas) in a moving fluid

[ here and everywhere LOG means NATURAL LOG]

Fig. 3.3 from Morris&Doak (M&D)

Start from LOG of current population size

The lower boundary: “quasi-extinction threshold”

WE can set this boundary wherever we wish, based on BIOLOGICAL and POLITICAL considerations (some disc. of this in Ch. 2)

What determines the likelihood of hitting the threshold is:

1. How fast (and in what direction) the position of the mean of the distribution moves

2. How fast the VARIANCE of the distribution increases over time

3. How far down the population must go from current (log) size to (log) threshold

4. How long into the future we are considering

VARIANCE – mathematical definition:

average SQUARED distance between the “particles” (or pop sizes) and their (arithmetic) mean

Eg for a set of n trajectories at time t

Var(N(t)) = (1/n) * sum(over i) of [ ( Ni(t) – Nbar(t) )^2 ]

If mean declines, CERTAIN to hit threshold (and more quickly the quicker the mean declines)

If mean increases quickly, LESS likely to hit threshold quickly

If var increases quickly, MORE likely to hit threshold quickly

Mean(t) = mu * t

Var(t) = sigma_squared * t (sigma_squared = sig2 below)

mu and sig2 are constants that determine RATE of increase in mean and variance, resp., over time

MU and SIGMA – Greek letters

[ Plot Mean and Var vs t for diff. values of mu (-.1 0 .1) and sig2 (positive only) ]

The effects of mu and sig2 on likelihood of hitting threshold most easily seen by plotting the QUASI-EXTINCTION TIME CUMULATIVE DISTRIBUTION FUNCTION or CDF vs time in the future (which we will abbreviate as G(T) )

Fig. 3.5 in M&D

G(T) is PROBABILITY that quasi-ext. threshold has been reached AT ANY TIME BETWEEN NOW AND FUTURE TIME T.

Because it is a probability, it must always lie between 0 and 1.

If mu > 0, G_max < 1 (because some pops go to infinity – NOT REALISTIC)

If mu < 0, G_max = 1

The equation that governs G(T) given Nc (current pop size), Nq, mu, and sig2 has been derived by physicists.

It is eq. 3.5 in M&D

where d = log(Nc) – log(Nq)

NOTE in fig. 3.5 that for any given mu, the prob. of extinct. increases more rapidly early on if sig2 is larger.

All of this is called the DIFFUSION APPROXIMATION for estimating extinction risk.

Advantages: easy to apply

Disadvantages: several assumptions that are often unrealistic (LATER)

NEXT:

Using the DA with real data

Intro:

The YGB

- isolated population – part of once-much-larger range

- long (44 yr.) record of counts of adult females (actually 3 yr. running sum of females with YOY cubs – 3 yrs. between births, so sum is approx. of total no. of adult females)

- 6 more years of data than in M&D

- counts only – no info. on pop structure in these data

- counts at dumps before, and by aerial survey after 1973

- fires in 1988 – did it change mu or sig2?

- legally protected from hunting, but contact with humans still a source of mortality; pressure to remove protection to have a SUCCESS STORY for the law (Endangered Species Act – delisted in 2007 and then Relisted in 2010 due to political pressure from conservation groups)

QUESTION: What is the current risk of extinction? Impt. to know before delisting.

Goal: compute extCDF for YGB, but to do so we must first estimate mu and sig2

TWO methods to do so.

But first, getting data into R (follows ygb.R, but do step by step in command window

Read csv data file into a DATA FRAME

data=read.csv("ygb_females_1959_2003.csv")

Cols. of data have names, Year and N, inherited from csv file

attach datafile – now can use col’s as variables

Use plot to plot the data

plot(Year,N,etc.)

1st, what is estimate of PER-CAPITA pop growth rate over 1 yr?

lam_t = N_t+1/N_t

can get all annual lam’s in 1 step:

n=length(N)

lam=N[-1]/N[-n]

Discuss above

NOTE: one fewer lam than censuses

Method 1:

mu – how much mean LOG pop size changes per year

therefore, (arith.) average of LOG( N(t+1) / N(t) ) is an estimate of mu

sig2 – how much VAR(log pop size) changes per year

therefore, variance of log( N(t+1) / N(t) ) is an estimate of sig2

compute mu and sig2 by “standard method”

mu=mean(log(lam)); mu

sig2=var(log(lam)); sig2

Now go to ygb.R

Using mu and sig2 to compute extCDF

writing your own functions in R – eg extcdf - stored in pva.functions.R

NOTES:

input variables must be given in correct order, or variable names must be given in function call

[ functions can define other internal functions, but they will not work outside the function ]

storing functions in other files and “source”-ing them

Result: very low prob.(extinction) at t=50

BUT

using only the best estimate of mu and sig2 ignores the fact that these are only estimates, and may miss the “true” values.

We can use the confidence intervals for mu and sig2 to put confidence limits on the CDF.

Procedure:

draw a random mu from a Normal distn. with mean mu_best and var. SEmubest (=sig2/q^.5); keep if within mu’s CI, repeat if not

draw a random number from a chi-squared distn. with 1 and q df and multiply it by sig2best/q; keep if within sig2’s CI, repeat if not

compute CDF and save any values (over t) that exceed past values, both above and below CDFbest

repeat many times, and plot limits above and below

use extprob.R

NOTE: wide CI around best CDF

MESSAGE: ext. risk of YGB could be as high as 16% at T=50

BUT what if there were a gap in the census – mean and var should change more over the gap than in the one-yr. intervals between other censuses – standard method inappropriate

DIFFUSION APPROXIMATION FOR PROBABILITY OF QUASI-EXTINCTION:

When censuses are taken every year, the easiest way to estimate the parameters mu and sigma^2 for the diffusion approximation are:

mu = mean(log(N[-1]/N[-n])) = mean(diff(log(N)))

sig2 = var(log(N[-1]/N[-n])) = var(diff(log(N)))

BUT... not appropriate if some intercensus intervals are longer than others (b/c pop should change more in such intervals)

Alternative: linear regression approach –

allows for diff. intercensus intervals

other advantages:

- easy confidence interval (at least for mu)

- can use regression tools to identify outliers

- can test for changes in mu and sigma^2 in different time periods (before/after dumps closed; before/after fires)

Estimating mu and sig2 by regression:

Regress x = sqrt(diff(Year)) vs y = diff(log(N))/x WITH ZERO INTERCEPT

Slope is estimate of mu

Mean of squared residuals around regression line is estimate of sig2

Plot what this looks like.

Doing a linear regression in R:

x=1:10

y=3+2*x+rnorm(10,0,1)

plot(x,y)

out=lm(y~x)

summary(out)

coef(out)

confint(out)

anova(out)

sig2=anova(out)[2,3]; sig2

s2=sum((out$resid)^2)/8; s2

CIsig2=(q-1)*sig2/qchisq( c(.975,.025),df=(q-1) )

IN-CLASS ASSIGNMENT, WORKING IN PAIRS:

1988 WAS THE YEAR OF THE YELLOWSTONE FIRES. IMAGINE THAT, BECAUSE OF THE FIRES, IT WAS NOT POSSIBLE TO DO THE GRIZZLY BEAR CENSUS THAT YEAR. USING THE YGB DATA, DELETE THE COUNT FROM 1988, ESTIMATE MU AND SIGMA^2 BY LINEAR REGRESSION, AND USE THOSE ESTIMATES TO PRODUCE A QUASI-EXTINCTION TIME CDF USING THE PROGRAM YGB.R.

other advantages of regression approach:

tests for outliers

tests for changes in mu and sig2 (e.g after 1983 fires)

confidence intervals on mu (can also be calculated directly from mu and sig2)

doing regression with YGB data: see more.ygb.stuff.R

Before running any of this prog – save standard estimates for comparison:

mu_s=mu

sig2_s=sig2

confidence limits on mu produced directly in bint

confidence limits on sig2 computed using chi2 distn.

[possibly skip or give overview]:

Two ways to look for outliers using regression output:

dffts

Rstudentized

Both indicate 1983 is odd – UNUSUALLY HIGH – consequences for estimated extinction risk?

1983 is not the year the dumps were closed, or the fire year. If a reason to discard it is determined, could delete this lambda and estimate mu and sig2 using the remaining lambdas.

Can also ask statistically if mu, sig2 change before/after fires or before/after dumps closed – see details in more.ygb.stuff. R

Uses of CDF

comparison of different pops w/ diff. Nc, mu, or sig2, diff thresholds [ M&D figs. 3.9,3.10 ]

Review and tests of assumptions:

I. Parameters mu and sig2 constant

Violations:

1) dens density dependence could change mu (and even sig2)

mu declines as N increases

mu could decline as N decreases

2) dem. stoch could change sig2

3) environmental trends could change both mu and sig2

II. No envtal autocorrelation

Whether lam was large last year has no effect on whether lam is lg. or small. this year

We’ll see how to test for and incorporate this next time

III. No extremely large or small values of lam (no bonanzas or catastrophes)

tests for outliers

how to include if found

IV. No observation error

- counts are assumed to be accurate – if not, they inflate sig2, making calculated ext. risk too high

HOMEWORK:

Read Ch. 4 (skip Ceiling Model) including Appendix

Next Big Topic: More complex count-based models:

- Density-Dependence (negative and positive, i.e., Allee effects)

- Environmental Autocorrelation and its interaction with D.D.

- Catastrophes and Bonanzas

(for later: Demographic Stochasticity)

I. An overview of models with negative density dependence

[ SKIP 1. Ceiling model (eq. 4.1 in M&D) ]

2. More realistic models with continuous change in lambda as N(t) increases

The DI model assumed loglam indep. of N(t), with mean mu that doesn’t change as N(t) changes. Var in loglam around mu is sig2, assumed to be caused by environmental variation.

[[

Why log?

1) log ( N(t+1)/N(t) ) = log ( N(t+1) ) – log( N(t) )

errors in estimating N(t+1) and N(t) have same effects on log lambda,

but errors in estimating denominator of N(t+1)/N(t) have much stronger effect than errors in estimating numerator

2) log ( N(t+1)/N(t) ) can go from –Inf to Inf, but N(t+1)/N(t) only goes from 0 to Inf (if N(t)>0). Log ratio more likely to be normal, and so tools that assume normal variation more appropriate.

]]

But more realistically, we might expect the mean log lambda (i.e., mu) to decline as N(t) increases. This is called “NEGATIVE density dependence”, because log lambda DECLINES.

Several patterns for this decline. Illustrated by the so-called “Theta Logistic Model”:

(on board)

N(t+1) = N(t) * exp{ r*( 1 – [N(t)/K]^theta ) }

Take logs of both sides:

log( N(t+1)/N(t) ) = r*( 1 – [N(t)/K]^theta )

Plot log population growth rate log ( N(t+1)/N(t) ) vs N(t)

show.theta.logistic.R

Patterns:

log lambda decreases linearly with N(t) (theta=1); this model is also known as the RICKER MODEL, widely used in fisheries modeling

log lambda decreases sharply only as N(t) approaches K (like ceiling model); theta>1

log lambda decreases sharply at first, then changes little as N(t)->K; theta0).

FIRST,

> how do we decide which of these patterns of density dependence best describes a given dataset?

> what are the consequences of negative density dependence for extinction risk?

To illustrate this, we will use a new data set. The Bay Checkerspot Butterfly, Euphydryas editha bayensis, the subject of a long-term population study by Paul Ehrlich’s laboratory at Stanford University in CA, USA.

Run checkerspot.m to observe counts (arith. and log. scales) and loglam vs Nt and logNt

II. Testing for density dependence using maximum likelihood and AIC

[ Probably skip to *****************, except these parts

Review of “maximum likelihood” parameter estimation with Normal errors (from Ch. 4 appendix of M&D)

Example: fitting the theta logistic function:

log( N(t+1)/N(t) ) = r*( 1 – [N(t)/K]^theta )

^^ Like a regression equation:

dep. variable: y[t] = log( N(t+1)/N(t) )

indep. var.: N(t)

params. or coefficients: r, K, theta

Rewrite eq. as

y[t] = f(p,N[t])

where y[t] is log growth rate in year t, the DEPENDENT VARIABLE, p=[r, K, theta] is a vector of parameter values, and f(p,N[t]) = p[1]*( 1 – ( N(t)/p[2] )^p[3] ) is the theta logistic function, where N[t] is the INDEPENDENT VARIABLE

What we are trying to do in maximum likelihood parameter estimation is to find the values of the parameters (or the “value” of the parameter vector p) that maximize the probability of observing the y[t]’s given the N[t]’s (and, of course, given the particular form of the model, f(p,N[t]) )

The Normal probability of seeing log growth rate y[t] at time t given a value of p and N[t] is

Pr{y[t] | p,N[t]}=[pic] (use show.normal.R to plot this)

where

[pic] is the average squared deviation between the observations and the prediction of the model (the “residual variance”).

We want to pick p so that f(p,Nt) is close to yt, because then Pr(y|p,N) will be maximum. But we want f(p,N) to be close to ALL yt’s given the Nt’s, so we may have to compromise.

The OVERALL probability of seeing ALL the data is the product of these probabilities over all times, keeping the model and the parameter values fixed as we cycle through the pairs of dep. and indep. variables. This overall probability is the LIKELIHOOD L of the observed data given the model equation and the value of p:

L = [pic]

where q is the number of pairs of values of the dep. and indep. variables (the “SAMPLE SIZE”)

Because this is a product of many small numbers (because prob’s are between 0 and 1), it will be very small. To prevent rounding errors, we take its log to get the LOG LIKELIHOOD

log L = [pic]

=[pic]

=[pic]

=[pic]

= [pic]

BUT…from the definition of [pic]:

[pic]

************************************

Therefore

[pic]

“THE LOG LIKELIHOOD FUNCTION”

So, all we need to compute the log likelihood function with Normal errors is the sample size q and the residual variance [pic], but as seen in the eq. for [pic], we need p and N[t] to compute it.

So what we do is to use a search algorithm (a minimization routine) to find the value of p, given the sequence of N’s and therefore the y’s (the log lambdas computed from the N’s) that MAXIMIZES THE LOG LIKELIHOOD. Because most routines are set up to find the minimum of a function, we will search for the value of p that MINIMIZES THE NEGATIVE LOG LIKELIHOOD. Because the NLL for normally distributed errors contains only Vr (and q, which is fixed), minimizing NLL is equivalent to minimizing Vr, or minimizing the sum of squared deviations between the data points and the model. For this reason, maximum likelihood fitting with Normal errors is equivalent to LEAST SQUARES PARAMETER ESTIMATION.

Now, to ask if there is negative DD in the data, and if so what form it takes, we need to fit several models AND TO CHOSE WHICH MODEL IS BEST. Using maximum likelihood estimation and AIC, there is a natural way to do this.

We’ll fit 3 models (which happen to be “nested” – simpler ones can be obtained by constraining parameters in more complex ones to particular values):

f(p,N[t]) No. parameters (incl. Vr):

The DI model: r 2

The Ricker model: r(1-N/K) 3

The theta-log model: r(1-(N/K)^theta 4

Easiest fitting methods differ (although we could use the most complex method for all models):

For DI:

the best r is simply the mean log lambda. Vr is the (biased) variance of the log lambdas (around r=mu)

Because we only need Vr and q to compute logL, we can do so directly once we have the variance of the log lambdas

using Nt from checkerspot.m:

loglamt=diff(log(Count));

r=mean(loglamt)

Vr=mean( (loglamt - r)^2 )

% NOTE: can also use Vr=(q-1)*var(loglamt)/q

For Ricker:

loglam is a linear function of N, so we can estimate r and K with a linear regression, where Vr is the residual variance (the mean squared deviation between the points and the linear regression line).

[ students can write code to do this following the syntax for “lm” in ygb.R ]

Finally, the theta-logistic model is nonlinear, so we must use a nonlinear fitting procedure to estimate its parameters. We use nlm (non-linear minimization):

Note: we could also use the R function “optim” to minimize the neg. log likelihood function directly

Finally, having obtained Vr for all 3 models, we can compute AICc (the corrected Akaike Information Criterion, where “correction” is for small sample size) for all three models.

The best model has the smallest AICc

AICc = −2*logL + 2* (p * q ) / (q – p – 1 )

logL is negative. The better the fit of a model, the larger the log likelihood (i.e., the less negative it is, so the smaller is the first term in AICc. In general, more parameters (higher p) should increase logL, and so decrease this first term. However, as p increases, the second term increases.

[ as q gets large relative to (p+1), the 2nd term becomes 2*p*q/q = 2*p , and AICc converges to AIC = −2*logL + 2*p]

Thus AICc (and AIC) attempts to achieve a balance between goodness-of-fit (measured by logL) and the number of parameters used to achieve that fit.

For JRC population, the AICc values are:

Model p logL AICc

1 DI 2 -41.26566 87.05306

2 Rick 3 -37.79878 82.68847

3 TLog 4 -37.10478 84.11433

SO, even though the TLog model has the highest (least negative) logL, it uses 1 more parameter to achieve this only-slightly-better fit than does the Ricker, so it has a higher AICc. In constrast, the Ricker has much lower logL than the DI, and has the lowest AICc of all models. Therefore we conclude that the model with linear negative DD (the Ricker) is the best model.

NEXT…

Simulating ext. prob. using the Ricker

We’ll need best estimate of r, K, Vr, and the last population size and the QET.

Using the script theta.logistic.R

> uses rnorm

> if N falls below nx, it is set to 0, because it will then always remain below nx

Results:

* tight overlap of lines means 50K trajectories is sufficient to characterize ext. risk FOR A GIVEN SET OF PARAMETERS – we have not incorporated parameter uncertainty, but could do so as we did using extprob.R (but we would then have to grapple with the fact that estimates of r and K are not independent)

* high risk of extinction for the best-fit parameters

INDEED, this population went extinct 10 yrs. after the last census. The high value of Vr may be why

NOTE: unlike in DI model, where Pr(ultimate extinction)0, with a DI model, it is always 1.

Consequences of negative DD in discrete time

Ricker fitted to checkerspot data shows high extinction risk.

But this is probably because the estimate of Vr is high, not because the pop. is predicted to oscillate (due to the 1 big peak in plot of fitted Ricker) – we might want to treat that as an outlier.

But with neg. density dependence, there can be another distinct extinction risk: population oscillations causing pop size to occasionally visit low numbers

Period-doubling bifurcations in the Ricker model

observe how trajectory and recruitment curve changes as r changes, using ricker.R

r = 1.9, 2.4, 2.6, 2.7, 3

Ricker model has “overcompensatory” (neg) DD –

- when above K, pop can decline below it in 1 time step

- when below K, pop can climb above it in 1 time step

depending on steepness of recruitment curve

Then use bifurc8.R to produce a bifurcation diagram

1. Positive density dependence: Allee effects

Models above assume that loglam only declines with increasing Nt. But it could also deline at DECREASING Nt, a phenomenon known as an

ALLEE EFFECT (after Warder Clyde Allee)

Potential causes:

- declining birth due to difficulties finding mates at low densities

- declining survival due to failure of group defense (including against abiotic conditions, as in conspecific nurse effects) or group foraging

A model with declining birth at low density (e.g. due to reduced mating) and declining survival at high density (due to resource limitation)

Birth = B(N) = a + (r−a) N/(A+N) [ Note: per-capita ]

if N=0, B=a

if N=Inf, B=a+(r−a) = r

when N=A, B=a + .5r - .5a = (a+r)/2 (birth is half way between its minimum and maximum values when N=A; A is the “half-saturation constant”)

Survival = S(N) =exp(-b*N)

if N=0, S=1

if N=Inf, S=0

Assume univoltine or annual life cycle:

Lambda = B * S = [a + (r-a) N/(A+N) ] exp(-bN)

(this is equivalent to eq. 4.12 in M&D if a=0, r=exp(r), and b=beta)

Plot B, S, Lambda, and N(t+1) vs N(t) PROGRAM allee.R

a = L(N=0)

if a>1, L(0) > 1 “WEAK” Allee effect

if a 2 SDs above/below mean

- Using extremes.R to incorporate outliers

replace “run-of-mill” values with extremes, with freq. determined from data

Issues:

* match single pos. (or neg.) outliers with complements?

* may need to modify extremes.R if there are BOTH outliers and density dependence (and other factors), as in checkerspot

Next Topic: Structured populations

Rationale:

We will begin by going backwards and removing several things we added to count-based models (environmental stochasticity and density dependence). After we see how the dynamics of structured populations behaves without these complications, we will put them back in. In addition, we will consider the effects of another force, demographic stochasticity.

Basic reason why structure matters:

Individuals don’t contribute equally to population growth

Example: semi-palmated sandpiper - PHOTO

REFERENCE: Hitchcock and Gratto-Trevor (Ecology 78:522-533, 1997)

Biological background:

individuals can live >3 yrs, can begin breeding at age 1, migrate to N. Canada, produce 1 nest per year. Common, but H> studied a declining population to identify causes of decline

Data on average vital rates

s0=.1293 # juvenile survival – prob(survival) from hatching to age 1

s1=.2543 # yr1 survival – prob(survival) from age 1 to age 2

s2=.593 # adult survival – 1yr. prob(surv) for all individuals age 2 or more

b1=.25 # prob(breed) at age 1

b2=.875 # prob(breed) at age 2

b3=.95 # prob(breed) at ages 3 and above

c=1.8625 # chicks per nest

Clearly, age strongly affects survival and likelihood of breeding.

Therefore separating individuals into age classes and keeping track of their separate contributions to next year’s population should improve our ability to predict/understand population growth rate ( a pop. of all juveniles will grow differently than a pop. of all adults)

Building a density-indep. model. As for the unstructured case, we’ll follow population growth over discrete 1 yr. intervals.

First, consider when we are censusing the population. That will determine what age (or size) classes we see.

Hitchcock and Gratto-Trevor censused just BEFORE nests produced, after birds returned to N. Canada on Spring migration – A PRE-BREEDING CENSUS

So they see 1YOs (born just after last census), 2YOs, and 3+YOs. Just after the census, new juveniles are born. The picture is this:

[pic]

Arrows are PER-CAPITA rates.

We could write equations for the change in ea. age class over 1 yr:

Let N1(t), N2(t), N3(t) == no. of 1, 2 and 3+ YOs, resp., in the population in yr. t

Sum up products of arrows times Ni(t):

N1(t+1) = b1*c*s0*N1(t) + b2*c*s0*N2(t) + b3*c*s0*N3(t)

N2(t+1) = s1*N1(t)

N3(t+1) = s2*N2(t) + s2*N3(t)

Using linear algebra, we can write this more concisely by separating the N terms from the rest of the r.h.s.:

n(t+1) = A * n(t)

where n(t) = [ N1(t) N2(t) N3(t)]’ and

A= [ b1 c s0 b2 c s0 b3 c s0; s1 0 0; 0 s2 s2]

Note: A contains all the PER-CAPITA effects, and n contains the numbers in ea. class

Using A as a table

cols = class this year

rows = class next year

elements: PER-CAPITA contributions from this year’s to next year’s population

Reviewing right mult of a matrix by a column vector - numerical example with starting vector [1 1 1]’

If break between classes, use life-cycle diagram to (re)illustrate the matrix/life history events – timing of annual census determines life stages “seen”:

[pic]

6 arrows above correspond to the 6 terms in the matrix, and represent one-year transitions

Put numbers in matrix (“POPULATION PROJECTION MATRIX”), and interpret matrix

cols = from life stage in yr. t

rows = to life stage in yr. t+1

entries are PER-CAPITA RATES

Rules for right-multiplying a matrix by a column vector – if dim(vec,2)=dim(mat,1)=dim(mat,2)

RESULT: another col. vec. of same dimension as original vec

If vector is [n1(t); n2(t); n3(t)], show that this gives us back the separate recursion equations for ea. life stage

Go to R to write a program to predict what the population will do when we repeatedly multiply the matrix A by the (changing) vector n

First, use sandpiper.matrix.R to show 3 ways to construct a proj. matrix in R

Do some matrix/vector multiplication in the command window

n=c(1,1,1) OR n = matrix(1,3,1)

n=A %*% n; n ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download