Portfolio 1: PVA
Nathan Brouwer - Manipulated by Isabella Peterson
9/20/2020
Data: Figure 5 (plus 10 years)
Make data vectors, calculate lambda, and put together data frame with all necessary data.
Census
The census period; an index from 1 to 39 of how many years of data have been collected. The following are different variations of creating a object contain numbers 1-39 representing the years of data collect during the census period.
census <- 1:39
census <- c(1:39)
census <- seq (1,39)
census <- seq(1, 39, by= 1) year t
The year: 1959 to 1997 (Dennis et al use 1959-1987)
year.t <- 1959:1997
year.t <- seq(1959,1997) # Another variation of the code above just using the seq() function to create the list of years from 1959 to 1997Population size (N)
Population size is recorded as the number of females with cubs from the census mentioned above.
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: example
Population growth rate is (or can be expressed as) a ratio of population sizes(N).
Enter the population size for each year
females.N.1959 <- 44 # Manually hard coding the number of females with cubs to the corresponding year in the census; same for line of code below
females.N.1960 <- 47Calculate the ratio of the 2 population sizes
lambda.59_60 <- females.N.1960/females.N.1959 # Hard coding the population growth rate with the populations sizes (N) from 1959-1960Access the population sizes by using bracket notation or vector indexing rather than hard coding
# Access the data using indexing
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.068182Calculate 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[NA] # Can't be done due to there being no data for 1958, census period began in 1959Population growth rate: vectorized
TASK Briefly describe (1-2 sentences) what this code is doing.
Prints out certain vector elements using vector indexing or bracket notation.
females.N[2:3] # Prints 2 and 3 from the female.N vector
#> [1] 47 46
females.N[1:2] # Prints the values of 1 and 2 from the female.N vector
#> [1] 44 47
females.N[2:3]/females.N[1:2] # Gives the two lambda values or populations growth rate
#> [1] 1.0681818 0.9787234This is similar to 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] # Sames as above just using a vector indexing to get a range of elements instead of a certain element
#> [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.0000000TASK What does this do? Briefly describe in 1 to 2 sentences why I am using length().
It uses length to get number of elements within the females.N vector aka 39 and assigning it to object len to limit hard coding 39 into chunk below.
len <- length(females.N) # Taking the length of the Vector (39) and its the same as above without hard coding
females.N[2:len]/females.N[1:len-1] # Similar to above, instead of using numeric values to get a range with vector indexing a object or variable is used
#> [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.0000000TASK What does this do? Briefly describe in 1 to 2 sentences what is different about this code chunk from the previous one.
Its a cleaner and more generalized form of the code chuck above. Instead of assigning the function length() to a object to be used within vector indexing the function itself is placed within the bracket notation.
females.N[2:length(females.N)]/females.N[1:length(females.N)-1] # Same as code above just generalized, the length is put inside the function
#> [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.0000000Negative 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] 10TASK What does this do? Briefly describe what the [-1] is doing.
[-1] is cutting of the first element of the vector through negative indexing
females.Ntemp[-1] # Cuts off the first element of the vector
#> [1] 47 46 44 46 45 46 40 39 39TASK How many lambdas can I calculate using the first 10 years of data?
females.Ntemp[2:10]/females.Ntemp[1:9] # Calculates 9 lambdas due to dropping of the first element of in the first females.Ntemp and the dropping of the last element in the second one
#> [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 39TASK Drop the second element
females.Ntemp[-2]
#> [1] 44 46 44 46 45 46 40 39 39TASK How do you drop the 10th element? Type in the code below.
females.Ntemp[-10]
#> [1] 44 47 46 44 46 45 46 40 39TASK How do you access the last element? Do this in a general way without hard-coding.
females.Ntemp[length(females.Ntemp)] # The length of the vector is also the number of the last element within the vector, in this case 10
#> [1] 39TASK 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[-10]
#> [1] 44 47 46 44 46 45 46 40 39
females.Ntemp[-length(females.Ntemp)]
#> [1] 44 47 46 44 46 45 46 40 39TASK 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)]Calculate 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?
Drops the first element of the vector.
What does females.N[-length(females.N)?
Drops the last element of the vector.
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 data frame
Create special columns
TASK What 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.)
Removes NA (Not available) values from lambda.i
lambda.i <- c(lambda.i, NA) # Removes NA from the vectorTASK Check the help file; what type of log does log() calculate (I forgot to put this question on the test!)
lambda_log <- log(lambda.i) #ln or natural logarithmsAssemble 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 data frame.
- summary() # Gives mean, median, mode, min and max from data frames
- head() # Shows top 6 rows of data frame
- tail() # Shows bottom 6 rows of data frame
- dim() # Gives the dimensions of the data frame
Examining 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(x = year.t,
xlab = "Year",
y= females.N,
ylab = "Populations Index (Females + Cubs)",
type = "b")
# Or
plot(females.N~year.t,
data= bear_N,
xlab = "Year",
ylab = "Populations Index (Females + Cubs)",
type= "b")
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
plot(females.N ~ year.t, data = bear_N,
type = "b",
ylab = "Population index (females + cubs)",
xlab = "Year")
abline(v = 1970)
abline(v = 1987, col="red") # Original extinction cut off
How do we determine if a population is likely to go extinct?
To model population dynamics we randomly pull population growth rates of lambda out of a hat.
hat_of_lambdas <- bear_N$lambda.i # Lambda is the population growth rate resulting from ratio of populations sizes(N). It is also a parameter that is variable
# hat_of lambdas is the object holding populations growth rates that will be drawn from random for modeling useis.na(hat_of_lambdas) # Searches for NA or missing data store in the object
#> [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 TRUEany(is.na(hat_of_lambdas) == TRUE) # If there are any NA a TRUE value should pop up
#> [1] TRUEDrop the NA
length(hat_of_lambdas)
#> [1] 39
hat_of_lambdas[39] # Element 39 is NA
#> [1] NA
hat_of_lambdas[-39] # Drops the 39t element
#> [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)] # More stable than negative indexing above
#> [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.0000000na.omit(hat_of_lambdas) # Removes NA from object
#> [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)] # Creates vector that drops the 39 element or the NA in the vectorHistogram of hat_of_lambdas
hist(hat_of_lambdas)
Random Sampling of Lambdas (population growth rates)
Takes a sample out of the the hat_of_lambdas vector for the populations growth rate for one year.
# Pulled a lambda value from the hat_of_lambdas vector
sample(x = hat_of_lambdas,
size = 1,replace = TRUE) #Populations growth for number(size) of years
#> [1] 0.975
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE) # Pulled lambda value is saved to vector lambda_rand.tCheck the bear.N data frame for intial population size
Observe the bear_N data frame through the following functions
head(bear_N) #Shows top 6 values of data frame
#> 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) # Shows bottom 6 values of data frame
#> 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
summary(bear_N) # Shows a summary of the data frame including its mean, median, mode, max and min
#> census year.t females.N lambda.i
#> Min. : 1.0 Min. :1959 Min. :33.00 Min. :0.8250
#> 1st Qu.:10.5 1st Qu.:1968 1st Qu.:39.00 1st Qu.:0.9452
#> Median :20.0 Median :1978 Median :44.00 Median :1.0000
#> Mean :20.0 Mean :1978 Mean :49.79 Mean :1.0281
#> 3rd Qu.:29.5 3rd Qu.:1988 3rd Qu.:57.00 3rd Qu.:1.1038
#> Max. :39.0 Max. :1997 Max. :99.00 Max. :1.3077
#> NA's :1
#> lambda_log
#> Min. :-0.19237
#> 1st Qu.:-0.05639
#> Median : 0.00000
#> Mean : 0.02134
#> 3rd Qu.: 0.09874
#> Max. : 0.26826
#> NA's :1
dim(bear_N) # Gives dimensions of data frame
#> [1] 39 5
N.1997 <- 99 # Assign 99 to N.1997 to get the starting pointOne Round of Population Simulation Through Hard Coding
Takes the sampled lambda and multiples it by the population size (N) from 1997 and assign it to object N.1998
1.22807*99
#> [1] 121.5789
lambda_rand.t*N.1997
#> [1] 104.0769
N.1998 <- lambda_rand.t*N.1997 Simulation the Hard Way
Manually pulling sample out of hat_of_lambdas (Populations growth rates) then multiplying it by the previous year’s populations and assigning to the currents years object.
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE) # Pulls out the lambda multiplies it by populations of previous year than puts populations size in new object
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.2004Plot Population
Makes a sequence of years for x-axis and combines the simulated population sizes into a object. Combine the two vectors into a data frame and makes a scatter plot showing the relationship between the simulated populations growth rate 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")
for() loop that plots a single point
Makes a for() loop to plot a single population size data point onto a scatter plot
# Initial Conditions(Populations Size)
N.1997 <- 99
N.initial <- 99
# Explore xlim = argument
plot(N.1997 ~ c(1997))
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50)) # xlim sets the size of the x-axis
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550)) # xlim sets the size of the x-axis, and ylim sets limit of y-axis
# For loop the hard way
N.current <- N.initial
# Start of for() loop
t <- 1 # t=time; t <- 1 is assigning # year to time
# Pulls lambda from hat
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
# Determine Population size
N.t <- N.current*lambda_rand.t
# Adds year to starting year
year.t <- 1997+t
# Plots the population
points(N.t ~ year.t) # Adds point to existing plot; Doesn't reset plot
# Updates N.current
N.current <- N.t50 point plotting for() loop
for () loop that pulls 50 samples from hat_of_lambdas and plots separate points onto the plot
# Makes a new plot with starting population size
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
# Starts at 1997 population size
N.current <- N.1997
# Changes t to 1:50 after running rest of code
for(t in 1:50){
# Pulls random lambda from the hat
lambda_rand.t <- sample(x = hat_of_lambdas,
size = 1,
replace = TRUE)
# Multiples 1997's population size by the random lambda(populations growth rate) and stores result into N.t
N.t <- N.current*lambda_rand.t
# Adds t (1:50) values to the starting year and stores it into year.t
year.t <- 1997+t
# Plots the point of the relations between the N.t and year.t
points(N.t ~ year.t)
# Resets N.t values 1997 for beginning of the loop
N.current <- N.t
}
Setting parameters for R plotting; Allows for plots to show up in a 3x3 grid within the plot viewer
par(mfrow = c(3,3), mar = c(1,1,1,1))for() loop for 50 data points with axis limits
Same as the for() loop above, but it makes a plot with limits on the x and y axis
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
}