Data: Figure 5 (plus 10 years)

Make data vectors, calculate lambda, and put together dataframe with all necessary data.

census

The census period; an index from 1 to 39 of how many years of data have been collected.

census <- c(1:39)

year t

The year: 1959 to 1997 (Dennis et al use 1959-1987)

year.t   <- c(1959:1997)

Population size

Population size is recorded as the number of females with cubs. This does not represent the actual full population size, but is correlate with it.

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)

Population growth rate: hard-coding example

Population growth rate (lambda) can be calculated as a ratio between the population size in two subsequent years. 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
lambda.59_60

## [1] 1.068182

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

# Access the data
females.N[1 ] #data in braket notation means the Nth of the data set

## [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

females.N.1960/females.N.1959==lambda.59_60 #==force Rstudio to make logic statement whether it is true that those two values equal to one another

## [1] TRUE

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?

females.N[1]

## [1] 44

#lambda.58_59 <- females.N[1]/females.N[ ]  this cannot be done since there is no data for 1958

Population growth rate: vectorized

TASK Briefly describe (1-2 sentence) what this code is doing.

2 Population growth rate (lambda) are calculated in this case. One is from the first to the second year, and another is from the second to the third year.

females.N[2:3]

## [1] 47 46

females.N[1:2]

## [1] 44 47

females.N[2:3]/females.N[1:2] #two lambda values are produced in this case

## [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) #to check the number of data points we have

## [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 What does this next line of code do? Briefly describe in 1 to 2 sentences why I am using length().

calculate all the lambda values throughout the years. Instead of hard coding the number of years of data(39), we use length() here to prevent from mis-counting the number of the data points we have.

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 What does this do? Briefly describe in 1 to 2 sentences what is different about this code chunk from the previous one.

In this code, we do not need to create an object called lens and we can just use length(), and everything is done on a single line of code.

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

females.Ntemp <- females.N[1:10]

Check - are there 10 numbers

length(females.Ntemp)

## [1] 10

TASK What does this do? Briefly describe what the [-1] is doing.

females.Ntemp[-1] prints 9 numbers. It is dropping off the first element of the vector.

females.Ntemp[-1]

## [1] 47 46 44 46 45 46 40 39 39

**TASK* How many lambdas can I calculate using the first 10 years of data? within 10 years of data, we can calculate 9 lambdas.

females.Ntemp[2:10]/females.Ntemp[1:9] #"Negative indexing" allows you to drop a specific element from a vector.

## [1] 1.0681818 0.9787234 0.9565217 1.0454545 0.9782609 1.0222222 0.8695652
## [8] 0.9750000 1.0000000

TASKDrop the the first element

females.Ntemp[-1]

## [1] 47 46 44 46 45 46 40 39 39

TASKDrop the second element

females.Ntemp[-2]

## [1] 44 46 44 46 45 46 40 39 39

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

females.Ntemp[-10]

## [1] 44 47 46 44 46 45 46 40 39

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

females.Ntemp[length(females.Ntemp)]

## [1] 39

TASKHow 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

TASKCalculate the first 9 lambdas.

#females.Ntemp[-1 ] means the data set without the first value and start from the second one
#females.Ntemp[-10] means the data set ends with the 9th value
lambda.i <- females.Ntemp[-1 ]/females.Ntemp[-10]

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

# Explanation are followed:
#All lambdas: 44 47 46 44 46 45 46 40 39 39

# Drop 1st    [-1]47 46 44 46 45 46 40 39 39 

# Drop last    44 47 46 44 46 45 46 40 39 [-10]

# Numerator        47 46 44 46 45 46 40 39 39   
#                 /  /  / ...
# Denominator   44 47 46 44 46 45 46 40 39 

# 47/44 is 1st lambda
# 46/47 is 2nd
# 44/46 is 3rd

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.

What does females.N[-1] do? In females.N data set, it means droping the first element and start from the second one
What does females.N[-length(females.N)? In females.N data set, it means droping the last element

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)]

Finish putting together dataframe

We have made 3 things

an index of the survey, 1 to 39
the raw time series of counts of the bears
a vector of lambda

Create special columns

TASKWhat does this code do? Why do I include NA in the code? (I didn’t cover this in lecture, so just type 1 line - your best guess. “I don’t know” is fine.)

If we had 10 years of data, we’d have 9 transitions.

if we want to put all data - counts (N) and lambdas into a single dataframe, the raw the vector of the raw data (counts) will be longer than the vector of lambdas. So we need to add an extra value to the lambdas to make them the same length. Otherwise, R throws an error.

lambda.i <- c(lambda.i,NA)

TASKCheck the help file; what type of log does log() calculate (I forgot to put this question on the test!)

#log() by default natural logarithms, which is not the log base 10
lambda_log <- log(lambda.i)

#log base 10 would be 
log10(lambda.i)

##  [1]  0.028645181 -0.009340026 -0.019305155  0.019305155 -0.009545318
##  [6]  0.009545318 -0.060697840 -0.010995384  0.000000000  0.032184683
## [11] -0.032184683  0.021719250 -0.010723865 -0.083546051  0.037788561
## [16] -0.024823584  0.059585690 -0.046996563 -0.012589127  0.048304680
## [21] -0.023481096  0.011899223  0.044582133 -0.021719250  0.116505569
## [26] -0.035472318  0.083776998 -0.074633618  0.096910013  0.034762106
## [31]  0.056318363 -0.030382629 -0.025935734 -0.057038501  0.089223184
## [36]  0.063386979  0.087150176  0.000000000           NA

Assemble the dataframe

bear_N <- data.frame(census,
                year.t,
                females.N,
                lambda.i, 
                lambda_log)

TASK

List 3 functions that allow you to examine this dataframe.

1.summary() 2.dim() 3.class()

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,                         # y ~ x
     data = bear_N,                              #data = the dataframe
     type = "b",                                 # b = plots both points and lines
     ylab = "Population index (females + cubs)", #ylab = y axis label
     xlab = "Year")

Bears love to eat trash. Yellowstone closed the last garbage dump in 1970 https://www.yellowstonepark.com/things-to-do/yellowstone-bears-no-longer-get-garbage-treats

Challenge questions

We will cover this all in the next lecture. Feel free to explore this code yourself.

CHALLENGE TASK

Plot a vertical line at 1970. Write a sentence or indicating if you think the population was impacted by this.

Plot lambda

CHALLENGE TASK

Make a histogram of population growth rates
Write a brief sentence describing the shape of these data.

Plot ln(lambda)

CHALLENGE TASK

Make a histogram of natural log population growth rates
Write a brief sentence describing the shape of these data.

Mean of log(lambda)

CHALLENGE TASK

Briefly describe what happens when you delete na.rm = T

mean(bear_N$lambda_lo, na.rm = T)

## [1] 0.02134027

Histogram with mean

In statistics the mean is often represented as the Greek letter “mu”. This can be represented as “u”.

CHALLENGE TASK Save the mean to an object called u

CHALLENGE TASK Make a histogram with the mean plotted on it

Density dependence

CHALLENGE TASK Make a graph that indicates if “density dependence” is occurring.

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

This will be the core topic of next lecture.

Thursday

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)

plot(females.N ~ year.t, data = bear_N, 
     type = "b",
     ylab = "Population index (females + cubs)",
     xlab = "Year")
abline(v = 1970)

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

**First, we create a vector. Then we asked wheather there are any NA values in two ways.

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

length(hat_of_lambdas)  #we checked the length of the vector

## [1] 39

hat_of_lambdas[39]   #check the last value

## [1] NA

hat_of_lambdas[-39]  #get rid of the last value by negative indexing

##  [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)]  #instead of hard-coding, get rid of the last value by negative indexing

##  [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) #get rid of NA in the dataset

##  [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)] #Create a vector of just lambdas, no NAs

hist(hat_of_lambdas) #create a histogram of the lambdas

Randomly Sampling Lambdas

we draw things out of the hat and the sample size is 1, and the sampling is with replacement

sample(x = hat_of_lambdas,  #the input is hat_of_lambda 
       size = 1,
       replace = TRUE)

## [1] 0.9787234

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

Initial Population Size

** population simulations with the final population size of bears in our data set.**

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

Population growth from 1997 to 1998

stimulate the population in 1998

1.22807*99

## [1] 121.5789

lambda_rand.t*N.1997

## [1] 106.6154

N.1998 <- lambda_rand.t*N.1997

Stimulate the population in the future

Steps of stimulating: 1. choose a random lambda value from the hat 2. multiple the lambda value with the population size to get the population size in the next year 3. repeat the step 1 and 2 to get the population size in many years

lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE) #choose a lambda value from a hat
N.1998 <- lambda_rand.t*N.1997              #calculate the population size in the next year

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

Stimulation plot

make a plot with data of population in 7 years and years

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")

Simplify the simulation

#TASK: ADD 1-2 sentences of notes here

#In this code we’ll skip the step of writing the code out many time. Instead, we’ll use a variable called “i” to keep track of how many years into the future we have gone, and initially assign it the value 1. We can then run the code that does the 1st year of simulated population change, then change i to 2. This is an outline of what will become a for() loop.

# Preliminaries
N.1997 <- 99
N.initial <- 99

# plot the initial population 
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)) #make x and y limit 

# The actual stimulation 

N.current <- N.initial

# set initial year
i <- 1             #if we change the i value by hand, then we will get a new graph
  
  # get a random lambda value 
  lambda_rand.i <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # calculate the population in the next year by using the random lambda value 
  N.i <- N.current*lambda_rand.i
  
  # determine the current year
  year.i <- 1997+i
  
  # make a graph 
  points(N.i ~ year.i)

  # store the value of the current population 
  N.current <- N.i

Stimulate a for()loop

#Instead of mannually changing the value of i each time, we use for()loop to handle the process of running the code, updating the year, and re-running the code.

so by using for() loop, i variable can be changed by R automatically from 1 from 1, to 2, to 3, to … the ending.number. ending.number is how many iterations or cycles you want the for() loop to run for. That is, how many time do you want the code to be repeated.

# the initial population in 1997 
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)) #set the x and y limit to allow all the info to appear on the graph

N.current <- N.1997

# make a for() loop
for(i in 1:50){                              # the stimulation run 50 times 
  lambda_rand.i <- sample(x = hat_of_lambdas, #draw a random lambda every time 
                          size = 1,
                          replace = TRUE)
  
  N.i <- N.current*lambda_rand.i        #the function of calculating the new population size
 
  year.i <- 1997+i  # we are making the year to match up with the changing i value in order to match up with the changing population size for each year
  #and the year is increasing every time by 1 year
  
  points(N.i ~ year.i) #making the plot match the stimulated population with the years
 
  N.current <- N.i #update the current population size
}

Set the plotting window to have a 3 by 3 grid.

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

i is a variable and will be changed 50 times in one run to get the population for the next 50 years. And if we run the code many times, we will get different result everytime, since the lambda is randomly picked and different

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

Plot a 10 x 10 grid of population

make a 10 x 10 grid which only plot one simulation.

par(mfrow = c(10,10), mar = c(0,0,0,0), xaxt = "n", yaxt = "n") #set up the plotting panel

run stimulation

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550), cex = 0.5)
abline(h = N.1997)
abline(h = 0, col = "red")
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, cex = 0.5)
  
  N.current <- N.t
}

ADD TITLE HERE

Add 1-2 sentence of notes here

## ADD TITLE HERE
num.years <- 50

## ADD TITLE HERE
N.vector <- rep(NA, num.years)

## ADD TITLE HERE
N.vector[1] <- N.1997

t<- 1
  # print the year
  print(t)

## [1] 1

  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t

ADD TITLE HERE

Add 1-2 sentence of notes here

## ADD TITLE HERE
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD TITLE HERE
#don't uncomment print() in notebook mode!
for(t in 1:num.years){
  # print the year
  #print(t)
  
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
  
}

ADD TITLE HERE

par(mfrow = c(1,1), mar = c(3,3,3,3), xaxt = "s", yaxt = "s")

df <- data.frame(N.vector, year = seq(1997, 1997+50))

plot(N.vector ~ year, data = df)

plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

ADD TITLE HERE

ADD 2-3 sentences of notes here

## ADD TITLE
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD TITLE
lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)


## ADD TITLE
t <- 1


  print(t)

## [1] 1

  lambda_rand.t <- lambdas_vector[t]
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
  t <- t+1
  N.vector[1:t]

## [1] 99.00000 96.84783

ADD TITLE HERE

## ADD NOTES
num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

## ADD NOTES
lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)


## ADD NOTES
for(t in 1:num.years){
  
  ## ADD NOTES
  lambda_rand.t <- lambdas_vector[t]
  
  ## ADD NOTES
  N.t <- N.vector[t]*lambda_rand.t
  
  ## ADD NOTES
  N.vector[t+1] <- N.t
  
}

## ADD NOTES
df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

How can I shorten the code?

lambdas_vector[t]
N.t

num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)

for(t in 1:num.years){
  lambda_rand.t <- lambdas_vector[t]
  
  N.t <- N.vector[t]*lambda_rand.t
  
  N.vector[t+1] <- N.t
  
}

df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

Shortened code

num.years <- 50
N.vector <- rep(NA, num.years)
N.vector[1] <- N.1997

lambdas_vector <- sample(x = hat_of_lambdas, size = num.years, replace = TRUE)

for(t in 1:num.years){
  N.vector[t+1] <- N.vector[t]*lambdas_vector[t]
  
}

df <- data.frame(N.vector, year = seq(1997, 1997+50))
plot(N.vector ~ year, data = df, ylim = c(0, max(df$N.vector)))

Profolio 1

Caiyi Wu

9/24/2020

Data: Figure 5 (plus 10 years)

census

year t

Population size

Population growth rate: hard-coding example

Population growth rate: vectorized

Negative indexing

Calcualte lambdas for all data

Finish putting together dataframe

Create special columns

Assemble the dataframe

Examing the population growth rates

Plotting the raw data

Challenge questions

Plot lambda

Plot ln(lambda)

Mean of log(lambda)

Histogram with mean

Density dependence

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

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

Randomly Sampling Lambdas

Initial Population Size

Population growth from 1997 to 1998

Stimulate the population in the future

Stimulation plot

Simplify the simulation

so by using for() loop, i variable can be changed by R automatically from 1 from 1, to 2, to 3, to … the ending.number. ending.number is how many iterations or cycles you want the for() loop to run for. That is, how many time do you want the code to be repeated.

Plot a 10 x 10 grid of population

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

ADD TITLE HERE

Shortened code

Parameteric sampling