In brouwern/popcycles: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

Code to implement the matrix population model of Wootton and Bell (1992). This code attempts to reproduce their original results and figures using the methods described in the paper. Hunter and Caswell (2005) describe an alternative way to approach the analysis but only carry out a portion of the analysis.

A note on terminology: metapopulation

The original paper refers to the model as a metapopulation model, as do Hunter and Caswell (2005). However, its generally considered best to reserve "metapopulation" for systems where extinction of subpopulations and recolonization vacated patches is common. A more appropriate term might be "spatially structured matrix population model."

Wooten and Bell 1992. A METAPOPULATION MODEL OF THE PEREGRINE FALCON IN CALIFORNIA: VIABILITY AND MANAGEMENT STRATEGIES. Ecological Applications 2: 307-321.

Key ecological details

TODO: add important details

Preliminaries

Packages

General functions for demography are contained in the popcycles package.

library(popcycles)

The magic package makes working with matrices easier.
TODO: check if I actually use this or if its just for Hunter and Caswell

#install.packages("magic")
# library(magic)

The Wootton and Bell model

Wootton and Bell (1994) provide the paramters for their model in Table 1 (pg 312).

TODO: build the table and render in markdown e.g. using pander()

Parameters

Source of parameters

TODO: brief description of how parameters obtained.

"Basic (unstructured) model"

## 
# demographic parameters
b <- 0.69 # no. female fledgling per paired female
y <- 0.36 # probability fledgling surviving to age 1
r <- 0.72 # probablity age 1 bird surviving to year 2 and acquring territory
s <- 0.77 # survival of territorial adult

# Note: 
r*y # = 0.2592

# Introduction parameters
## Note: Authors use I, which is an R function
## I add 0 or 50 to the object name.
I0  <- 0
I50 <- 50

Initial population sizes

## territorial females
N.n <- 60
N.s <- 29

## non-territorial females
F.n <- 19
F.s <- 20

Metapopulation model

## 
# demographic parameters
b.n <- 0.71 # no. female fledgling per paired female
b.s <- 0.53 # no. female fledgling per paired female

y.n <- 0.36 # probability fledgling surviving to age 1
y.s <- 0.36 # probability fledgling surviving to age 1

r.n <- 0.72 # probablity age 1 bird surviving to year 2 and acquring territory
r.s <- 0.72 # probablity age 1 bird surviving to year 2 and acquring territory

s.n <- 0.77 # survival of territorial adult
s.s <- 0.77 # survival of territorial adult

m   <- 0.27 # probability of a nonterritorial female migrating into the other subpopulation

#fraction of captive-reared fledglings introduced into the northern subpopulation
f0.0 <- 0.0
f0.5 <- 0.5
f1.0 <- 1.0

Mean Matrix

A mean matrix averaging the Northern and Southern subpopulations.

Equation 1, pg 309, Wootton and Bell 1992

The general form of the matrix is

TODO: put this in latex

0 y*b r s

# mean matrix
Mbar <- matrix(data = c(0, y*b,
                        r, s),
            byrow = T, 
            nrow = 2)

# round y*b to look for rounding errors
Mbar.round <- matrix(data = c(0, round(y*b,2),
                              r, s),
            byrow = T, 
            nrow = 2)

# vector to represent introductions
Ix50 <- c(y*I50, 0)

Meta-population matrix

Equation 3 pg 310. They use h, which is 1 -m

c = probability of surviving migration; assumed to be 1

A compact description of the model is available in a follow-up paper by the author's: file:///Users/nlb24/Downloads/appendixA.htm

NOTE: in some versions of the model there is no dispersal between the northern and southern subpopulations. This results in a reducible model, which is not good mathematically, though ecologically realistic. I believe lambda for this matrix is equal to the average of lambda of the two submatrices, but I need to check. Reducible models may also be non-ergodic. Function is_irreducible() checks for irreproducibility, while is_ergodic checks for non-ergodicity. Both functions were written by Iain Stott.

Stott et al. (2010) explain the important of reproduciblity and ergodicity:

"To be useful for predictive or prospective analyses, PPM models should generally be irreducible (the associated life cycle graph contains the necessary transition rates to facilitate pathways from all stages to all other stages) and therefore ergodic (whatever initial stage structure is used in the population projection, it will always exhibit the same stable asymptotic growth rate)" Stott et al. https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/j.2041-210X.2010.00032.x

TODO: additional functions in Stott et al may be useful. Also add their test matrices to test/demo the functions

c<- 1
Mmeta <- matrix( data = c(0.00,       y.n*b.n, 0.00,      0.00,
                          r.n*(1-m),  s.n,     r.n*m*c,   0.00,
                          0.00,       0.00,    0.00,      y.s*b.s,
                          r.s*m*c,    0.00,    r.s*(1-m), s.s),
           nrow = 4,
           byrow = T)

Projecting the model

Initial conditions

# Initial pop vector for mean model
Nbar.init <- c(F.n+F.s, N.n+N.s)
sum(Nbar.init)

# Full Initial pop vector for meta model
Nmeta.init <- c(F.n, N.n, F.s, N.s)
sum(Nmeta.init)

Mean matrix, no re-introduction

x <- rep(NA, 125)
pop.historyMbar <- data.frame(year = 1:length(x),
                                Fx = x,
                                 Nx = x)
pop.historyMbar[1,2:3] <- Nbar.init

Nbar.i <- Nbar.init

# Project using for() loop
for(i in 2:nrow((pop.historyMbar)) ){

  # nultiple mean matrix Mbar by current popualtion vector Nbar.i
  Nbar.i <- Mbar%*%Nbar.i

  # Store current popualtion vector Nbar.i 
  pop.historyMbar[i, 2:3] <- Nbar.i
}

Metapop - no re-introduction

x <- rep(NA, 125)
pop.historyMmeta <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Nn = x,
                                 Ntot = x)
pop.historyMmeta[1,2:5] <- Nmeta.init

Nmeta.i <- Nmeta.init

# project using for loop
for(i in 2:nrow((pop.historyMmeta)) ){
  Nmeta.i <- Mmeta%*%Nmeta.i
  pop.historyMmeta[i, 2:5] <- Nmeta.i
}

pop.historyMmeta$Ntot <- apply(pop.historyMmeta[,c(3,5)],MARGIN = 1,sum)

tail(pop.historyMmeta)

Mean matrix, with introduction

TODO: could create a function that takes the arguments of a matrix, starting population vector, and reintroduction amount (0, >0) and runs the projects. Not necessary but might simplify code further down in the model

x <- rep(NA, 125)
pop.histMbar2 <- data.frame(year = 1:length(x),
                                Fx = x,
                                 Nx = x)
pop.histMbar2[1,2:3] <- Nbar.init

#introduction vector
I.vector <- c(y*I50/2,0)

Nbar.i <- Nbar.init
for(i in 2:nrow((pop.histMbar2)) ){

  # Project population
  Nbar.i <- Mbar%*%Nbar.i 

  # Add introduced birds
  Nbar.i <- Nbar.i  + I.vector

  # Store current population state
  pop.histMbar2[i, 2:3] <- Nbar.i
}

Metapop - with re-introduction

25:25 allocation of offspring to North vs. South subpopulation. f = 0.5

x <- rep(NA, 125)
pop.historyMmeta_25.25 <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Ns = x,
                                 Ntot = x)
pop.historyMmeta_25.25[1,2:5] <- Nmeta.init

#introduction vector
I.vector_25.25 <- c(f0.5*y*I50/2, 
                    0, 
                    (1-f0.5)*y*I50/2, 
                    0)


Nmeta.i <- Nmeta.init
for(i in 2:nrow((pop.historyMmeta_25.25)) ){

  Nmeta.i <- Mmeta%*%Nmeta.i + I.vector_25.25

  pop.historyMmeta_25.25[i, 2:5] <- Nmeta.i
}

pop.historyMmeta_25.25$Ntot <- apply(pop.historyMmeta_25.25[,c(3,5)],
                                     MARGIN = 1,
                                     sum)

Plot results - Figure 3

This qualitatively matches Figure 3 of Wootton and Bell 1992

Note - they plot only breeders, not unpaired individuals

Data in table 1 indicates initial breeders = 89; Figure 3 shows breeders = ~ 99 Figure 4 indicates ~98, matching figure 3 Figure 5 indicates ~89 (panel B+C), matching my calcs

Could be which year they plot as year 1! Figure 3, they appear to plot year 2 as initial

plot(Nx ~ year, pop.historyMbar, pch = 7,
     ylim = c(0, 150),
     xlim = c(1,100),
     main = "Figure 3")
points(Nx ~ year, pop.histMbar2,   pch = 9)
points(Ntot ~ year, pop.historyMmeta, pch = 15, col = 2)
points(Ntot ~year, pop.historyMmeta_25.25, pch = 18, col = 2, cex =1.5)
abline(h = Nbar.init[2])
abline(v = 2)

Metapop - re-introduction with varying allocation

25:25

f0.5 = 25:25 allocation of offspring, run above f = 0.5

50:0

x <- rep(NA, 125)
pop.historyMmeta_50.0 <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Nn = x,
                                 Ntot = x)
pop.historyMmeta_50.0[1,2:5] <- Nmeta.init

#introduction vector
I.vector_50.0 <- c(y*I50/2, 0, 0, 0)


Nmeta.i <- Nmeta.init

for(i in 2:nrow((pop.historyMmeta_50.0)) ){
  Nmeta.i <- Mmeta%*%Nmeta.i + I.vector_50.0
  pop.historyMmeta_50.0[i, 2:5] <- Nmeta.i
}

pop.historyMmeta_50.0$Ntot <- apply(pop.historyMmeta_50.0[,c(3,5)],MARGIN = 1,sum)

0:50

x <- rep(NA, 125)
pop.historyMmeta_0.50 <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Nn = x,
                                 Ntot = x)
pop.historyMmeta_0.50[1,2:5] <- Nmeta.init

#introduction vector
I.vector_0.50 <- c(0, 0, y*I50/2, 0)


Nmeta.i <- Nmeta.init
for(i in 2:nrow((pop.historyMmeta_0.50)) ){
  Nmeta.i <- Mmeta%*%Nmeta.i + I.vector_0.50
  pop.historyMmeta_0.50[i, 2:5] <- Nmeta.i
}

pop.historyMmeta_0.50$Ntot <- apply(pop.historyMmeta_0.50[,c(3,5)],MARGIN = 1,sum)

head(pop.historyMmeta_50.0)
head(pop.historyMmeta_0.50)

Plot results - Figure 4

This qualitatively matches Figure 4 of Wootton and Bell 1992, except for variation likely due to ambiguity about what initial population sizes were

They appear to plot year 5 as the starting point, perhpas to get rid of initial transient dynamics.

plot(Ntot ~ year, pop.historyMmeta_25.25,
     ylim = c(85, 105),
     xlim = c(2,100),
     pch = 7,
     main = "Figure 4")
points(Ntot ~year,
       pop.historyMmeta_50.0, pch = 18, col = 2, cex = 1.3)
points(Ntot ~year,
       pop.historyMmeta_0.50, pch = 15, col = 2)
abline(v = 2)

Results

Results for models without immigration reported in Table 2

Lambda

In table 2 they report population growth rate "R" as 0.9574.

Different between my results may be due to rounding error, analytically routine for eigenvalue calculation (they did this in 1994...) or rounding error if they calculated lambda numerically by simulating population dynamics.

eigen(Mbar)         # 0.956903
eigen(Mbar.round)   # 0.9579092

calc_lam(Mbar)

The value reported for their original single matrix representing the metapopulation matches the value reported for the vec-perm approach in Hunter and Cawell 2005.

calc_lam(Mmeta) # 0.9433222 vs 0.9438 in Table 2

Sensitivity

General form of matrix 0 y*b r s

Values my sensitivities / elasticities are within rounding error for r and s.

TODO: do analysis/set up function to calcualte lower-level sensitivites y*b

calc_S(Mbar)

Note: their elasticities sum to >1.0. For matix elements, sum(e.ij) should be 1. Is this also true for elasticities of lower-level parameters?

E.Mbar.table2 <- c(0.16,0.16,0.16,0.68)
sum(E.Mbar.table2)

calc_S(Mmeta)

Note: elasticities of lower-level parametrs >1.0. See note above.

sum(c(0.11,0.06,0.11,0.06,0.11,0.57,0.42,-0.04))

Metapop model with updated adult survival

They assessed their model against observed population counts from the previous decade and concluded their estimate of adult survival was inappropriate for the North population. They updated it to 0.91 based on various calculations I haven't yet reproduced because I don't have the original abundance data they used (could be extracted from figures)

Update matrix

s.n <- s.n91 <- 0.91
s.s <- s.s80 <- 0.80 #is this updated or left at 0.77?
Mmeta.Sn.91 <- matrix( data = c(0.00,       y.n*b.n, 0.00,      0.00,
                                r.n*(1-m),  s.n,     r.n*m*c,   0.00,
                                0.00,       0.00,    0.00,      y.s*b.s,
                                r.s*m*c,    0.00,    r.s*(1-m), s.s),
           nrow = 4,
           byrow = T)

Build their ceiling model

They implement their ceiling model by swapping in an entirely different matrix when they hit the carrying capacity for the north (I haven't identified yet if they say what they do for the south - maybe they just stop the simulation / cap things)

See equation 5 (?)

Note: This new matrix introduces floaters, which are age 1 and stay in stage class 1 (cell 1,1) r.n*(1-m) term for recruitment in north (cell 2,1) gets hard coded to 0

The main ceiling model is structure like this:

d.n <- s.n
Mmeta.ceiling <- matrix( data = c(d.n*(1-m), y.n*b.n, d.n*m*c,   0.00,
                                  0.00,      s.n,     0.00,     0.00,
                                  0.00,      0.00,    0.00,      y.s*b.s,
                                  r.s*m*c,   0.00,    r.s*(1-m), s.s),
           nrow = 4,
           byrow = T) 

New vector defined for recruitment

(1-s.n) is the number of territory individuals that die (1-s.n)*Tn is the number territories that open up due to death

As part of the recruitment vector (1-s.n)*Tn is the number of "floaters"/first years added to the territory holder population.

-1(1-s.n)Tn is the number of "floaters"/first years pulled out from that class and upgraded.

Not this doesn't include any introductions - assumes introduces cease when population reaches carrying capacity?

T.n <- 100
T.s <- 120
I.ceiling <- c(-1*(1 - s.n)*T.n, (1-s.n)*T.n, 0, 0)

Run ceiling model

As simple function to help debugging.
TODO: as noted above, we could write a general function that carries out the anlysis entirely

ceiling_function <- function(){
  x <- rep(NA, 200)
pop.historyMmeta_ceiling <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Ns = x,
                                 Ntot = x)
pop.historyMmeta_ceiling[1,2:5] <- Nmeta.init



  Nmeta.i <- Nmeta.init
for(i in 2:nrow((pop.historyMmeta_ceiling)) ){

  # if(i == 8){
  #   browser()
  # }
  #  



  N.n.current <- Nmeta.i[2] 
  if(N.n.current < T.n){
  Nmeta.i <- Mmeta.Sn.91%*%Nmeta.i 
  }

  if(N.n.current >= T.n+1){
   Nmeta.i <- Mmeta.ceiling%*%Nmeta.i  + I.ceiling


  }

    #Create hard cap at 100
  # if(Nmeta.i[2] > 100){
  #   Nmeta.i[2] <- 100
  # }


  # extinction
  if(i == 100){
    Nmeta.i[1:2] <- 0
  }

  pop.historyMmeta_ceiling[i, 2:5] <- Nmeta.i

}
  return(pop.historyMmeta_ceiling)
}

Results qualitatively similar. Their model doesn't appear to set a hard boundary at 100 breeding birds in the North.

If I put one in, it creates output with very different equilibrium density for southern pop than what hey show or what I get without the hard cap.

Their plot shows different initial population sizes, and trimming off the bit of transience at the beginning doesn't totally bring it in line.

In their model the southern population also reaches the same stable population level (~90) before and after the population crash at t = 100

# x <- rep(NA, 200)
# pop.historyMmeta_ceiling <- data.frame(year = 1:length(x),
#                                  Fn = x,
#                                  Nn = x,
#                                  Fs = x,
#                                  Ns = x,
#                                  Ntot = x)
# pop.historyMmeta_ceiling[1,2:5] <- Nmeta.init




pop.historyMmeta_ceiling <- ceiling_function()
pop.historyMmeta_ceiling$Ntot <- apply(pop.historyMmeta_ceiling[,c(3,5)],1,sum)
pop.historyMmeta_ceiling$Nnorth <- apply(pop.historyMmeta_ceiling[,c(2,3)],1,sum)
pop.historyMmeta_ceiling$Nsouth <- apply(pop.historyMmeta_ceiling[,c(4,5)],1,sum)


#pop.historyMmeta_ceiling$Ntot <- apply(pop.historyMmeta_ceiling[,c(3,5)],MARGIN = 1,sum)

head(pop.historyMmeta_ceiling,40)
#tail(pop.historyMmeta_ceiling)

plot(Nn ~ year, data = pop.historyMmeta_ceiling, pch = 7,
     xlim = c(0,200),
     ylim = c(0,150))
abline(h = 100, col = 2)
points(Ns ~ year, data = pop.historyMmeta_ceiling, pch = 18, cex = 1.2)
abline(v = 2)

Ceiling model with no migration

Mmeta.Sn.91.nomig <- matrix( 
                       data = c(0.00,       y.n*b.n, 0.00,      0.00,
                                r.n,        s.n,     0.00,      0.00,
                                0.00,       0.00,    0.00,      y.s*b.s,
                                0.00,       0.00,    r.s,       s.s),
           nrow = 4,
           byrow = T)

They also run a ceiling model version (Figure 7) where migration is not allowed.

Mmeta.ceiling.nomig <- matrix( 
                         data = c(d.n,       y.n*b.n, 0.00,     0.00,
                                  0.00,      s.n,     0.00,     0.00,
                                  0.00,      0.00,    0.00,     y.s*b.s,
                                  0.00,      0.00,    r.s,      s.s),
           nrow = 4,
           byrow = T) 

Run ceiling model

As function to help debugging

ceiling_function <- function(){
  x <- rep(NA, 100)
pop.historyMmeta_ceiling <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Ns = x,
                                 Ntot = x)
pop.historyMmeta_ceiling[1,2:5] <- Nmeta.init



  Nmeta.i <- Nmeta.init
for(i in 2:nrow((pop.historyMmeta_ceiling)) ){

  # if(i == 8){
  #   browser()
  # }
  #  



  N.n.current <- Nmeta.i[2] 
  if(N.n.current < T.n){
  Nmeta.i <- Mmeta.Sn.91.nomig%*%Nmeta.i 
  }

  if(N.n.current >= T.n+1){
   Nmeta.i <- Mmeta.ceiling.nomig%*%Nmeta.i  + I.ceiling


  }

    #Create hard cap at 100
  # if(Nmeta.i[2] > 100){
  #   Nmeta.i[2] <- 100
  # }


  pop.historyMmeta_ceiling[i, 2:5] <- Nmeta.i

}
  return(pop.historyMmeta_ceiling)
}

pop.historyMmeta_ceiling_nomig <- ceiling_function()

pop.historyMmeta_ceiling_nomig$Ntot <- apply(pop.historyMmeta_ceiling_nomig[,c(3,5)],1,sum)

Results qualitatively similar to figure 7 Do they have mig-vs no mig labels backwards?
w/o migration population should grow FASTER This is what mine shows,

They approach asymptote for southern population more smoothly

## 7 A
plot(Nn ~ year, data = pop.historyMmeta_ceiling[1:99,], pch = 7,
     xlim = c(0,100),
     ylim = c(0,150),
     ylab = "Breeding pairs",
     main = "Figure 7A",
     xaxs="i",
     yaxs="i")
abline(h = 100, col = 2)
points(Nn ~ year, data = pop.historyMmeta_ceiling_nomig, pch = 15, cex = 1.2)
points(Ns ~ year, data = pop.historyMmeta_ceiling[1:99,], pch = 9, cex = 1.2)

points(Ns ~ year, data = pop.historyMmeta_ceiling_nomig, pch = 18, cex = 1.2)

abline(v = 2)

## 7 A
# "no sink" = no migration, north only
plot(Nn ~ year, data = pop.historyMmeta_ceiling_nomig, pch = 7,
     xlim = c(0,100),
     ylim = c(0,275),
     ylab = "Total population",
     main = "Figure 7B",
     xaxs="i",
     yaxs="i")

# "No migration" =  Sink present, no migration, N =  north + south, 
points(Ntot ~ year, data = pop.historyMmeta_ceiling_nomig, pch = 18, cex = 1.2)

# "sink present" = sink presnet, with migration (same as figure 7A), N = north + south
points(Ntot ~ year, data = pop.historyMmeta_ceiling[1:99,], pch = 15, cex = 1.2)

Figure 8

plot(Nn ~ year, data = pop.historyMmeta_ceiling[1:99,], pch = 7,
     xlim = c(0,100),
     ylim = c(0,200),
     type = "l",
     ylab = "Number of femalse",
     main = "Figure 8A",
     xaxs="i",
     yaxs="i",
     main = "Figure 8")
abline(h = 100, col = 2)
points(Nnorth ~ year, data = pop.historyMmeta_ceiling[1:99,], type = "l")

plot(Ns ~ year, data = pop.historyMmeta_ceiling[1:99,], pch = 7,
     xlim = c(0,100),
     ylim = c(0,200),
     type = "l",
     ylab = "Number of femalse",
     main = "Figure 8B",
     xaxs="i",
     yaxs="i")
abline(h = 100, col = 2)
points(Nsouth ~ year, data = pop.historyMmeta_ceiling[1:99,], type = "l")

Figure 8c

"Simulations of the population without floaters (i.e., assuming that birds that did not recruit to the population at age 2 died) suggested that the migration of nonterritorial birds (> 1 year old) from the source to the sink population strongly affected the abundance of territorial birds in the sink (Fig. 8C). Including floatersin the model elevated the expected number of breeding birds in the south from 39 to 94 pairs." (pg 318)

From caption: "The effect on sink population sizes of non breeding birds older than 1 yr of age."

Ceiling model with no floaters

# d.n*(1-m) term of element 1,1 replaced with 0
# d.n*m*c term of of element 1,3 replaced with 0
# should element 4,3 be changed?  r.s*(1-m)

Mmeta.ceiling.nofloat <- matrix( data = c(0.00,      y.n*b.n, 0.00,      0.00,
                                  0.00,      s.n,     0.00,      0.00,
                                  0.00,      0.00,    0.00,      y.s*b.s,
                                  r.s*m*c,   0.00,    r.s*(1-m), s.s),
           nrow = 4,
           byrow = T) 

Run ceiling model with no floaters

ceiling_function_no_floaters <- function(){
  x <- rep(NA, 100)
pop.historyMmeta_ceiling <- data.frame(year = 1:length(x),
                                 Fn = x,
                                 Nn = x,
                                 Fs = x,
                                 Ns = x,
                                 Ntot = x)
pop.historyMmeta_ceiling[1,2:5] <- Nmeta.init



  Nmeta.i <- Nmeta.init
for(i in 2:nrow((pop.historyMmeta_ceiling)) ){

  # if(i == 8){
  #   browser()
  # }
  #  



  N.n.current <- Nmeta.i[2] 
  if(N.n.current < T.n){
  Nmeta.i <- Mmeta.Sn.91.nomig%*%Nmeta.i 
  }

  if(N.n.current >= T.n+1){
   Nmeta.i <- Mmeta.ceiling.nofloat%*%Nmeta.i  + I.ceiling


  }

    #Create hard cap at 100
  # if(Nmeta.i[2] > 100){
  #   Nmeta.i[2] <- 100
  # }


  pop.historyMmeta_ceiling[i, 2:5] <- Nmeta.i

}
  return(pop.historyMmeta_ceiling)
}

pop.historyMmeta_ceiling_nofloat <- ceiling_function_no_floaters()

pop.historyMmeta_ceiling_nofloat$Ntot <- apply(pop.historyMmeta_ceiling_nofloat[,c(3,5)],1,sum)

plot(Ns ~ year, data = pop.historyMmeta_ceiling[1:99,], pch = 7,
     xlim = c(0,100),
     ylim = c(0,100),
     ylab = "Breeding pairs",
     main = "Figure 7A",
     xaxs="i",
     yaxs="i")
points(Ns ~ year, data = pop.historyMmeta_ceiling_nofloat)

brouwern/popcycles documentation built on June 10, 2021, 3:23 p.m.