Goals:

  1. To attempt to fit a competition model to actual data, namely Gause’s famous two-paremecium system.
  2. To practice loading data, using functions, and performing numerical experimentation in R
  3. There are 5 exercises in this lab to submit by midnight next Monday, April 1.

Note - you can copy-paste code from this lab into your R-studio, but do so carefully and place it into a script, which you can save later as an .R file.

Figure 1 (Gause 1934)

Step 1. Load data

There are two datasets posted on Blackboard or here and here: single.csv and mixture.csv.

Load these into R, using using the FILE/IMPORT DATADET point and click or … better … using the read.csv() command as follows:

single <- read.csv("single.csv")
mixture <- read.csv("mixture.csv")

Q1. What do these datasets represent?

Step 2. Plot these data

There are multiple ways to do this. For example, both of the following do the same thing:

plot(single$Day, single$caudatum)
plot(caudatum~Day, data = single)

You can make plots prettier in different ways. Compare the following:

plot(caudatum~Day, data = single, type = "l")
plot(caudatum~Day, data = single, type = "o")
plot(caudatum~Day, data = single, type = "o", col = "blue")

The plot() function makes a new plot, but you can also add points to an existing plot with the points() or lines() function. For example, see the following code:

plot(aurelia~Day, type="o", data = single, col = "darkorange")
lines(caudatum~Day, data = single, type = "o", col = "darkblue")

Q2. Use R to load the mixture.R dataset and plot those two curves - in blue and orange - on the same figure as well.

Step 3. Load some functions

Now - download and open the fittingfunctions.R file, from Blackboard or here. Run the entire code in that file, either selecting all (Ctrl-A) and running (Ctrl-Enter), clicking on the Source button under Run up top or (simple yet fancy), run the source command as follows:

source("fittingfunctions.R")

There are a bunch of functions in here that allow you to experiment with logistic and competition models. For example, to fit the logistic curve to one of the paramecium data sets, just do:

fitLogistic(y = "aurelia", time = "Day", data = single)
## 
## Formula: Y ~ N.logistic(X, N0, K, r0)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## N0   1.2780     0.7533   1.697    0.112    
## K  222.6742     5.0593  44.013  < 2e-16 ***
## r0   0.7123     0.0828   8.602 5.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.65 on 14 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 3.322e-06
##   (6 observations deleted due to missingness)

If you save the output of this function (as an R object), you can draw the curve on your figure as well using the linesLogistic() function:

au.fit <- fitLogistic(y = "aurelia", time = "Day", data = single)
plot(aurelia~Day, type="o", data = single)
linesLogistic(au.fit)

Q3: Make the fit line fat and red by adding the following to your linesLogistic function: linesLogistic(au.fit, lwd = 2, col = "red")

Q4: Perform the same fit for P. caudatum and report the estimates.

Step 4. Experiment with the Competition Simulator

The simulateCompetition() function numerically generates the 2 species process we discussed in class, i.e.:

\[{dN_1 \over dt} = r_1 N_1 \left(1-{N_1 \over K_1} - \alpha {N_2\over K_1}\right)\] \[{dN_2 \over dt} = r_2 N_2 \left(1-{N_2 \over K_2} - \beta { N_1\over K_2}\right)\]

Where everything is familiar, except there is also an \(\alpha\) which cuts \(N_1\) down based on how many \(N_2\) there are, and a \(\beta\) which cuts \(N_2\) down based on how many \(N_1\) there are.

The way it works is as follows. You first enter a bunch of parameter values in, for example:

pars <- list(K1 = 100, K2 = 100, 
             r1 = 1, r2 = 1,
             alpha = 2, beta = 1,
             n1 = 20, n2 = 10)

Then, you can pass these parameters all at once into a simulation function, and save the output into an object, here called mysim:

mysim <- simulateCompetition(tmax = 24, tstart = 2, pars)

Look at the structure of mysim - by clicking on it, or typing View(mysim) or trying either of these commands:

head(mysim)
##    Day       N1       N2
## 1 2.00 20.00000 10.00000
## 2 2.01 20.12000 10.07000
## 3 2.02 20.24020 10.14030
## 4 2.03 20.36058 10.21089
## 5 2.04 20.48115 10.28179
## 6 2.05 20.60190 10.35298
str(mysim)
## 'data.frame':    2201 obs. of  3 variables:
##  $ Day: num  2 2.01 2.02 2.03 2.04 2.05 2.06 2.07 2.08 2.09 ...
##  $ N1 : num  20 20.1 20.2 20.4 20.5 ...
##  $ N2 : num  10 10.1 10.1 10.2 10.3 ...

You can plot this, just as above. But I helpfully provided a function that does this for you:

simPlot(mysim)

Furthermore, you can plot this on top of a data set, which is nice. For example:

simPlot(mysim, data = mixture, species1 = "aurelia", species2 = "caudatum")

Obviously, this is not a very good fit! I also wrote a function which quantifies how lousy of a fit this is using \(R^2\). The closer \(R^2\) is to 1, the better. That function works like this:

getR2(mysim, data = mixture, species1 = "aurelia", species2 = "caudatum")
##        N1        N2 
## -4.767976 -1.190730

Yikes! Those are really bad numbers. They’re so far from 1, they’re negative.

Q5 or The Challenge: Experiment with this system and see if you can find values of \(\alpha\) (effect of caudatum on aurelia) and \(\beta\) (effect of aurelia on caudatum) which more or less captures the competition process observed in the data. How high can you get the \(R^2\) value to go?

Note - to do the above you basically just have to copy paste the following lines of code, adjust the numbers in the pars object, and keep experimenting. But use the following values as fixed: \(K_1 = 210\), \(K_2 = 235\), \(N_1 = N_2 = 5\). That leaves \(\alpha\), \(\beta\), \(r_1\) and \(r_2\) to play with:

pars <- list(K1 = 222, K2 = 202, 
             r1 = 0.71, r2 = 0.96,
             alpha = 0, beta = 0,
             n1 = 1.28, n2 = 0.2)

mysim <- simulateCompetition(tmax = 24, tstart = 2, pars)
simPlot(mysim, data = mixture, species1 = "aurelia", species2 = "caudatum")
getR2(mysim, data = mixture, species1 = "aurelia", species2 = "caudatum")

How close can you get yout \(R^2\) values to 1?

Bonus on logistic

These functions, by the way, work great on the logistic data as well, where there is NO competition:

pars <- list(K1 = 222, K2 = 202, 
             r1 = 0.71, r2 = 0.96,
             alpha = 0, beta = 0,
             n1 = 1.28, n2 = 0.2)

mysim <- simulateCompetition(tmax = 24, tstart = 0, pars, dt = 0.001)
simPlot(mysim, data = single, species1 = "aurelia", species2 = "caudatum")

getR2(mysim, data = single, species1 = "aurelia", species2 = "caudatum")
##        N1        N2 
## 0.9815410 0.9806087