### R code from vignette source 'curve_1.Rnw'

###################################################
### code chunk number 1: curve_1.Rnw:7-11
###################################################
options(prompt="> ", continue="+ ", digits=7, width=70, show.signif.stars=T)
par(mfrow=c(1,1))
rm(list=ls())
setwd("~/data/Journal/Monographs/StatsMonographs/CurveFit/")


###################################################
### code chunk number 2: curve_1.Rnw:64-65
###################################################
set.seed(1485)


###################################################
### code chunk number 3: curve_1.Rnw:68-76
###################################################
len <- 24
x <- runif(len)
y <- x^3 + rnorm(len, 0, .06)
ds <- data.frame(x=x, y=y)
str(ds)
plot(y ~ x, main="Known cubic, with noise")
s <- seq(0, 1, length=100)
lines(s, s^3, lty=2, col="green")


###################################################
### code chunk number 4: curve_1.Rnw:112-114
###################################################
m <- nls(y ~ I(x^power), data=ds, start=list(power=1), trace=T)
class(m)


###################################################
### code chunk number 5: curve_1.Rnw:125-127
###################################################
summary(m)
summary(m)$coefficients


###################################################
### code chunk number 6: curve_1.Rnw:141-148
###################################################
power <- round(summary(m)$coefficients[1],3)
power.se <- round(summary(m)$coefficients[2],3)
plot(y ~ x, main="Fitted power model", sub="Blue: fit; green: known")
s <- seq(0, 1, length=100)
lines(s, s^3, lty=2, col="green")
lines(s, predict(m, list(x=s)), lty=1, col="blue")
text(0, 0.5, paste("y =x^ (", power, " +/- ", power.se, ")", sep=""),pos=4)


###################################################
### code chunk number 7: curve_1.Rnw:164-167
###################################################
(RSS.p <- sum(residuals(m)^2))
(TSS <- sum((y - mean(y))^2))
1-(RSS.p/TSS)


###################################################
### code chunk number 8: curve_1.Rnw:173-174
###################################################
1 - sum((x^3 - y)^2)/TSS


###################################################
### code chunk number 9: aic
###################################################
AIC(m)


###################################################
### code chunk number 10: curve_1.Rnw:239-248
###################################################
rhs <- function(x, b0, b1) { b0 + x^b1}
m.2 <- nls(y ~ rhs(x, intercept, power) , data=ds, start=list(intercept=0, power=2), trace=T)
summary(m.2)
plot(ds$y ~ ds$x, main="Fitted power model, with intercept", sub="Blue: fit; magenta: fit w/o intercept; green: known")
abline(h=0, lty=1, lwd=0.5)
lines(s, s^3, lty=2, col="green")
lines(s, predict(m.2, list(x=s)), lty=1, col="blue")
lines(s, predict(m, list(x=s)), lty=2, col="magenta")
segments(x, y, x, fitted(m.2), lty=2, col="red")


###################################################
### code chunk number 11: curve_1.Rnw:257-261
###################################################
(RSS.pi <- sum(residuals(m.2)^2))
(r2.pi <- (1 - (RSS.pi/TSS)))
(r2.p <- 1 - (RSS.p/TSS))
(r2.3 <- 1 - sum((x^3 - y)^2)/TSS)


###################################################
### code chunk number 12: adj-r2
###################################################
(n <- dim(ds)[1])
(p <- length(coefficients(m.2)))
(r2.adj.pi <- 1 - (((n-1)/(n-p)) * (1-r2.pi)))
r2.pi
r2.pi - r2.adj.pi


###################################################
### code chunk number 13: curve_1.Rnw:301-302
###################################################
anova(m.2, m)


###################################################
### code chunk number 14: curve_1.Rnw:318-320
###################################################
AIC(m.2)
AIC(m)


###################################################
### code chunk number 15: curve_1.Rnw:349-350
###################################################
min(y)


###################################################
### code chunk number 16: curve_1.Rnw:357-364
###################################################
offset <- 0.1
ds$ly <- log(ds$y+offset)
m.l <- lm(ly ~ x, data=ds)
summary(m.l)
plot(ds$ly ~ ds$x, xlab="x", ylab="log(y+.1)", main="Log-linear fit")
abline(m.l)
text(0, 0.4, pos=4, paste("log(y) = ",round(coefficients(m.l)[1],3),"+",round(coefficients(m.l)[2],3)))


###################################################
### code chunk number 17: curve_1.Rnw:374-377
###################################################
par(mfrow=c(2,2))
plot(m.l)
par(mfrow=c(1,1))


###################################################
### code chunk number 18: curve_1.Rnw:391-400
###################################################
m.e <- nls(y ~ I(exp(1)^(a + b*x)), data=ds, start=list(a=0, b=1), trace=T)
summary(m.e)$coefficients
a <- round(summary(m.e)$coefficients[1,1],4)
b <- round(summary(m.e)$coefficients[2,1],4)
plot(y ~ x, main="Fitted exponential function", sub="Blue: fit; green: known")
s <- seq(0, 1, length=100)
lines(s, s^3, lty=2, col="green")
lines(s, predict(m.e, list(x=s)), lty=1, col="blue")
text(0, 0.5, paste("y =e^ (", a, " + ", b, " * x)", sep=""),pos=4)


###################################################
### code chunk number 19: curve_1.Rnw:405-410
###################################################
RSS.p
(RSS.e <- sum(residuals(m.e)^2))
TSS
1- RSS.p/TSS
1- RSS.e/TSS


###################################################
### code chunk number 20: curve_1.Rnw:417-419
###################################################
AIC(m)
AIC(m.e)


###################################################
### code chunk number 21: curve_1.Rnw:439-442
###################################################
f.lrp <- function(x, a, b, t.x) {
  ifelse(x > t.x, a + b*t.x, a + b*x)
}


###################################################
### code chunk number 22: curve_1.Rnw:452-459
###################################################
f.lvls <- seq(0,120,by=10)
a.0 <- 2; b.0 <- 0.05; t.x.0 <- 70
test <- data.frame(x = f.lvls, y = f.lrp(f.lvls, a.0, b.0, t.x.0))
test <- rbind(test, test, test)
set.seed <- 1040
test$y <- test$y + rnorm(length(test$y), 0,.2)
str(test)


###################################################
### code chunk number 23: curve_1.Rnw:468-473
###################################################
plot(test$y ~ test$x, main="Linear response and plateau yield response", xlab="Fertilizer added", ylab="Crop yield")
(max.yield <- a.0+b.0*t.x.0)
lines(x =c(0, t.x.0, 120), y = c(a.0, max.yield, max.yield), lty=2)
abline(v=t.x.0, lty=3)
abline(h=max.yield, lty=3)


###################################################
### code chunk number 24: curve_1.Rnw:479-482
###################################################
test$rep <- as.factor(rep(1:3, each=length(test$y)/3))
str(test)
table(test$rep)


###################################################
### code chunk number 25: curve_1.Rnw:488-489
###################################################
by(test$y, test$rep, mean)


###################################################
### code chunk number 26: curve_1.Rnw:498-501
###################################################
m.lrp <- nls(y ~ f.lrp(x, a, b, t.x), data=test, start=list(a=0, b=.1, t.x=50), trace=T,control=list(warnOnly=T, minFactor=1/2048))
summary(m.lrp)
coefficients(m.lrp)


###################################################
### code chunk number 27: r2.lrp
###################################################
(RSS.lrp <- sum(residuals(m.lrp)^2))
(TSS <- sum((test$y - mean(test$y))^2))
(r2.m.lrp <- 1-(RSS.lrp/TSS))
(n <- dim(test)[1])
(p <- length(coefficients(m.lrp)))
(r2.adj.m.lrp <- 1 - (((n-1)/(n-p)) * (1-r2.m.lrp)))
r2.m.lrp - r2.adj.m.lrp


###################################################
### code chunk number 28: curve_1.Rnw:524-535
###################################################
plot(test$y ~ test$x, main="Linear response and plateau yield response", xlab="Fertilizer added", ylab="Crop yield")
(max.yield <- a.0+b.0*t.x.0)
lines(x =c(0, t.x.0, 120), y = c(a.0, max.yield, max.yield), lty=2, col="blue")
abline(v=t.x.0, lty=3, col="blue")
abline(h=max.yield, lty=3, col="blue")
(max.yield <- coefficients(m.lrp)["a"] + coefficients(m.lrp)["b"]*coefficients(m.lrp)["t.x"])
lines(x =c(0, coefficients(m.lrp)["t.x"], 120), y = c(coefficients(m.lrp)["a"], max.yield, max.yield), lty=1)
abline(v=coefficients(m.lrp)["t.x"], lty=4)
abline(h=max.yield, lty=4)
text(120, 4, "known true model", col="blue", pos=2)
text(120, 3.5, "fitted model", col="black", pos=2)


###################################################
### code chunk number 29: lrp.AIC
###################################################
m.quad <- lm(y ~ I(x^2) + x, data=test)
summary(m.quad)
summary(m.quad)$adj.r.squared
r2.adj.m.lrp
r2.adj.m.lrp - summary(m.quad)$adj.r.squared
AIC(m.quad)
AIC(m.lrp)