data(crops)
ind <- 6
xv <- crops$Time[crops$Code == ind]
yv <- crops$DM[crops$Code == ind]
xlab0 <- "Time (d)"
ylab0 <- "Dry mass (g)"
dev.new()
plot(xv, yv, cex=1.5, cex.lab=1.5, cex.axis=1.5, xlab=xlab0, ylab=ylab0)
# Define the beta sigmoid model (bsm)
bsm <- function(P, x){
P <- cbind(P)
if(length(P) !=4 ) {stop("The number of parameters should be 4!")}
ropt <- P[1]
topt <- P[2]
tmin <- P[3]
tmax <- P[4]
tailor.fun <- function(x){
x[x < tmin] <- tmin
x[x > tmax] <- tmax
return(x)
}
x <- tailor.fun(x)
ropt*(x-tmin)*(x-2*tmax+topt)/(topt+tmin-2*tmax)*(
(x-tmin)/(topt-tmin) )^((topt-tmin)/(tmax-topt))
}
# For the original beta sigmoid model
ini.val0 <- c(60, 30, seq(0, 10, 20), 100)
fit1 <- fitIPEC( bsm, x=xv, y=yv, ini.val=ini.val0, xlim=NULL, ylim=NULL,
xlab=xlab0, ylab=ylab0, fig.opt=TRUE,
control=list(trace=FALSE, reltol=1e-20, maxit=50000) )
fit1$par
w <- rep(1/as.numeric(tapply(yv, xv, var)), tapply(yv, xv, length))
fit2 <- fitIPEC( bsm, x=xv, y=yv, ini.val=ini.val0, weights=w, xlim=NULL,
ylim=NULL, xlab=xlab0, ylab=ylab0, fig.opt=TRUE,
control=list(trace=FALSE, reltol=1e-20, maxit=50000) )
fit2$par
dev.new()
xp <- seq(0, 120, len=2000)
yp <- bsm(P=fit2$par, x=xp)
xv2 <- as.numeric(tapply(xv, xv, mean))
yv2 <- as.numeric(tapply(yv, xv, mean))
sd2 <- as.numeric(tapply(yv, xv, sd))
Up <- yv2+sd2
Low <- yv2-sd2
plot( xv2, yv2, xlab=xlab0, ylab=ylab0, cex.lab=1.5,
cex.axis=1.5, xlim=c(0,120), ylim=c(-5, 100), type="n" )
lines( xp, yp, col=4 )
points( xv2, yv2, pch=1, cex=1.5, col=2 )
for(i in 1:length(Up)){
lines(c(xv2[i], xv2[i]), c(Low[i], Up[i]), col=6)
}
Run the code above in your browser using DataLab