#  Epälineaariset mallit (Ripley MAS with S-PLUS, 1999, 8. luku)

library(MASS)
library(nls)
#postscript(file="ch08.ps", width=8, height=6, pointsize=9)
#options(width=65, digits=5)

# Esimerkki

data(wtloss)
attach(wtloss)
# alter margin 4; others are default
oldpar <- par(mar=c(5.1, 4.1, 4.1, 4.1))
#plot(Days, Weight, type="p", ylab="Weight (kg)")
plot(Days,Weight, type="p", xlab="Päivät", ylab="Paino (kg)")
Wt.lbs <- pretty(range(Weight*2.205))
axis(side=4, at=Wt.lbs/2.205, lab=Wt.lbs, srt=90)
#mtext("Weight (lb)", side=4, line=3)
mtext("Paino (lb)", side=4, line=3)
par(oldpar) # restore settings
detach()


# Epälineaarisen regressiomallin estimointi

wtloss.st <- c(b0=90, b1=95, th=120)
wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th),
   data = wtloss, start = wtloss.st, trace = T)
wtloss.fm

# Menetelmäfunktioita

summary(wtloss.fm)
deviance(wtloss.fm)
vcov(wtloss.fm)   # Parametrien kovarianssimatriisi
lines(Days, predict(wtloss.fm))

#An `nls' object is a type of fitted model object.  It has methods
#     for the generic functions `coef', `formula', `resid', `print',
#     `summary', and `fitted'.

# Esimerkiksi 
fitted(wtloss.fm)
# Ks. myös help("predict")

#################

##########################
# Puromycin
##########################
library(nls)
     data(Puromycin)
#	attach(Puromycin)	
     plot(rate ~ conc, data = Puromycin, las = 1,
          xlab = "Substrate concentration (ppm)",
          ylab = "Reaction velocity (counts/min/min)",
          pch = as.integer(Puromycin$state),
          col = as.integer(Puromycin$state),
          main = "Puromycin data and fitted Michaelis-Menten curves")
# help("plot")
# help("par")
     ## simplest form of fitting the Michaelis-Menten model to these data
     fm1 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
                subset = state == "treated",
                start = c(Vm = 200, K = 0.05), trace = TRUE)
     fm2 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
                subset = state == "untreated",
                start = c(Vm = 160, K = 0.05), trace = TRUE)
     summary(fm1)
     summary(fm2)
     ## using partial linearity
     fm3 <- nls(rate ~ conc/(K + conc), data = Puromycin,
                subset = state == "treated", start = c(K = 0.05),
                algorithm = "plinear", trace = TRUE)
     ## using a self-starting model
     fm4 <- nls(rate ~ SSmicmen(conc, Vm, K), data = Puromycin,
                subset = state == "treated")
     summary(fm4)
     ## add fitted lines to the plot
     conc <- seq(0, 1.2, length = 101)
     lines(conc, predict(fm1, list(conc = conc)), lty = 1, col = 1)
     lines(conc, predict(fm2, list(conc = conc)), lty = 2, col = 2)
     legend(0.8, 120, levels(Puromycin$state),
            col = 1:2, lty = 1:2, pch = 1:2)


##########################