##Data

census <- 1:39
year.t <- 1959:1997
females.N <- c(44,47,46,44,46,
               45,46,40,39,39,
               42,39,41,40,33,
               36,34,39,35,34,
               38,36,37,41,39,
               51,47,57,48,60,
               65,74,69,65,57,
               70,81,99,99)
lambda.i <- females.N[-1]/females.N[-length(females.N)]
lambda.i <- c(lambda.i,NA)
lambda_log <- log(lambda.i)
bear_N <- data.frame(census, year.t, females.N, lambda.i, lambda_log)

Population growth rate: example

Population growth rate is…

Enter the population size for each year

females.N.1959 <-44 
females.N.1960 <-47

Calculate the ratio of the 2 population sizes

lambda.59_60 <- females.N.1960/females.N.1959

Access the population sizes by using bracket notation rather than hard coding

# Access the data
females.N[1 ]
#> [1] 44
females.N[ 2]
#> [1] 47

# store in objects
females.N.1959 <- females.N[1 ]
females.N.1960 <- females.N[2 ]

# confirm the output
females.N.1960/females.N.1959
#> [1] 1.068182

Calculate lambda using bracket notation

lambda.59_60 <- females.N[2 ]/females.N[1 ]

The first year of data is 1959. What is lambda for 1958 to 1959? Null will occur because the year 1958 isnt in the data

#females.N[1]
#lambda.58_59 <- females.N[1]/females.N[ ]

Population growth rate: vectorized

TASK

To create vectors that contain the bear population data in order to carry out a lambda calculation with less hardcoding. this allows us to carry out the same equation for multiple datas instead of hand writing all the data in

females.N[2:3]
#> [1] 47 46
females.N[1:2]
#> [1] 44 47

females.N[2:3]/females.N[1:2]
#> [1] 1.0681818 0.9787234

This is similar t the previous code chunk, just using all of the data (no need to describe)

length(females.N)
#> [1] 39
females.N[2:39]/females.N[1:38]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000

TASK the first code line creates a variable (len) that contains the length of the female population (the census), the 2nd line calculates the lamda using vectors. instead of hard coding the number of censuses we can use len-1 in order to subtract from the total census (39) into 38.

len <- length(females.N)
females.N[2:len]/females.N[1:len-1]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000

TASK This code calculates all the lambdas for the bear population using 1 line. instead of assigning the length of the vector into a variable we directly code it into the formula (same with len-1 which will become length(females.N)-1 )


females.N[2:length(females.N)]/females.N[1:length(females.N)-1]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000

Negative indexing

Make a short vector to play with; first 10 years

females.N[1:10]
#>  [1] 44 47 46 44 46 45 46 40 39 39
#>  [1] 44 47 46 44 46 45 46 40 39 39
females.N[seq(1,10)]
#>  [1] 44 47 46 44 46 45 46 40 39 39
#>  [1] 44 47 46 44 46 45 46 40 39 39

females.Ntemp <- females.N[seq(1,10)]

Check - are there 10 numbers

females.Ntemp
#>  [1] 44 47 46 44 46 45 46 40 39 39

TASK

Removing the first index of a vector

females.Ntemp[-1]
#> [1] 47 46 44 46 45 46 40 39 39
#> [1] 47 46 44 46 45 46 40 39 39

TASK How many lambdas can I calculate using the first 10 years of data? answer: 9

females.Ntemp[2:10]/females.Ntemp[1:9]
#> [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#> [8] 0.9750000 1.0000000
#> [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#> [8] 0.9750000 1.0000000

“Negative indexing” allows you to drop a specific element from a vector.

TASK Drop the the first element

females.Ntemp[-1]
#> [1] 47 46 44 46 45 46 40 39 39

TASK Drop the second element

females.Ntemp[-2]
#> [1] 44 46 44 46 45 46 40 39 39

TASK

How do you drop the 10th element? Type in the code below.

females.Ntemp[-10]
#> [1] 44 47 46 44 46 45 46 40 39

TASK How do you access the last element? Do this in a general way without hard-coding.

females.Ntemp[length(females.Ntemp)]
#> [1] 39

TASK How do DROP the last element? Do this in a general way without hard-coding. By general, I mean in a way that if the length of the vector females.Ntemp changed the code would still drop the correct element.

females.Ntemp[-length(females.Ntemp)]
#> [1] 44 47 46 44 46 45 46 40 39

TASK Calculate the first 9 lambdas.

lambda.i <- females.Ntemp[-1]/females.Ntemp[-10]

Converting between these 2 code chunks would be a good test question : )

lambda.i <- females.Ntemp[-1]/females.Ntemp[-length(females.Ntemp)]

Calcualte lambdas for all data

TASK

Below each bulleted line describe what the parts of the code do. Run the code to test it.

TASK Calculate lambdas for all of the data

females.N[-1]
#>  [1] 47 46 44 46 45 46 40 39 39 42 39 41 40 33 36 34 39 35 34 38 36 37 41 39 51
#> [26] 47 57 48 60 65 74 69 65 57 70 81 99 99
females.N[-length(females.N)]
#>  [1] 44 47 46 44 46 45 46 40 39 39 42 39 41 40 33 36 34 39 35 34 38 36 37 41 39
#> [26] 51 47 57 48 60 65 74 69 65 57 70 81 99

lambda.i <- females.N[ -1]/females.N[-length(females.N) ]

TASK

List 3 functions that allow you to examine this dataframe.

1.head 2.tail 3.summary

Examing the population growth rates

Plotting the raw data

TASK

  • Plot a time series graph of the number of bears (y) versus time (x)
  • Label the y axis “Population index (females + cubs)”
  • Label the x axis “Year”
  • Change the plot to type = “b” so that both points and dots are shown.
plot(females.N ~ year.t, data = bear_N, 
     type = "b",
     ylab = "Population index (females + cubs)",
     xlab = "Year")

## How do we determine if a population is likely to go extinct?

to model the population dynamics we randomly pull population growth rates out of a hat

hat_of_lambdas <- bear_N$lambda.i
is.na(hat_of_lambdas)
#>  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [37] FALSE FALSE  TRUE
any(is.na(hat_of_lambdas) == TRUE)
#> [1] TRUE

Drop the NA

length(hat_of_lambdas)
#> [1] 39
hat_of_lambdas[39]
#> [1] NA
hat_of_lambdas[-39]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000
hat_of_lambdas[-length(hat_of_lambdas)]
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000
na.omit(hat_of_lambdas)
#>  [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
#>  [8] 0.9750000 1.0000000 1.0769231 0.9285714 1.0512821 0.9756098 0.8250000
#> [15] 1.0909091 0.9444444 1.1470588 0.8974359 0.9714286 1.1176471 0.9473684
#> [22] 1.0277778 1.1081081 0.9512195 1.3076923 0.9215686 1.2127660 0.8421053
#> [29] 1.2500000 1.0833333 1.1384615 0.9324324 0.9420290 0.8769231 1.2280702
#> [36] 1.1571429 1.2222222 1.0000000
#> attr(,"na.action")
#> [1] 39
#> attr(,"class")
#> [1] "omit"
hat_of_lambdas <- hat_of_lambdas[-length(hat_of_lambdas)]
hist(hat_of_lambdas)

Random sampling of lambda

in this chunk of code we pull random lambda from the hat and use the 2nd line of code (316) in order to save that random lambda we pulled into a object.

# pulls out a random lambda from the hat
sample(x = hat_of_lambdas, size = 1,
       replace = TRUE) #replace = true is the same as turning on recycling
#> [1] 0.9714286
#save the lambda to an object/variable
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)

starting initial population size

the code below examins the top and bottom of the dataframe (head and tail) while the last line assigns the number of bears we start the simulation at

head(bear_N)
#>   census year.t females.N  lambda.i  lambda_log
#> 1      1   1959        44 1.0681818  0.06595797
#> 2      2   1960        47 0.9787234 -0.02150621
#> 3      3   1961        46 0.9565217 -0.04445176
#> 4      4   1962        44 1.0454545  0.04445176
#> 5      5   1963        46 0.9782609 -0.02197891
#> 6      6   1964        45 1.0222222  0.02197891
tail(bear_N)
#>    census year.t females.N  lambda.i lambda_log
#> 34     34   1992        65 0.8769231 -0.1313360
#> 35     35   1993        57 1.2280702  0.2054440
#> 36     36   1994        70 1.1571429  0.1459539
#> 37     37   1995        81 1.2222222  0.2006707
#> 38     38   1996        99 1.0000000  0.0000000
#> 39     39   1997        99        NA         NA
N.1997 <- 99 #last available population

One round of population simulation

calculates the number of bears after 1 year of simulation

1.22807*99
#> [1] 121.5789
lambda_rand.t*N.1997
#> [1] 96.525
N.1998 <- lambda_rand.t*N.1997

Simulation the hard way ;C

Calculates the simulated number of bears for the following years manually as if we are creating a table, somewhat like excel

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1998 <- lambda_rand.t*N.1997

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.1999 <- lambda_rand.t*N.1998

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2000 <- lambda_rand.t*N.1999

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2001 <- lambda_rand.t*N.2000

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2002 <- lambda_rand.t*N.2001

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2003 <- lambda_rand.t*N.2002

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2004 <- lambda_rand.t*N.2003

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
N.2005 <- lambda_rand.t*N.2004

All rounds of population simulation automated

puts the amounts of years into a vector and simulate the future generation bears, lastly we make a graph to represent our new simulation

year <- seq(1997, 2004)
N.rand <- c(N.1998,N.1999,N.2000,N.2001,N.2002,N.2003,N.2004,N.2005)
df.rand <- data.frame(N.rand, year)
plot(N.rand ~ year, data = df.rand, type = "b")

##For loops

we plot the simulation of a new year and move around where they are located on a graph which is then plotted.

#initial population

N.1997 <- 99
N.initial <- 99

# Explore xlimit arguments which sets where we start the population and ploting
plot(N.1997 ~ c(1997))

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50))

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))


#for loop hard way
N.current <- N.initial

# this is where the forloop would be located
t <- 1
  
  # Grab a lambda
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # determin new population size
  N.t <- N.current*lambda_rand.t
  
  # tell R which year the growth is
  year.t <- 1997+t
  
  # adds points to a plot
  points(N.t ~ year.t)

  
  # update N.current
  N.current <- N.t

we plot 1 round of simulation for 50 years using a for loop


# create new plot
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))

N.current <- N.1997

# creates for loop and set its conditions
for(t in 1:50)
  {
  
  # pick a random lambda out of the hat
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # update the population size
  N.t <- N.current*lambda_rand.t
  
  # updates the year
  year.t <- 1997+t
  
  # adds a new point to the graph
  points(N.t ~ year.t)
  
  # updates the current population
  N.current <- N.t
}

Goofy r plotting code/magic

par(mfrow = c(3,3), mar = c(1,1,1,1))

simulates the population growth for the next 50 years, but we can make grids with this due to Dr Brouwer’s magic abilities

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
N.current <- N.1997
for(t in 1:50){
  
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  N.t <- N.current*lambda_rand.t
  
  year.t <- 1997+t
  
  points(N.t ~ year.t)
  
  N.current <- N.t
}