1. (2 point). The runif() function draws randomly from a uniform distribution. What would the average value be for a vector created with the following code:

Question1_Vector<-runif(10,2,4) #Hint: Think about how and if this changes if you run it multiple times?

The expected average value for a uniform distribution between 2 and 4 is 3

2. (4 point). Start with the constant growth model from the example script (found in the chunk below). Add in immigration and emmigration into the model. Immigration should range from 0 to 4. Emmigration should range from 1 to 2. Create a plot that tracks total population, total births, immigration, and emmigration. BE SURE TO INCLUDE: a legend, x axis, y axis

Hint: You will have to create a data frame as the DF_Stoich from the example script will not work, if you change the name of the data frame to anything other than DF_Stoic, you will have to update that in the code below.

Be sure you have removed “#” infront of any lines of code you want to run

#for(i in 1:50){
#  
#  b<-round(runif(1,2,5),0) #births from 2 to 5, round it
#  d<-round(runif(1,0,3),0) #deaths from 0-3
#  g<-b-d #growth is birth - deaths
#  
#  DF_Stoich[1,1]<-10 #This could be written DF_Stoich$Population[1], set initial population to 10
#  DF_Stoich[1,2]<-0 #set initial births
#  DF_Stoich[1,3]<-0 #set initial deaths
#  DF_Stoich[1,4]<-0 #set initial growth
#  
#  DF_Stoich[i+1,1]<-DF_Stoich[i,1]+g #This could be written DF_Stoich$Population[i+1],
#                                     #Population at i+1 is population at i + growth
#  DF_Stoich[i+1,2]<-DF_Stoich[i,2]+b #This tracks total births
#  DF_Stoich[i+1,3]<-DF_Stoich[i,3]+d #This tracks total deaths
#  DF_Stoich[i+1,4]<-DF_Stoich[i,4]+g #This tracks total growth (which is just the difference between)
#                                     #Initial population and grow
#} #This ends our for loop
n <- 50
DF_Stoich <- data.frame(
  Population = numeric(n + 1),
  Births = numeric(n + 1),
  Deaths = numeric(n + 1),
  Growth = numeric(n + 1),
  Immigration = numeric(n + 1),
  Emigration = numeric(n + 1)
)

DF_Stoich[1, ] <- c(10, 0, 0, 0, 0, 0)  # Initial values

for (i in 1:n) {
  b <- round(runif(1, 2, 5), 0)
  d <- round(runif(1, 0, 3), 0)
  im <- round(runif(1, 0, 4), 0)
  em <- round(runif(1, 1, 2), 0)
  g <- b - d + im - em
  
  DF_Stoich[i + 1, ] <- c(
    DF_Stoich[i, 1] + g,
    DF_Stoich[i, 2] + b,
    DF_Stoich[i, 3] + d,
    DF_Stoich[i, 4] + g,
    DF_Stoich[i, 5] + im,
    DF_Stoich[i, 6] + em
  )
}

# Basic plot
plot(DF_Stoich$Population, type = "l", col = "black", ylim = c(0, max(DF_Stoich$Population)), 
     xlab = "Time (Iterations)", ylab = "Count", main = "Population Growth with Immigration and Emigration")
lines(DF_Stoich$Births, col = "blue")
lines(DF_Stoich$Deaths, col = "red")
lines(DF_Stoich$Immigration, col = "green")
lines(DF_Stoich$Emigration, col = "purple")
legend("topright", legend = c("Population", "Births", "Deaths", "Immigration", "Emigration"), 
       col = c("black", "blue", "red", "green", "purple"), lty = 1)

3. (4 point). In the code below, an IF statement has been added into the for loop (lines 79-90). What is this IF statement doing? How could an IF statement like this be useful in modeling population changes.

Hint, plotting this data may be helpful in answering this question.

#Create Data Frame
Question_3_Data_Frame<-data.frame("Population" =c(rep(NA,50)),
                      "Births" =c(rep(NA,50)),
                      "Death" =c(rep(NA,50)),
                      "Growth" =c(rep(NA,50)))

#Run for loop 24 times, since the first value is defined by us, we only need
#to generate the next 24 values
for(i in 1:50){
  
  b<-round(runif(1,1,3),0) #births from 2 to 5, round it
  d<-round(runif(1,0,5),0) #deaths from 0-3
  g<-b-d #growth is birth - deaths
  
  Question_3_Data_Frame[1,1]<-10 #This could be written Question_3_Data_Frame$Population[1], set initial population to 10
  Question_3_Data_Frame[1,2]<-0 #set initial births
  Question_3_Data_Frame[1,3]<-0 #set initial deaths
  Question_3_Data_Frame[1,4]<-0 #set initial growth
  
  Question_3_Data_Frame[i+1,1]<-Question_3_Data_Frame[i,1]+g #This could be written Question_3_Data_Frame$Population[i+1],
  #Population at i+1 is population at i + growth
  Question_3_Data_Frame[i+1,2]<-Question_3_Data_Frame[i,2]+b #This tracks total births
  Question_3_Data_Frame[i+1,3]<-Question_3_Data_Frame[i,3]+d #This tracks total deaths
  Question_3_Data_Frame[i+1,4]<-Question_3_Data_Frame[i,4]+g #This tracks total growth (which is just the difference between) Initial population and grow
  
  if(Question_3_Data_Frame[i,1]<1){
    Question_3_Data_Frame[i,1]<-0
    Question_3_Data_Frame[i+1,1]<-0
    
    Question_3_Data_Frame[i,2]<-0
    Question_3_Data_Frame[i+1,2]<-0
    
    Question_3_Data_Frame[i,3]<-0
    Question_3_Data_Frame[i+1,3]<-0
    
    Question_3_Data_Frame[i,4]<-0
    Question_3_Data_Frame[i+1,4]<-0
  }
  
  
} #This ends our for loop

The IF statement ensures that if the population drops below 1, all future values remain at zero, preventing negative population values. This mimics real life extinction scenarios.

Q4. (15 Points). Take the following steps to incorperate natural disasters or tother major events into our models? Please include the code for all parts instead of one section of code that incorperates everything as I will give partial credit even if the final chunk of code has errors.

  • Generate two different strings representing time. We will be run this model for 3 years but our growth rates and other changes of interest will occur at the month by month level so we need to track what month it is. Using a numeric string, create a vector that tracks the month over 3 years. Hint, this should have a length of 36 (1 point)
# Numeric vector tracking months over 3 years
time_numeric <- rep(1:12, 3)
print(time_numeric)  # Should have a length of 36
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1
## [26]  2  3  4  5  6  7  8  9 10 11 12
# Character string version of the same time representation
time_character <- rep(c("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
                        "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"), 3)
print(time_character)  # Also length of 36
##  [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"
## [13] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"
## [25] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"
  • Create a data frame that will hold 100 population runs over the 3 years. Make sure the length of this data frame matches the length of your time variable (1 point)
# Define the number of months and iterations
num_months <- length(time_numeric)  # Should be 36
num_iterations <- 100  

# Create an empty data frame to store population values
pop_data <- as.data.frame(matrix(NA, nrow=num_months, ncol=num_iterations))  
print(dim(pop_data))  # Check dimensions (36, 100)
## [1]  36 100
  • Name the columns of this data frame “Model_Iteration_1” through “Model_Iteration_100” (1 point)
# Naming the columns appropriately
colnames(pop_data) <- paste0("Model_Iteration_", 1:num_iterations)
print(head(colnames(pop_data)))  # Check column names
## [1] "Model_Iteration_1" "Model_Iteration_2" "Model_Iteration_3"
## [4] "Model_Iteration_4" "Model_Iteration_5" "Model_Iteration_6"
  • Create a for loop(s) that fills tracks 100 model runs. Carrying capacity (k) should equal 50. Starting population (N0) should equal 15. Growth rate should range from -5% to 7%. (2 points)
for (j in 1:num_iterations) {
  N <- 15  # Initial population
  k <- 50  # Carrying capacity
  growth_rate <- runif(num_months, -0.05, 0.07)  # Growth rate for each month

  for (i in 1:num_months) {
    N <- min(k, N + round(N * growth_rate[i]))  # Apply growth
    pop_data[i, j] <- N  # Store population
  }
}
  • This population lives in an area where sometimes very large storms occur! Storms only happen in April and May (month 4 and 5). There is a 10% chance that in any April or May, a storm will kill 30 individuals in the population. Use if statements to add this storms into the model (5 points)
for (j in 1:num_iterations) {
  N <- 15  
  k <- 50  
  growth_rate <- runif(num_months, -0.05, 0.07)  

  for (i in 1:num_months) {
    # Check if the current month is April (4) or May (5)
    if (time_numeric[i] == 4 | time_numeric[i] == 5) {
      if (runif(1) < 0.1) {  # 10% chance of storm
        N <- max(0, N - 30)  # Ensure population does not go below zero
      }
    }
    N <- min(k, N + round(N * growth_rate[i]))  
    pop_data[i, j] <- N  
  }
}

HINT: You will need one if statement that determines whether or not a storm can occur (this if statement will determine if it is april or may). You will need a second if statement (nested in the first if statement) to determine whether the storm happened. I have included an example in the chunk below of an if statement similar to what you may need. However, this will obvious need to be altered

  • In many cases, the storm may reduce the population to 0. In this case there will be no more reproduction or population growth. Incorperate an if statement that sets population to 0 for all time periods after a crash Hint, look at the code in question 3 (2 points)
for (j in 1:num_iterations) {
  N <- 15  # Starting population
  k <- 50  # Carrying capacity
  growth_rate <- runif(num_months, -0.05, 0.07)  # Growth rate per month

  for (i in 1:num_months) {
    
    if (N > 0) {  # Only update if population is still alive
      
      # Storm check (only in April and May)
      if (time_numeric[i] == 4 | time_numeric[i] == 5) {
        if (runif(1) < 0.1) {  # 10% chance of storm happening
          N <- max(0, N - 30)  # Reduce population by 30 but not below zero
        }
      }

      # Growth only happens if N is still > 0 after storm
      if (N > 0) {
        N <- min(k, N + round(N * growth_rate[i]))  # Apply logistic growth
      }
      
    }  # End of if (N > 0) check

    pop_data[i, j] <- N  # Store updated population in data frame
  }
}
  • Plot all model runs form your model as well as the average population growth (3 points)
plot(time_numeric, pop_data[,1], type='l', ylim=c(0, max(pop_data, na.rm=TRUE)), 
     xlab='Month', ylab='Population', col='gray', main="Population Model with Storms")

for (j in 2:num_iterations) {
  lines(time_numeric, pop_data[,j], col='gray')
}

# Add the average population growth in red
lines(time_numeric, rowMeans(pop_data, na.rm=TRUE), col='red', lwd=2)  

# Add a legend for clarity
legend("topright", legend=c("Individual Runs", "Average Population"), col=c("gray", "red"), lwd=2)

-BONUS (2 points!) Write a line of code that automatically counts how many times the population crashes

# Create vector for month
time <- c(rep(1:12, 3))  # Repeat months 1-12 for 3 years, creating a vector of length 36

# Create dataframe to track changes
num_months <- length(time)  # Total number of months (should be 36)
num_iterations <- 100  # We will run 100 different model simulations

pop_data <- matrix(NA, nrow=num_months, ncol=num_iterations)  # Create an empty matrix for population data
colnames(pop_data) <- paste0("Model_Iteration_", 1:num_iterations)  # Name the columns to track each model iteration

# Counter for population crashes
crash_count <- 0  # Variable to count the number of times population crashes

# For Loops to populate dataframe
for (j in 1:num_iterations) {  # Loop over 100 model runs
  N <- 15  # Starting population size
  k <- 50  # Carrying capacity
  growth_rate <- runif(num_months, -0.05, 0.07)  # Generate random growth rate between -5% and 7%

  crashed <- FALSE  # Track if this model run has crashed

  for (i in 1:num_months) {  # Loop over months
    
    # Incorporate if statement into model
    # The example if loop tests to see if the month is April (month 4) or May (month 5)
    if (N > 0) {  # Only update population if it is still alive
      
      if (time[i] == 4 | time[i] == 5) {  # Check if the current month is April or May
        storm_chance <- runif(1)  # Draw a random number between 0 and 1
        
        if (storm_chance < 0.1) {  # 10% chance that a storm occurs
          N <- max(0, N - 30)  # If a storm happens, decrease population by 30 (but not below 0)
        }
      }
      
      # Apply growth only if the population is still alive
      if (N > 0) {
        N <- min(k, N + round(N * growth_rate[i]))  # Growth model with carrying capacity
      } else if (!crashed) {  
        # If population reaches 0 for the first time in this run
        crashed <- TRUE  # Mark this iteration as crashed
        crash_count <- crash_count + 1  # Increase total crash count
      }
      
    }  # End of if (N > 0) check

    pop_data[i, j] <- N  # Store updated population in data frame
  }
}

# The | symbol (found below backslash on your keyboard) represents the logical test OR. 
# The following statement will print "TRUE" if x equals 2 OR if x equals 5
x <- 5
if (x == 2 | x == 5) {
  print("TRUE")
}
## [1] "TRUE"
# Print the number of crashes across all model runs
print(paste("Total population crashes:", crash_count))
## [1] "Total population crashes: 38"