Embedded Files
Ecological simulation can take many forms and we will not get even close to touching on them all. Simulation can be a powerful tool for understanding ecological dynamics and interactions, but it can also be an asset before you even go out and collect data!
The script for class can be downloaded here: Week 10 R script
Some quick review and useful tools
Some quick review and useful tools
For loops redux
We learned about for loops back in week 5, remember that? Well for loops are a commonly used tool in simulation. You will see them in various simulation codes. But it doesn't stop there. You might also see a for loop in a for loop... What the what??? Basically, the outer for loop runs until it hits the inner for loop and repeats. Lets work though this.
for(i in 1:3){ cat("outer loop ", i, "\n") for (j in 1:5) { cat("inner loop ", j, "\n") } }
Here is another example.
mns <- c() for (i in 10:30) { sims <- c() for (j in 1:100) { x <- sum(rnorm(i, 80, 10)) sims <- c(sims, x) } mns <- c(mns, mean(sims)) } mns
Functions
We could also think about using a function to perform the same task.
sims <- function(i) { x <- c() for (j in 1:100) { x <- sum(rnorm(i, 80, 10)) sims <- c(sims, x) } return(x) }
We can feed this function into a for loop. Using a combination of for loops and functions can really help you keep your sanity if you are embarking on a complex simulation.
mns <- c() for (i in 10:30) { mns <- c(mns, mean(sims(i))) } mns
Some basic simulation
Some basic simulation
Lets go back to Week 8 (ohhh so long ago!) and think about the model of lamprey density. Doing our literature review we have found there is a hypothesized association between depth and velocity. It would be prudent to do a bit of simulation to see how many sites we might need to visit to reliably estimate the strength of this association.
First we will simulate our covariates, depth and velocity.
n <- 10 depth <- runif(n, 0.22, 2.2) velocity <- runif(n, 0.03, 68) B0 <- 0.237 # THESE VALUES FROM A BIT OF PILOT DATA B1 <- -0.0063 # THESE VALUES FROM A BIT OF PILOT DATA B2 <- 0.125856 # THESE VALUES FROM A BIT OF PILOT DATA sim_density <- B0 + B1 * depth + B2 * velocity + rnorm(n, mean = 0, sd = 1)
Another way to add error is as follows (since this approach is how you can simulate other distribution such as the Poisson or binomial).
pred <- B0 + B1 * depth + B2 * velocity + rnorm(n, mean = 0, sd = 1) sim_density <- rnorm(n, mean = pred, sd = 1)
Lets look at what we created!
par(mfrow = c(2, 1), mar = c(4, 3, 0, 0), oma = c(1, 2, 1, 1)) plot(depth, sim_density, ylab = "", las = 1, xlab = "Depth (m)") plot(velocity, sim_density, ylab = "", las = 1, xlab = "Velocity (m/s)") mylab <- expression(paste("Density (lamprey meter "^-2, ")", sep = "")) mtext(side = 2, mylab, outer = TRUE)
sim_dat <- data.frame(depth = depth, velocity = velocity, density = sim_density) fit <- lm(density ~ depth + velocity, sim_dat) summary(fit)
We can extract the standard errors and calculate the coefficient of variation as you learned in week 8.
stderrs <- sqrt(diag(vcov(fit))) ests <- coef(fit) abs(stderrs/ests) * 100
Lets see what happens if we up n to 100
n <- 100 depth <- runif(n, 0.22, 2.2) velocity <- runif(n, 0.03, 68) sim_density <- B0 + B1 * depth + B2 * velocity + rnorm(n, mean = 0, sd = 1) sim_dat <- data.frame(depth = depth, velocity = velocity, density = sim_density) fit <- lm(density ~ depth + velocity, sim_dat) summary(fit) stderrs <- sqrt(diag(vcov(fit))) ests <- coef(fit) abs(stderrs/ests) * 100
Lets see what happens if we up n to 200
n <- 200 depth <- runif(n, 0.22, 2.2) velocity <- runif(n, 0.03, 68) sim_density <- B0 + B1 * depth + B2 * velocity + rnorm(n, mean = 0, sd = 1) sim_dat <- data.frame(depth = depth, velocity = velocity, density = sim_density) fit <- lm(density ~ depth + velocity, sim_dat) summary(fit)
stderrs <- sqrt(diag(vcov(fit)))
ests <- coef(fit) abs(stderrs/ests) * 100
We can take some of what we just did and roll it into a function that can be used in a bigger simulation.
sim <- function(n = 20) { depth <- runif(n, 0.22, 2.2) velocity <- runif(n, 0.03, 68) B0 <- 0.237 B1 <- -0.0063 B2 <- 0.125856 sim_density <- B0 + B1 * depth + B2 * velocity + rnorm(n, mean = 0, sd = 1) fit <- lm(sim_density ~ depth + velocity) betas <- coef(fit) return(betas) } sim(n = 30)
Here we take that function and use a nested for loop to first set the sample size and then run 100 replicates for that sample size and collect the parameter estimates.
sampleSize <- seq(10, 100, by = 10) output <- data.frame() for (i in 1:length(sampleSize)) { for (j in 1:100) { betas <- sim(n = sampleSize[i]) output_app <- data.frame(n = sampleSize[i], B0 = betas[1], B1 = betas[2], B2 = betas[3]) output <- rbind(output, output_app) } } head(output)
We can plot those results and see what happened in our simulation. You remember how to plot right?
par(mfrow = c(3, 1), mar = c(3, 1, 0, 1), oma = c(1, 3, 1, 1)) boxplot(B0 ~ n, output) abline(h = 0.237) boxplot(B1 ~ n, output) abline(h = -0.0063) boxplot(B2 ~ n, output) abline(h = 0.125856) mtext(side = 2, "Estimate", outer = TRUE, line = 1.5)
Here is an example for non-normal data. In this case we are simulating some counts assuming a Poisson distribution. The key here is how the Poisson distributed data is generated from the predictions using rpois()
#sample size n <- 1000 #regression coefficients beta0 <- 1 beta1 <- 0.2 beta2 <- -0.5 #generate covariate values x <- runif(n=n, min=0, max=1.5) x2<- runif(n=n, min=-3,max=3) #compute mu's mu <- exp(beta0 + beta1 * x + beta2*x2) #generate Y-values y <- rpois(n=n, lambda=mu) #data set data <- data.frame(y=y, x=x) head(data) glm(y~x + x2,data, family="poisson")
An exponential population model
An exponential population model
The exponential model requires a single parameter, the intrinsic growth rate r. In addition, you need to give the population a place to start, the initial population size. Lets make this a bit more concrete with an example. An invasive species was recently introduced on Marys Peak. The Forest Service has estimated the abundance to be 2 plants per hectare. The intrinsic growth rate of 0.1 year-1 was estimated during a nearby invasion. Assuming that this invasion will be similar, the Forest Service would like to know how large the population will be after 10 years. Just to be clear the model parameters are as follows:
r = 0.1 year-1
N0 = 2 plants hectare-1
Here is how we would go about doing this simulation using R. First we need to assign the objects that represent r, N0, and the number of years to simulate as follows.
N0 <- 2 r <- 0.1 nyears <- 10
Next we are going to create an object to hold the simulated population densities. This is easily done by making a vector of NAs that is 11 elements long. Why 11? Well this is to allow us to set the initial population as the first element and then have 10 more elements to hold the simulated population.
N <- rep(NA, nyears + 1) N <- N0 for (i in 1:nyears) { N[i + 1] <- r * N[i] + N[i] } N
And we can plot the simulation output…
plot(c(0:10), N, ylab = "Population size (plants/hectare)", xlab = "Year", las = 1, type = "b")
There are several different ways to approach this problem in R, all producing the same results. Here are a couple ways you may see code specified if you happen to be Googling around.
The approach below simulates from year 2 to 11 and simulates the population given the previous population size (i.e, Nt-1).
N <- rep(NA, nyears + 1) N[1] <- N0 for (i in 2:(nyears + 1)) { N[i] <- r * N[i - 1] + N[i - 1] } N
The approach below uses the collection operator instead of creating a vector. One thing to consider is that this method can be very efficient in terms of computing time.
N <- N0 for (i in 1:nyears) { Nt1 <- r * N[i] + N[i] N <- c(N, Nt1) } N
Logistic model
Logistic model
The exponential model is not a good approximation of population dynamics beyond the initial rapid colonization. The addition of a carrying capacity effectively allows the population to approach a maximum size. The equation is similar to the exponential model with the addition of an additional parameter K.
Building on the previous example, previous studies indicate that populations in similar areas never exceed 40 plants hectare-1. Just to be clear the model parameters are as follows:
r = 0.1 year-1
N0 = 2 plants hectare-1
K = 40
N0 <- 2 r <- 0.1 nyears <- 50 K <- 40 N <- rep(NA, nyears + 1) N[1] <- N0 for (i in 1:nyears) { N[i + 1] <- r * N[i] * ((K - N[i])/K) + N[i] } N
Lets look at what we have going on here with this logistic model.
plot(c(1:51),N,ylab="Abundance",xlab="years")
Breaking the for loop
Sometimes in simulations you need to know when a condition occurs. Continuing with the invasion example, the Forest Service will apply an herbicide treatment when the population exceeds 30 plants per acre, but they need to know when this will happen given the model (If you want to know how to make decisions like these in a structured way check out Jim's 544 class or his book here).
A for loop can be stopped using the break argument.
N0 <- 2 r <- 0.1 nyears <- 50 K <- 40 N <- rep(NA, nyears + 1) N[1] <- N0 for (i in 1:nyears) { N[i + 1] <- r * N[i] * ((K - N[i])/K) + N[i] if (N[i + 1] >= 30) break } N
You can figure out how many years it took to exceed 30 plants hectare-1 using the length function omiting the NAs (of course you remember na.omit() right?).
length(na.omit(N))
Age/stage structured population model
Fecundity: Age 1 = 0, Age 2 = 2, Age 3 = 3, Age 4 = 4
Survival: Age 1 to 2: 0.6, Age 2 to 3: 0.3, Age 3 to 4: 0.24
# AN AGE OR STAGE MATRIX OF FECUNDITIES AND SURVIVALS A<- matrix(c(0,2,3,4, 0.6,0,0,0, 0,0.3,0,0, 0,0,0.24,0),nrow=4, byrow=T) # A MATRIX TO HOLD THE FORECASTS OF THE POP population<- matrix(0, nrow=4, ncol=10) # SEED THE POPULATION population[,1]<- c(400, 200, 100, 40) # FORECAST THE STRUCTURED POPULATION for(i in 2:10){ population[,i]<- A%*%population[,i-1] } # CALCULATE THE TOTAL POPULATION AT EACH TIME STEP totalPopulation<- apply(population, 2, sum) plot(totalPopulation)
plot(c(1:10),population[1,],ylim=c(0,8000),type='l') points(c(1:10),population[2,],type='l') points(c(1:10),population[3,],type='l') points(c(1:10),population[4,],type='l')
# LAMBDA FOR MATRIX A lambda<-as.numeric(eigen(A)$values[1]) w<-as.numeric(eigen(A)$vectors[,1]) w<-w/sum(w) # STABLE AGE DISTRIBUTION FOR MATRIX A w v<-as.numeric(eigen(t(A))$vectors[,1]) v<-v/v[1]; # REPRODUCTIVE VALUE FOR MATRIX A v s<-outer(v,w)/sum(v*w) # SENSITIVITY FOR MATRIX A s e<-(s*A)/lambda # ELASTICITIES FOR MATRIX A e
Exercises
Exercises
1) Simulate a population for a 10 year period assuming the population is described by a logistic function with the following parameters:
r = 0.15 year-1
N0 = 6 plants hectare-1
K = 49
Do the following: simulate the population for 100 years and report the population size at year 89 and provide a plot of the population over time.
2) Using the same model described above do the following:
Simulate the population for a 10 year period adding normally distributed error with a mean = 0 and a standard deviation of 1.
Repeat the simulation 10 times and plot these 10 replicates on a single plot
3) Comment on anything you wish we would have covered that was not covered or you thought could deserve more depth.
Page updated
Google Sites
Report abuse