Functions

bulleted lists in markdown

mean()
var()

population of Grizzly bears doomed bc downward trend in bears

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 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…

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 can be expressed as a ratio of population sizes.

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

lambda.59_60 <- females.N.1960/females.N.1959
(47-44)/44

## [1] 0.06818182

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

44*1.068182

## [1] 47.00001

females.N[2] #females.N.1960

## [1] 47

females.N[1] #females.N.1959

## [1] 44

females.N.1959 <- 44
females.N.1960 <- 47
females.N.1959 <- females.N[2]
females.N.1960 <- females.N[1]

females.N.1960/females.N.1959

## [1] 0.9361702

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

39 years of data but only 38 lambdas because each requires a transition from a year (ex. 1959 to 1960)

Population Growth rate: vectorized

What does this do?

females.N[2:3]

## [1] 47 46

females.N[c(2,3)]

## [1] 47 46

#gives you second and third years of data then next is first and second

females.N[1:2]

## [1] 44 47

females.N[2:3]/females.N[1:2]

## [1] 1.0681818 0.9787234

Calculate all lambdas

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

#every single lambda/transition from that graph in the paper
#we're dividing the first year by the second year and we can't divide by 0 year so start (compares 2 vs 1, 3 vs 2, etc.)

What does this do?

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

2:39 is the same as seq(2,39)

2:length(females.N)

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

2:39

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

seq(2,39)

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

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]
females.Ntemp <- females.N[seq(1,10)]
females.Ntemp <- females.N[c(1,2,3,4,5,6,7,8,9,10)]

females.Ntemp[-1]

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

# -1 removes the first value

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

Negative indexing allows you to drop a specific element from a vector

The first element

females.Ntemp[]

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

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

If you want to drop multiple elements, do the following

females.Ntemp[-c(1:2)] #vector from 1 to 2 then add the negative so it drops

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

females.Ntemp[-c(1:2)]

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

females.Ntemp[-c(1,2)]

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

females.Ntemp[-1:-2] #looks weird but also works

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

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

females.Ntemp[-1]

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

females.Ntemp[-10] #drops 10th

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

#how do we write generalized code that drops the last element of a vector?
#see below

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

females.Ntemp[-10] #drops the last element bc only 10 elements

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

#or can say: determine the length then length will correspond to last element of the vector then negative indexing will drop it
# (Drop more than one use the colon; seen previously)
females.Ntemp[-length(females.Ntemp)]

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

Calculate the first 9 lambdas.

#have ten numbers here but want to calculate 9 lambdas. drop first value then compare with second value then drop last value

#44 47 46 44 46 45 46 40 39 39

#44 47 46 44 46 45 46 40 39 39 
#   44 47 46 44 46 45 46 40 39 39 #values on end can't use for calculations

#47 46 44 46 45 46 40 39 39 numerator
#44 47 46 44 46 45 46 40 39  denominator

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

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

Calculate lambdas for all data

** TASK **

Below each line describe what the parts of the code do. Run the code to test it. * What does females.N[-1] do? * What does females.N[-length(females.N)]?

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

Create special columns

** TASK **

What does this code do?

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

** TASK **

WHat type of log does log() calculate?

lambda_log <- log(lambda.i)

R uses log as natural log; to use log 10 you need to use log10

Assemble dataframe

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

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

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

hist(bear_N$lambda.i)

Plot ln(lambda)

Challene Task * make a histogram of natural log population growth rates * write a brief sentence describing the shape of these data.

hist(bear_N$lambda_log)

Mean of log(lambda)

CHALLENGE TASK Briefly describe what happens whne you delete na.rm=T

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

## [1] 0.02134027

Histogram with Mean

In statistics the mean is often represented with mu. This can be represented as “u”. Challenge Task Save the mean to an object called u

u <- mean(bear_N$lambda_log, na.rm = T)
u

## [1] 0.02134027

Challenge Task Make a historam with the mean plotted on it

hist(bear_N$lambda_log)
abline(v = u)

Density dependence

Challenge task Make a graph that indicates if “density dependence” is occurring.

plot(lambda.i ~ females.N, data = bear_N, 
     type = "p",
     ylab = "Population growth Rate",
     xlab = "Females with cubs")
abline(v = 1970)

Thursday; Bear Data Continued

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

use computational tools to answer that question
simulate population (invent population and predict into the future what’ll happen with that population)
take into account that things vary
simulate population dynamics to learn how variable/stochastic things are
take key parameters and characterize how variable that things is
Essentially, to model population dynamics, we randomly pull population growth rates out of a hat

hist(bear_N$lambda.i)

#take numbers put them in a hat and draw them out randomly, draw out the lambda values next

hat_of_lambdas <- bear_N$lambda.i
#NA shows up and it's hard to deal with; missing data; # NA and NAN = null; Null removes an attribute

is.na() tells you T/F if something is NA

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

#gives you logical statements, not NA

Are there any values that are true?

any(is.na(hat_of_lambdas) == TRUE)

## [1] TRUE

length(hat_of_lambdas)

## [1] 39

hat_of_lambdas[39]

## [1] NA

hat_of_lambdas[-39] #dropped NA value

##  [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)] #KNOW THIS; general coding

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

omit the NA/null values

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"

hist(hat_of_lambdas)

Random Sampling of Lambdas

Every time you pull a value from the hat, you’ll end up with a different value! The values pulled will determine how much a population is growing

#Pulled a lambda from the hat
sample(x=hat_of_lambdas, size = 1, replace = TRUE)

## [1] 1.212766

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

#pull lambda from hat and save to object
#start with a new hat every time so that values are replaced after each pull
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1, replace = TRUE)
lambda_rand.t

## [1] 1.068182

Get initial population size

**Good-to-know codes to obtain general information about the dataset

head(bear_N) #show top 6 values

##   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) #bottom 6 values

##    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 #start off with 99 bears

summary(bear_N)

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

## [1] 39  5

N.1997 <- 99

One round of population simulate

Add 1-2 sentences notes here

1.22807*99  #1.22807 was the random value we pulled from the hat and that's the popln growth rate so if we want to find out how many bears in 1998 so multiply by 99 bears

## [1] 121.5789

lambda_rand.t*N.1997

## [1] 105.75

N.1998 <- lambda_rand.t*N.1997 #population size for next year

Simulation the hard way

Almost like Excel, only worse

Pull out a random lambda and multiply it by current population size

#hypothetical scenario for 50 years in the future
#1997 to 1998
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1, replace = TRUE)
N.1998 <- lambda_rand.t*N.1997

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

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

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

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

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

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

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

#carrying capacity prevents populations from getting too big

Plot Population Change

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

Start with initial population size

# Initial Conditions
N.1997 <- 99
N.initial <- 99

# explore xlim = argument (xlim sets x axis)
plot(N.1997 ~ c(1997)) #plots a single point

plot(N.1997 ~ c(1997), xlim = c(1997, 1997+50)) #orient this point all the way to the left

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

# Not a for () loop for loop
# for loop the hard way 

N.current <- N.initial

#This is where the for () loop would be
t <- 4

#randomly select a lambda
lambda_rand.t <- sample(x = hat_of_lambdas, size = 1, replace = TRUE)

#determine pop size
N.t <- N.current*lambda_rand.t

#tell R what year we're at
year.t <- 1997+t

#plot the points, where my population went
#points() updates an existing graph
points(N.t - year.t)

#update N.current
N.current <- N.t
#Change the year on line 484 from t <- 1 to t <- 2

#make a new plot with starting population size
plot(N.1997 ~c(1997), xlim = c(1997, 1997+50), ylim = c(0,550))

#start at 1997's pop size
N.current <- N.1997

#instead of having to add 1 to t each time and adding individual increments, we use this for() loop 
for (t in 1:50){
  #add comment here
  lambda_rand.t <- sample(x= hat_of_lambdas,
                          size = 1,
                          replace= TRUE)
  #Add comment here
  N.t <- N.current*lambda_rand.t
  
  #Add comment here
  year.t <- 1997+t
  
  #Add comment here
  points(N.t ~ year.t)
  
  #Add comment here
  N.current <- N.t
}

Dont read into things that aren’t there. Might think there is a carring capacity at 400.

PATTERN DEPENDENCE DOESNT MEAN A BIOLOGICAL PATTERN

Goofy R plotting code/magic # 3 x 3 panel of graphs (9 graphs to see variability)

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

Simulation: sometimes population will increase, increase and die off, stay fairly stable

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
}

10x10 plot

par(mfrow = c(10,10), mar = c(0,0,0,0), xaxt = "n", yaxt = "n")

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
}

num.years <- 50

N.vector <- rep(NA, num.years)

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

Markalanda In Class Code started 9/15/2020 Lecture_Bear Data

Shaniya Markalanda

9/17/2020

Functions

bulleted lists in markdown

population of Grizzly bears doomed bc downward trend in bears

Data: Figure 5 (plus 10 years)

Population size

Population growth rate: example

Population Growth rate: vectorized

Calculate lambdas for all data

Finish putting together dataframe

Create special columns

Assemble dataframe

Plot ln(lambda)

Mean of log(lambda)

Histogram with Mean

Density dependence

Thursday; Bear Data Continued

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

Random Sampling of Lambdas

Get initial population size

One round of population simulate

Simulation the hard way

Plot Population Change

For Loop()

Dont read into things that aren’t there. Might think there is a carring capacity at 400.

PATTERN DEPENDENCE DOESNT MEAN A BIOLOGICAL PATTERN

10x10 plot