“First Law” of Population Dynamics (like Newton’s …
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.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related searches
- examples of the first law of motion
- newton s first law illustration
- newton s first law of motion formula
- newton s first law formula
- first law of planetary motion
- johannes kepler first law of planetary motion
- first law of planetary motion drawing kepler
- newton s first law worksheet answers
- newton s first law equation
- newton s first law definition and explanations
- newton s first law of motion
- newton s first law of motion sentence