model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dtN <- r*N*(1 - N/K)
    return(list(c(dtN)))  
  }) 
} 

p <- c(r=1,K=100)    # Estimates close to Gause
s <- c(N=2)
run(18)              # Numerical solution of Logistic eq.

aurelia <- read.table("data/aurelia1.txt", header=TRUE)
caudatum <- read.table("data/caudatum1.txt", header=TRUE)

timePlot(aurelia,draw=points)
timePlot(caudatum,draw=points)
names(aurelia) <- c("time", "N")    # Rename the column into N 
names(caudatum) <- c("time", "N")   # Rename the column into N 

free <- c("r","K","N")
fA <- fit(aurelia,free=free,main="Aurelia",legend=FALSE)
fC <- fit(caudatum,free=free,main="Caudatum",legend=FALSE)
summary(fA)
summary(fC)

# Since estimating N(0) is difficult we use the data to fix N(0)

A0 <- aurelia[1,2]
C0 <- caudatum[1,2]
free <- c("r","K")
s <- c(N=A0); fA <- fit(aurelia,free=free,main="Aurelia",legend=FALSE)
s <- c(N=C0); fC <- fit(caudatum,free=free,main="Caudatum",legend=FALSE)
summary(fA)
summary(fC)

# Acquire confidence ranges by bootstrapping the data 100x 

s <- c(N=A0); p["r"] <- fA$par["r"]; p["K"] <- fA$par["K"]
fA <- fit(aurelia,free=free,bootstrap=100,main="Aurelia",legend=FALSE)
s <- c(N=C0); p["r"] <- fC$par["r"]; p["K"] <- fC$par["K"]
fC <- fit(caudatum,free=free,bootstrap=100,main="Caudatum",legend=FALSE)

# Since the confidence ranges on r overlap, fit both data sets together with and 
# without a shared growth rate r.
# Fixed contains the two different initial conditions taken from the data.

fB4 <- fit(list(aurelia,caudatum),fixed=list(N=c(A0,C0)),differ=c("r","K"),add=TRUE,legend=FALSE,main="4 free parameters")
fB3 <- fit(list(aurelia,caudatum),fixed=list(N=c(A0,C0)),free="r",differ="K",add=TRUE,legend=FALSE,main="3 free parameters")

# Test whether or not the improvement is significant with an F test
# p1 and p2 define the number of free parameters and n the number of data points

ftest=function(ssr1,p1,ssr2,p2,n) {
  if (p2 > p1) {
    df1 <- p2-p1
    df2 <- n-p2
    f <- ((ssr1-ssr2)/df1)/(ssr2/df2)
    cat("F[",df1,",",df2,"] =", f, ": P = ",1-pf(f,df1,df2),"\n")
  } else
    ftest(ssr2,p2,ssr1,p1,n)
}

ftest(fB4$ssr,4,fB3$ssr,3,nrow(aurelia)+nrow(caudatum))