Lecture #18: Control Structures III

Tip #6: Use `break` to set up a conditional to break out of loop early

# create a simple random growth population model function

##################################################
# FUNCTION: ran_walk
# stochastic random walk 
# input: times = number of time steps
#        n1 = initial population size (= n[1])
#        lambda = finite rate of increase
#        noise_sd = sd of a normal distribution with mean 0
# output: vector n with population sizes > 0 
#         until extinction, then NA 
#------------------------------------------------- 
library(ggplot2)
ran_walk <- function(times=100,n1=50,lambda=1.00,noise_sd=10) {
                n <- rep(NA,times)  # create output vector
                n[1] <- n1 # initialize with starting population size
                noise <- rnorm(n=times,mean=0,sd=noise_sd) # create noise vector
                for(i in 1:(times-1)) {
                  n[i + 1] <- lambda*n[i] + noise[i]
                  if(n[i + 1] <=0) {
                    n[i + 1] <- NA
                    cat("Population extinction at time",i-1,"\n")
                    break}
                }

return(n)
}

# explore paramaters in plot function
qplot(x=1:100,y=ran_walk(),geom="line")

## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to
## see where this warning was generated.

## Population extinction at time 35

qplot(x=1:100,y=ran_walk(noise_sd=0),geom="line")

qplot(x=1:100,y=ran_walk(lambda=0.92,noise_sd=0),geom="line")

Extensions of the model for realistic populations

discrete integers to represent counts of individuals (use round())
extinction of sexually reproducing population if all same sex (use unisexExtinct == runif(1) <= 2*(0.5)^n
add environmental noise that is not a random walk (add rnorm(0,1))
add measurement error (add rnorm(0,1)) for 0s without extinction

Using double for loops

m <- matrix(round(runif(20),digits=2),nrow=5)
# loop over rows
for (i in 1:nrow(m)) { # could use for (i in seq_len(nrow(m)))
  m[i,] <- m[i,] + i
} 
print(m)

##      [,1] [,2] [,3] [,4]
## [1,] 1.14 1.02 1.49 1.38
## [2,] 2.51 2.69 2.89 2.26
## [3,] 3.93 3.61 3.12 3.69
## [4,] 4.02 4.62 4.21 4.19
## [5,] 5.44 5.09 5.24 5.75

# Loop over columns
m <- matrix(round(runif(20),digits=2),nrow=5)
for (j in 1:ncol(m)) {
  m[,j] <- m[,j] + j
}
print(m)

##      [,1] [,2] [,3] [,4]
## [1,] 1.28 2.42 3.76 4.80
## [2,] 1.78 2.02 3.48 4.64
## [3,] 1.30 2.23 3.68 4.78
## [4,] 1.96 2.50 3.32 4.36
## [5,] 1.67 2.79 3.67 4.57

# Loop over rows and columns
m <- matrix(round(runif(20),digits=2),nrow=5)
for (i in 1:nrow(m)) {
  for (j in 1:ncol(m)) {
    m[i,j] <- m[i,j] + i + j
  } # end of column j loop
} # end or row i loop
 print(m)

##      [,1] [,2] [,3] [,4]
## [1,] 2.97 3.96 4.46 5.04
## [2,] 3.57 4.42 5.70 6.38
## [3,] 4.04 5.18 6.31 7.35
## [4,] 5.99 6.46 7.82 8.39
## [5,] 6.06 7.93 8.14 9.30

Writing functions for equations and sweeping over parameters

# S = cA^z species area function, but what does it look like??

##################################################
# function: SpeciesAreaCurve
# creates power function relationship for S and A
# input: A is a vector of island areas
#        c is the intercept constant
#        z is the slope constant
# output: S is a vector of species richness values
#------------------------------------------------- 
species_area_curve <- function(A=1:5000,c= 0.5,z=0.26){
  
  S <- c*(A^z)
return(S)
}
head(species_area_curve())

## [1] 0.5000000 0.5987394 0.6653061 0.7169776
## [5] 0.7598051 0.7966899

##################################################
# function: species_area_plot
# plot species area curves with parameter values
# input: A = vector of areas
#        c = single value for c parameter
#        z = single value for z parameter

# output: smoothed curve with parameters in graph
#------------------------------------------------- 
species_area_plot <- function(A=1:5000,c= 0.5,z=0.26) {
plot(x=A,y=species_area_curve(A,c,z),type="l",xlab="Island Area",ylab="S",ylim=c(0,2500))
mtext(paste("c =", c,"  z =",z),cex=0.7) 
#      return()
}
species_area_plot()

Now build a grid of plots!

# global variables
c_pars <- c(100,150,175)
z_pars <- c(0.10, 0.16, 0.26, 0.3)
par(mfrow=c(3,4))
for (i in seq_along(c_pars)) {
  for (j in seq_along(z_pars)) {
    species_area_plot(c=c_pars[i],z=z_pars[j])
  }
}

par(mfrow=c(1,1))

Looping with `while` or `repeat`

# looping with for
cut_point <- 0.1
z <- NA
ran_data <- runif(100)
for (i in seq_along(ran_data)) {
 z <- ranData[i]
 if (z < cut_point) break
}
print(z)

# looping with while

z <- NA
cycle_number <- 0
while (is.na(z) | z >= cut_point) {
  z <- runif(1)
  cycle_number <- cycle_number + 1
  
}
print(z)
print(cycle_number)

# looping with repeat

z <- NA
cycle_number <- 0

repeat {
  z <- runif(1)
cycle_number <- cycle_number + 1
  
  if (z <= cut_point) break
}
print(z)
print(cycle_number)

# add code for cycle number
# try setting limit to 0.001

#

Using the `expand.grid()` function to create a dataframe with parameter combinations

expand.grid(c_pars,z_pars)

##    Var1 Var2
## 1   100 0.10
## 2   150 0.10
## 3   175 0.10
## 4   100 0.16
## 5   150 0.16
## 6   175 0.16
## 7   100 0.26
## 8   150 0.26
## 9   175 0.26
## 10  100 0.30
## 11  150 0.30
## 12  175 0.30

##################################################
# function: sa_output
# Summary stats for species-area power function
# input: vector of predicted species richness 
# output: list of max-min, coefficient of variation 
#------------------------------------------------- 
sa_output <- function(S=runif(1:10)) {

sum_stats <- list(s_gain=max(S)-min(S),s_cv=sd(S)/mean(S))

return(sum_stats)
}
sa_output()

## $s_gain
## [1] 0.9750248
## 
## $s_cv
## [1] 0.4131853

# Build program body with a single loop through 
# the parameters in modelFrame

# Global variables
area <- 1:5000
c_pars <- c(100,150,175)
z_pars <- c(0.10, 0.16, 0.26, 0.3)

# set up model frame
model_frame <- expand.grid(c=c_pars,z=z_pars)
model_frame$s_gain <- NA
model_frame$s_cv <- NA
print(model_frame)

##      c    z s_gain s_cv
## 1  100 0.10     NA   NA
## 2  150 0.10     NA   NA
## 3  175 0.10     NA   NA
## 4  100 0.16     NA   NA
## 5  150 0.16     NA   NA
## 6  175 0.16     NA   NA
## 7  100 0.26     NA   NA
## 8  150 0.26     NA   NA
## 9  175 0.26     NA   NA
## 10 100 0.30     NA   NA
## 11 150 0.30     NA   NA
## 12 175 0.30     NA   NA

# cycle through model calculations
for (i in 1:nrow(model_frame)) {
  
  # generate S vector
  temp1 <- species_area_curve(A=area,
                            c=model_frame[i,1],
                            z=model_frame[i,2])
  # calculate output stats
  temp2 <- sa_output(temp1)
  # pass results to columns in data frame
  model_frame[i,c(3,4)] <- temp2
  
}
print(model_frame)

##      c    z    s_gain       s_cv
## 1  100 0.10  134.3673 0.09109192
## 2  150 0.10  201.5509 0.09109192
## 3  175 0.10  235.1428 0.09109192
## 4  100 0.16  290.6926 0.13905495
## 5  150 0.16  436.0389 0.13905495
## 6  175 0.16  508.7121 0.13905495
## 7  100 0.26  815.6557 0.21070232
## 8  150 0.26 1223.4835 0.21070232
## 9  175 0.26 1427.3975 0.21070232
## 10 100 0.30 1187.3333 0.23699746
## 11 150 0.30 1780.9999 0.23699746
## 12 175 0.30 2077.8333 0.23699746

Parameter sweeping redux with `ggplot` graphics

library(ggplot2)

area <- 1:5 #keep this very small initially
c_pars <- c(100,150,175)
z_pars <- c(0.10, 0.16, 0.26, 0.3)

# set up model frame
model_frame <- expand.grid(c=c_pars,z=z_pars,A=area)
model_frame$S <- NA


# loop through the parameters and fill with SA function

for (i in 1:length(c_pars)) {
  for (j in 1:length(z_pars)) {
    model_frame[model_frame$c==c_pars[i] & model_frame$z==z_pars[j],"S"] <-   species_area_curve(A=area,c=c_pars[i],z=z_pars[j])
  }
}


for (i in 1:nrow(model_frame)) {
  model_frame[i,"S"] <- species_area_curve(A=model_frame$A[i],
                   c=model_frame$c[i],
                   z=model_frame$z[i])
}
# print(modelFrame) # check by printing a data frame with limited parameter values

library(ggplot2)

p1 <- ggplot(data=model_frame)
p1 + geom_line(mapping= aes(x=A,y=S)) +
  facet_grid(c~z)

p2 <- p1
p2 + geom_line(mapping=aes(x=A,y=S,group=z)) +
  facet_grid(.~c)

p3 <- p1
p3 + geom_line(mapping=aes(x=A,y=S,group=c)) +
  facet_grid(z~.)

Lecture #18: Control Structures III

Nicholas J. Gotelli

6 April 2021

Tip #6: Use `break` to set up a conditional to break out of loop early

Extensions of the model for realistic populations

Using double for loops

Writing functions for equations and sweeping over parameters

Now build a grid of plots!

Looping with `while` or `repeat`

Using the `expand.grid()` function to create a dataframe with parameter combinations

Parameter sweeping redux with `ggplot` graphics

Lecture #18: Control Structures III

Nicholas J. Gotelli

6 April 2021

Tip #6: Use break to set up a conditional to break out of loop early

Extensions of the model for realistic populations

Using double for loops

Writing functions for equations and sweeping over parameters

Now build a grid of plots!

Looping with while or repeat

Using the expand.grid() function to create a dataframe with parameter combinations

Parameter sweeping redux with ggplot graphics

Tip #6: Use `break` to set up a conditional to break out of loop early

Looping with `while` or `repeat`

Using the `expand.grid()` function to create a dataframe with parameter combinations

Parameter sweeping redux with `ggplot` graphics