Title image Fall 2018

Toadstools and Faeries

The purpose of this project is to bring together all of the pieces--vectors, matrices, for loops, if statements, and functions--into a single ecological modeling application.

Deep Background: each year toadstools in the forest produce faeries in the spring. In the fall, each faerie plants two toadstool seeds before dissolving in the final rays of the sun as it sets on the first day of winter. If there is a hard winter (because of climate change) then fewer of the toadstool seeds germinate in the spring. If it is a milder winter, more of the toadstools survive.

The purpose of this project is to model the 200-year survival rate of toadstools given different probabilities of a hard winter. In order to do this, we will create the following functions.

The final step will be to use the image function to view the resulting data. By adjusting some parameters of the toadstool-faerie model, we can see the effects of changing different parameters.


Project Tasks

    Open R Studio or a text editor and Create a new R Script or R Markdown file (your choice) called project7. Put your name at the top of the file and save it to an appropriate working directory.

  1. Model variables: place the following variable declarations at the top of your code file. These are the parameters for the model. By modifying their values, you can test the effects of different parameters on the model.
    # model variables
    gbl_initial_pop_toadstools <- 500 # initial population of toadstools
    gbl_chance_faerie <- 0.8          # the chance a toadstool will create a faerie in the spring
    gbl_num_seeds_per_faerie <- 2     # the number of seeds a faerie can carry
    gbl_chance_seed <- 0.7            # the chance a seed will germinate
    gbl_hard_winter <- 0.7            # the percent of seeds that survive a hard winter
    gbl_mild_winter <- 0.9            # the percent of seeds that survive a mild winter
    gbl_max_toadstools <- 1000        # the max number of toadstools (carrying capacity)
    gbl_min_toadstools <- 10          # the minimum number of toadstools for a viable population
    
    # simulation variables
    gbl_NSims <- 100                  # number of simulations to execute
    gbl_NYears <- 200                 # number of years to simulation
    gbl_chance_of_hard_winter = c( 1/10.0, 1/9.0, 1/8.0, 1/7.0, 1/6.0, 1/5.0, 1/4.0, 1/3.0 )
  2. Create the function simYear, which should have two parameters: num_toadstools, hard_winter_chance. It should calculate the number of toadstools remaining after the winter and return that value. The template for the function is below.
    # calculates the number of toadstools after simulating a year
    # returns the population size at the end of the year
    simYear <- function( num_toadstools, hard_winter_chance ) {
      
      # calculate the number of faeries created
      # code: assign to num_faeries the product of num_toadstools and gbl_chance_faerie
      
      # calculate the number of seeds that are planted
      # code: assign to num_seeds the product of the num_faeries and the gbl_num_seeds_per_faerie
    
      # figure out if it was a hard or mild winter and how many toadstools survive and germinate
      # code: assign to this_year a single value from runif
      # code: if this_year is less than hard_winter_chance
          # code: assign to spring_toadstools the product of num_seeds, gbl_hard_winter, and gbl_chance_seed
      # code: else
          # code: assign to spring_toadstools the product of num_seeds, gbl_mild_winter, and gbl_chance_seed
      
      # bounds check the number of toadstools
      # code: if spring_toadstools is greater than gbl_max_toadstools
          # code: assign to spring_toadstools the value in gbl_max_toadstools
    
      # code: if spring_toadstools is less than gbl_min_toadstools
          # code: assign to spring_toadstools the value 0  
      
      # code: return spring_toadstools
    }

    Once you have written the above function, you can test it with the following code.

    print( simYear( 500, 1.0 ) )
    print( simYear( 500, 0.0 ) )

    It should print out the values 392 and 504.

  3. Create the function runSim, which should have two parameters: num_years, chance_hard_winter. The function should use the simYear function to simulate the specified number of years (num_years) and store each year's population size into a vector. The function should return a vector with num_years population values. The template is below.
    # runs a simulation for the specified number of years and chance of a hard winter
    # returns a vector with the population size for each year of the simulation
    runSim <- function(num_years, chance_hard_winter) {
      # code: assign to year_results a new double vector of size num_years 
      # code: assign pop_toadstools the value in gbl_initial_pop_toadstools
      # code: for each year of the simulation (1:num_years)
          # code: assign to pop_toadstools the result of calling simYear with the 
                  arguments pop_toadstools and chance_hard_winter
          # code: assign to the next position in the year_results vector the value 
                  pop_toadstools (use your loop variable as the index)
    
      # code: return year_results
    }

    Once you have written the above function, you can test it with the following code.

    print( round(runSim( 10, 0.0 ),0) )
    print( round(runSim( 10, 1.0 ),0) )
    set.seed(42)
    print( round(runSim( 10, 0.2 ),0) )

    You should get the following results

    [1] 504 508 512 516 520 524 529 533 537 541
    [1] 392 307 241 189 148 116  91  71  56  44
    [1] 504 508 512 516 520 524 529 414 418 421
  4. Create a function runManySims, which should have three parameters: num_sims, num_years, chance. The function should first execute num_sims simulations using the runSim function, storing the results in a matrix (one row for each simulation). It should then loop over each year and calculate the percentage of simulations with a population greater than 0, storing the result in a vector. The function should then return the vector of percentages.
    # runs a group of simulations given number of simulations, number of years, and chance of a hard winter
    # returns a single vector with the probability of survival for each year
    runManySims <- function( num_sims, num_years, chance ) {
      # calculate the matrix of num_sims v. num_years
      # code: assign to results a new matrix with num_years columns and num_sims rows
      # code: for i in the number of simulations (1:num_sims)
          code: assign to the ith row the result of calling runSim, passing in num_years and chance as the arguments
      
      # for each year, calculate the percent that is not zero
      # code: assign to prob_of_survival a new double vector of length num_years
      # code: for i in the number of years (1:num_years)
          # code: assign to the ith element of prob_of_survival the sum of the results 
                  in column that are greater than zero  which is then divided by num_sims
                  You can use the expression results[,i] > 0 to get a logical vector whose sum
                  is the number of values in column i greater than 0.
    
      # code: return prob_of_survival
    }	

    Once you have written the above function, you can test it with the following code.

    set.seed(42)
    print( round( runManySims( 100, 30, 0.5 ), 2 ) )

    You should get the following results

    [1] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
    [20] 0.99 0.99 0.99 0.99 0.98 0.96 0.94 0.87 0.81 0.73 0.64
  5. Create a function calcProbSurvival, which should have three parameters: num_sims, num_years, chance_values. The function should execute the runManySims function for all of the values in the chance_values vector (chances of a hard winter), storing the results in a matrix. The function should return the matrix.
    # runs a group of simulations, evaluating survival probabilities for each year
    # for a set of chance of hard winter values. 
    # Returns a matrix with columns being chance of hard winter and rows being year
    calcProbSurvival <- function( num_sims, num_years, chance_values ) {
      # code: assign to prob_survival a matrix with num_years rows and as many columns as the length of chance_values
      # code: assign to cur_col the value 1
      # code: for each chance in chance_values
          # code: assign to colum cur_col of the prob_survival matrix the result of 
                  executing runManySims with the arguments num_sims, num_years, and chance
          # code: assign to cur_col the value of cur_col plus 1
    
      # code: return prob_survival
    }

    Once you have written the above function, run it using the following code.

    all_results <- calcProbSurvival( gbl_NSims, gbl_NYears, gbl_chance_of_hard_winter )
  6. Use the image function to plot the matrix all_results. Use 1:gbl_NYears as the x argument, use gbl_chance_of_hard_winter as the y argument, and use all_results as the z argument. Add appropriate x and y labels and a title.

Report

Answer the following questions. Submit your answers as a plain text file or PDF in your handin directory. Put your name at the top of the file. No credit will be given for any other format or for a file without a name.


Handin

Create a project7 folder inside the Private folder in your Courses directory. Put your R Script file, and your text/PDF file with your answers into the project7 folder. If you used R Markdown for the project, put your answers and the code all in the same file.