PVA Assignment 1: for 9/15/2020 and 9/17/2020

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 <- 1:39 # The ':' operator can be used in some cases as a shorthand for loop sequence
census <- seq(1, 39)

year t

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

year.t <- 1959:1997

Population size

Population size is recorded as the number of females with cubs each year. Females are used to track bear and the grizzly population, because they factor in the birthing rates for the total bear population and they can count cubs given birth to each year.

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 often denoted as lambda. The population growth rate typically takes the change of individuals of species in year + 1 / individuals of species in year This can be donated as N(t+1)/N(t). In this case we are comparing the bear populations in 1960 (47) and 1959 (44).

Enter the population size for each year

females.N.1959 <- females.N[1] # Population size: 44 
females.N.1960 <- females.N[1 + 1] # Population size: 47, I also used the 1 + 1 to say that it's the following year.

females.N.1960

## [1] 47

females.N.1959

## [1] 44

Calculate the ratio of the 2 population sizes

lambda.59_60 <- females.N.1960/females.N.1959 # This is equivalent to 47/44 => 1.068182


lambda.59_60 # This means the annual growth rate from 1959 to 1960 was about 7%

## [1] 1.068182

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

# lambda is 1.068182

Calculate lambda using bracket notation

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

## [1] NA

# lambda is 1.068182

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

# Since lambda would need the bear population from 1958, lambda would not exist in this case.

Population growth rate: vectorized

females.N[2:3] This is getting the population data of bears from years 1960 to 1961 females.N[1:2] This is getting the population data of bears from years 1959 to 1960 females.N[2:3]/females.N[1:2] This is taking the data of bears from year 1960 and 1959 then dividing them to get the lambda.58_59 value above, it’s then doing the same operation but will the 1961 and 1960 year.

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

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

len <- length(females.N) This is getting the length as a numeric value of females.N vector, which is 99. females.N[2:len]/females.N[1:len-1] This is equivalent to females.N[2:39]/females.N[1:38], but using a less hard coded method.

What does this do? Briefly describe in 1 to 2 sentences why I am using length().

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

This is equivalent to ^^ code. However instead of using the object len, we are directly inserting the length(females.N). What does this do? Briefly describe in 1 to 2 sentences what is different about this code chunk from the previous one.

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] # Hard Coded first ten years

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

females.Ntemp <- females.N[-c(11 : length(females.N))]

females.N[1:10] == females.Ntemp # The options above are the same in output

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Check - are there 10 numbers

length(females.Ntemp) == 10

## [1] TRUE

# Prints TRUE if length of females.Ntemp is equal to 10

-1 does NOT point to the last index, this strange, but it removes the first element in this vector, and then gets items 2:length(females.Ntemp)

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

females.Ntemp[-1]

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

length(females.Ntemp[2:length(females.Ntemp)]) == length(females.Ntemp[-1])

## [1] TRUE

The first 10 years will generate 9 lambdas, since we cannot get a lambda for 1959 and 1958. How many lambdas can I calculate using the first 10 years of data?

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

length(females.Ntemp[2:10]/females.Ntemp[1:9]) == 9

## [1] TRUE

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

Dropping the first element can be done by using -1.

Drop the the first element

females.Ntemp[-1]

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

Dropping the second element can be done using -2

Drop the second element

females.Ntemp[-2]

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

Dropping the 10th element can be done by using -10

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

We can access the last element with the length(), by stating length(females.Ntemp).

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

females.Ntemp[length(females.Ntemp)] # 39 in this case

## [1] 39

This can be done using the same logic as above, with with using negative indexing. We can either do this in two different ways below, by using a vector or not 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

females.Ntemp[-c(length(females.Ntemp))]

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

females.Ntemp[-length(females.Ntemp)] == females.Ntemp[-c(length(females.Ntemp))]

## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

We can do this by dividing the first 10 years + 1 year / the last 10 years, by negative indexing The above models the N(t + 1)/ N(t) equation. Calculate the first 9 lambdas.

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

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

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

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

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

Calcualte lambdas for all data

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

What does females.N[-1] do? This removes the first item and gets the vector from 2: length of the vector females.N[-1]
What does females.N[-length(females.N)?

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) ]
lambda.i

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

Finish putting together dataframe

Create special columns

I’m not to sure but, my guess is that it could “break” the code, since NA cannot be used in mathematical operations typically, maybe it’s used as a placeholder.

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

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

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

log() takes the natural log (ln) of a function. If we want regular log base 10, I think it’s log10()

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

lambda_log <- log(2.72) # Should be about 1, since we our comparing e^e
lambda_log <- log(lambda.i)

Assemble the dataframe

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

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
## 7       7   1965        46 0.8695652 -0.13976194
## 8       8   1966        40 0.9750000 -0.02531781
## 9       9   1967        39 1.0000000  0.00000000
## 10     10   1968        39 1.0769231  0.07410797
## 11     11   1969        42 0.9285714 -0.07410797
## 12     12   1970        39 1.0512821  0.05001042
## 13     13   1971        41 0.9756098 -0.02469261
## 14     14   1972        40 0.8250000 -0.19237189
## 15     15   1973        33 1.0909091  0.08701138
## 16     16   1974        36 0.9444444 -0.05715841
## 17     17   1975        34 1.1470588  0.13720112
## 18     18   1976        39 0.8974359 -0.10821358
## 19     19   1977        35 0.9714286 -0.02898754
## 20     20   1978        34 1.1176471  0.11122564
## 21     21   1979        38 0.9473684 -0.05406722
## 22     22   1980        36 1.0277778  0.02739897
## 23     23   1981        37 1.1081081  0.10265415
## 24     24   1982        41 0.9512195 -0.05001042
## 25     25   1983        39 1.3076923  0.26826399
## 26     26   1984        51 0.9215686 -0.08167803
## 27     27   1985        47 1.2127660  0.19290367
## 28     28   1986        57 0.8421053 -0.17185026
## 29     29   1987        48 1.2500000  0.22314355
## 30     30   1988        60 1.0833333  0.08004271
## 31     31   1989        65 1.1384615  0.12967782
## 32     32   1990        74 0.9324324 -0.06995859
## 33     33   1991        69 0.9420290 -0.05971923
## 34     34   1992        65 0.8769231 -0.13133600
## 35     35   1993        57 1.2280702  0.20544397
## 36     36   1994        70 1.1571429  0.14595391
## 37     37   1995        81 1.2222222  0.20067070
## 38     38   1996        99 1.0000000  0.00000000
## 39     39   1997        99        NA          NA

dim(bear_N)

## [1] 39  5

str(bear_N)

## 'data.frame':    39 obs. of  5 variables:
##  $ census    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ year.t    : int  1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 ...
##  $ females.N : num  44 47 46 44 46 45 46 40 39 39 ...
##  $ lambda.i  : num  1.068 0.979 0.957 1.045 0.978 ...
##  $ lambda_log: num  0.066 -0.0215 -0.0445 0.0445 -0.022 ...

colnames(bear_N)

## [1] "census"     "year.t"     "females.N"  "lambda.i"   "lambda_log"

rownames(bear_N)

##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
## [16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
## [31] "31" "32" "33" "34" "35" "36" "37" "38" "39"

List 3 functions that allow you to examine this dataframe.

dim() - Shows the dimensions of the data frame by rows and cols, in that order.
str() - Shows the overall structure of the data set, really good for analyzing columns output in the rows
colnames() or rownames() - gets the name of each cols or rows depending on which is used.

Examing the population growth rates

Plotting the raw data

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, xlab = "Year", ylab = "Population index (females + cubs)", 
     main = "Number of Bears versus Time", 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

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

This will be the core topic of next lecture.

I think we could do a simulation of the bear population using the sample(). We could use the lambda values and then try to run simulations projecting the bear populations randomly increasing of decreasing by each year.

PVA for 9/17/2020

Data we will be working with for 9/17/2020 Lecture

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?

We can try to predict the future population of bears by using simulation. In simulation we will randomly choose a population growth rate (lambda) and multiply by the last year of data(1997) and keep looping with the previous year and multiplying with a random lambda until we get to the year we want to simulate too.

hat_of_lambdas <- bear_N$lambda.i
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

is.na(hat_of_lambdas) # Are there any NA present in this vector?

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

# TRUE

any(is.na(hat_of_lambdas) == TRUE) # This will be TRUE, since there is an NA

## [1] TRUE

Drop the NA We can drop the NA by dropping the last element from our vector w/ negative indexing.

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"

We can re-update the hat of lambdas by getting rid of the last index which has the NA.

hat_of_lambdas <- hat_of_lambdas[-length(hat_of_lambdas)]
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

This is a histogram of the hat of lambda values. The histogram does appear to be normal in distribution.

hist(hat_of_lambdas)

Let’s Pick a Lambda out of the Hat!

Using the sample() function we are going to use the vector hat_of_lambdas as our sample size, pick 1 lambda value and then if we were going to pick another lambda we would recycle our last value.

# sample(sample as x, size is how many times we will pick an item, 
# replace is do we replace if we select another item?)
sample(x = hat_of_lambdas, size = 1,replace = TRUE)

## [1] 0.9324324

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

Head and Tails with data frames in R

The head function will produce some of the first rows of items in bear_N. The tail function will produce some of the last rows of items in bear_N.

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

Predicting Growth Rate in 1998, by using a random lambda.

We are going to select a random lambda with the object lambda_rand.t then we are going to multiply that by the last year’s recorded population to try to predict the number of bears for 1998.

lambda_rand.t*N.1997

## [1] 99

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

Predicting future data the ha(R)d way

Essentially we are going to do what we did above, but with the value we get from N.1998, we are then going to do the same calculation with lambda, lambda *random bear population prior to get the next years random bear population.

Get the actual bear data from 1997, then choose random lambda and multiply to get simulated bear population in 1998. using the random simulated bear population from 1998, then take another random lambda and multiply by 1998 to get 1999, this process is continued to year 2005.

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

Implementing simulated population of bears from year 1997 - 2005

The year sequence and the random populations of years 1998 - 2005 produced above are used to create a data frame. We then use the data frame along with our x-axis: years and y-axis N.rand to plot the population of bears per each year.

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

Let’s make some plot limits

Using xlim and ylim to limit the amount of years to 2047 and the total population of a given year to 550.

# Our initial population of bears

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

# Plot of population of bears vs. year 1997 w/ our x and y limits, this will only plot the one point
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))

# Setting our current population of bears to the initial bears or 99
N.current <- N.initial

# Setting our year iterator t to 1.
t <- 1
  
  # Getting a random lambda from the pop. of bears recorded.
  lambda_rand.t <- sample(x = hat_of_lambdas, size = 1,replace = TRUE)
  
  # Generating the next year with our current bear population 99 and random lambda above.
  N.t <- N.current*lambda_rand.t
  
  # Increases the year by 1, or the first simulated year of data for 1998.
  year.t <- 1997+t
  
  # gets points from current year population versus a given year
  # points is often used to add more points into a plot
  # In this case it adds another set of a point to the graph
  points(N.t ~ year.t)

  # Updates the current population of bears to N.t, which is 
  # simulated population of bears from 1998.

  # We need to use for loop to get all the points in graph
  N.current <- N.t

We can go through the same process above. Instead of recording each either through objects, we can use the for loop to loop through years 1997 - 2047 using t to track the year.

# Starts plot with the first point of y: N.1997 <- 99 and x: c(1997) <- 1997
plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
# abline(v = 2047) # Added the vertical lines to see the domain clearer by years 
# abline(v = 1997) 

N.current <- N.1997

# Uses for loop to loop through the 50 years of simulated data
for(t in 1:50){
  
  # Generates the random lambda value to use
  lambda_rand.t <- sample(x = hat_of_lambdas, 
                          size = 1,
                          replace = TRUE)
  
  # Generates the simulated year of population of bears
  N.t <- N.current*lambda_rand.t
  
  # This stores our years, starts at 1997, then increases
  # by 1 year until we reach 2047.
  year.t <- 1997+t
  
  # adds new points to the plot, this is from our simulated 
  # population of bears by the year we simulated the data from.
  
  points(N.t ~ year.t)
  
  # Update current population of bears, so that N.t is a new
  # simulated number of bears 
  N.current <- N.t
}

Par sets graphical parameters. In this case, the mfrow makes it so we can look at 3 x 3 graphs in the plot viewer in the code chunk below. While mar, sets the margin size.

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

This is the same code as two chunks above. The only difference is now with the par(), we are looking at 3 x 3 simulations with those margin specs listed above.

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50), ylim = c(0, 550))
N.current <- N.1997
for(t in 1:50){ # for loop will loop through the simulated 50 years
  
  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
}

Assignment 1 PVA

Jordan Walbridge

9/22/2020

PVA Assignment 1: for 9/15/2020 and 9/17/2020

Data: Figure 5 (plus 10 years)

census

year t

Population size

Population growth rate: 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

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

PVA for 9/17/2020

Data we will be working with for 9/17/2020 Lecture

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

Let’s Pick a Lambda out of the Hat!

Head and Tails with data frames in R

Predicting Growth Rate in 1998, by using a random lambda.

Predicting future data the ha(R)d way

Implementing simulated population of bears from year 1997 - 2005

Let’s make some plot limits