inst/doc/gauseR_examples.R

In gauseR: Lotka-Volterra Models for Gause's 'Struggle for Existence'

## -----------------------------------------------------------------------------
# load package
require(gauseR)

# load data
data(gause_1934_book_f22)

logistic_data<-gause_1934_book_f22[gause_1934_book_f22$Treatment=="Pa",]

plot(Volume_Species2~Day, logistic_data)

## -----------------------------------------------------------------------------
# calculate time-lagged abundance
lagged_data<-get_lag(x=logistic_data$Volume_Species2, time = logistic_data$Day)

# calculate per-capita growth rate
lagged_data$dNNdt <- percap_growth(x = lagged_data$x, laggedx = lagged_data$laggedx, dt = lagged_data$dt)

# plot relationship
plot(dNNdt~laggedx, lagged_data,
     xlab="Lagged Abundance (N)",
     ylab="Per-capita Growth Rate (dN/Ndt)",
     xlim=c(0, 250), ylim=c(0, 1))
abline(h=0, v=0, lty=3)

# fit model to relationship
mod<-lm(dNNdt~laggedx, lagged_data)
abline(mod, lwd=2, col=2)

# label parameters
arrows(25, 0.6, 1, 0.8, length = 0.1, lwd=1.5)
text(25, 0.6, "y-intercept: r", pos=1)

arrows(200, 0.4, 232, 0.01, length = 0.1, lwd=1.5)
text(200, 0.4, "x-intercept: K", pos=3)

arrows(80, predict(mod, newdata=data.frame(laggedx=80)),
         130, predict(mod, newdata=data.frame(laggedx=80)),
         length = 0.1, lwd=1.5)
arrows(130, predict(mod, newdata=data.frame(laggedx=80)),
       130, predict(mod, newdata=data.frame(laggedx=130)),
       length = 0.1, lwd=1.5)
text(130, predict(mod, newdata=data.frame(laggedx=80)), "slope: aii", pos=3)

## -----------------------------------------------------------------------------
plot(Volume_Species2~Day, logistic_data, xlab="Day", ylab="P. aurelia Std. Volume")

timelst<-seq(0, 25, by=0.1) # sequence of time values to plot
prediction<-get_logistic(time = timelst, N0 = 0.6, r = 0.8, K=230)

lines(timelst, prediction, lwd=2)

## -----------------------------------------------------------------------------
prediction_short<-get_logistic(time = logistic_data$Day, N0 = 0.6, r = 0.8, K=230)

test_goodness_of_fit(observed = logistic_data$Volume_Species2, predicted = prediction_short)

## -----------------------------------------------------------------------------
# load data from competition experiment
data(gause_1934_book_f32)

# keep all data - no separate treatments exist for this experiment
predatorpreydata<-gause_1934_book_f32

# get time-lagged observations for each species
prey_lagged<-get_lag(x = predatorpreydata$Individuals_Prey, time = predatorpreydata$Day)
predator_lagged<-get_lag(x = predatorpreydata$Individuals_Predator, time = predatorpreydata$Day)

# calculate per-capita growth rates
prey_dNNdt<-percap_growth(x = prey_lagged$x, laggedx = prey_lagged$laggedx, dt = prey_lagged$dt)
predator_dNNdt<-percap_growth(x = predator_lagged$x,
      laggedx = predator_lagged$laggedx, dt = predator_lagged$dt)

# fit linear models to dNNdt, based on average
# abundances between current and lagged time steps
prey_mod_dat<-data.frame(prey_dNNdt=prey_dNNdt, prey=prey_lagged$laggedx,
      predator=predator_lagged$laggedx)
mod_prey<-lm(prey_dNNdt~prey+predator, data=prey_mod_dat)

predator_mod_dat<-data.frame(predator_dNNdt=predator_dNNdt,
      predator=predator_lagged$laggedx, prey=prey_lagged$laggedx)
mod_predator<-lm(predator_dNNdt~predator+prey, data=predator_mod_dat)

# model summaries
summary(mod_prey)
summary(mod_predator)

# extract parameters
# growth rates
r1 <- unname(coef(mod_prey)["(Intercept)"])
r2 <- unname(coef(mod_predator)["(Intercept)"])

# self-limitation
a11 <- unname(coef(mod_prey)["prey"])
a22 <- unname(coef(mod_predator)["predator"])

# effect of Pa on Pc
a12 <- unname(coef(mod_prey)["predator"])
# effect of Pc on Pa
a21 <- unname(coef(mod_predator)["prey"])

# run ODE:
# make parameter vector:
parms <- c(r1, r2, a11, a12, a21, a22)
initialN <- c(4, 0.1)
out <- deSolve::ode(y=initialN, times=seq(1, 17, length=100), func=lv_interaction, parms=parms)
matplot(out[,1], out[,-1], type="l",
        xlab="time", ylab="N", col=c("black","red"), lty=c(1,3), lwd=2, ylim=c(0, 60))
legend("topright", c("Pc", "Dn"), col=c(1,2), lwd=2, lty=c(1,3))

# now, plot in points from data
points(predatorpreydata$Day, predatorpreydata$Individuals_Predator , col=2)
points(predatorpreydata$Day, predatorpreydata$Individuals_Prey, col=1)

## -----------------------------------------------------------------------------
# Data for the optimizer:
# Must have a column with time steps labeled 'time', and
# columns for each species in the community.
opt_data<-data.frame(time=predatorpreydata$Day, Prey=predatorpreydata$Individuals_Prey,
    Predator=predatorpreydata$Individuals_Predator)

# Save the signs of the parameters -
# optimizer works in log space, so these
# must be specified separately
parm_signs<-sign(parms)

# parameter vector for optimizer -
# must be a vector with, first, the
# starting abundances in log space,
# and second, the parameter values,
# again in log space
pars<-c(log(initialN), log(abs(parms)))

# run optimizer
optout<-optim(par = pars, fn = lv_optim, hessian = TRUE,
             opt_data=opt_data, parm_signs=parm_signs)

# extract parameter vector:
parms <- exp(optout$par[-c(1:2)])*parm_signs
initialN <- exp(optout$par[1:2])

out <- deSolve::ode(y=initialN, times=seq(1, 17, length=100), func=lv_interaction, parms=parms)
matplot(out[,1], out[,-1], type="l",
        xlab="time", ylab="N", col=c("black","red"), lty=c(1,3), lwd=2, ylim=c(0, 60))
legend("topright", c("Pc", "Dn"), col=c(1,2), lwd=2, lty=c(1,3))

# now, plot in points from data
points(predatorpreydata$Day, predatorpreydata$Individuals_Predator , col=2)
points(predatorpreydata$Day, predatorpreydata$Individuals_Prey, col=1)

## -----------------------------------------------------------------------------
#load competition data
data("gause_1934_science_f02_03")

#subset out data from species grown in mixture
mixturedat<-gause_1934_science_f02_03[gause_1934_science_f02_03$Treatment=="Mixture",]

#extract time and species data
time<-mixturedat$Day
species<-data.frame(mixturedat$Volume_Species1, mixturedat$Volume_Species2)
colnames(species)<-c("P_caudatum", "P_aurelia")

#run wrapper
gause_out<-gause_wrapper(time=time, species=species)

## -----------------------------------------------------------------------------
#load competition data
data("gause_1934_science_f02_03")

#subset out data from species grown in mixture
mixturedat<-gause_1934_science_f02_03[gause_1934_science_f02_03$Treatment=="Mixture",]

#extract time and species data
time<-mixturedat$Day
species<-data.frame(mixturedat$Volume_Species1, mixturedat$Volume_Species2)
colnames(species)<-c("P_caudatum", "P_aurelia")

#run wrapper
gause_out<-gause_wrapper(time=time, species=species)

# number of species
N<-ncol(gause_out$rawdata)-1
# parameters
pars_full<-c(gause_out$parameter_intervals$mu)
# data.frame for optimization
fittigdata<-data.frame(y=unlist(gause_out$rawdata[,-1]),
                       time=gause_out$rawdata$time,
                       N=N)

yest<-ode_prediction(pars_full, time=fittigdata$time, N=fittigdata$N)
plot(fittigdata$y, yest, xlab="observation", ylab="prediction")
abline(a=0, b=1, lty=2)

#example of how to apply function, using nls()
mod<-nls(y~ode_prediction(pars_full, time, N),
           start = list(pars_full=pars_full),
           data=fittigdata)
summary(mod)