[R] contour plot for non-linear models

[Josep Perarnau i Codina]

Hello all,

I've tried to reproduce the contour plot that appears in the book of
Venables and Ripley, at page 255. Is a F-statistic surface and a
confidence region for the regression parameters of a non-linear model.
It uses the stormer data that are in the MASS package. 

I haven't been able to reproduce the plot either in R ( version 1.5 )
and S.  It makes the axes and it puts the labels, but it doesn't shows
the contour lines.

Could you help me resolving this problem?

This is the code that appears in the book and I've been testing.

> fm0 <- lm(Wt*Time ~ Viscosity + Time -1, data= stormer)
> b0 <- coef(fm0); names(b0) <- c("b1", "b2")
> storm.fm <- nls(Time ~ b1*Viscosity/(Wt-b2), data = stormer, start=b0,
trace = T)

> bc <- coef(storm.fm)
> se <- sqrt(diag(vcov(storm.fm)))
> dv <- deviance(storm.fm)

> par(pty = "s")
> b1 <- bc[1] + seq(-3*se[1], 3*se[1], length = 51)
> b2 <- bc[2] + seq(-3*se[2], 3*se[2], length = 51)
> bv <- expand.grid(b1,b2)
> ssq <- function(b)
+ sum((Time - b[1]*Viscosity/(Wt-b[2]))^2)
> dbetas <- apply( bv, 1 , ssq)
> attach(stormer)
> cc <- matrix(Time - rep(bv[,1],rep(23,2601)) *
+ Viscosity/(Wt - rep(bv[,2], rep(23, 2601))), 23)
> dbetas <- matrix(drop(rep(1,23) %% cc^2), 51)
> fstat <- matrix( ((dbetas -dv)/2) / (dv/21) , 51, 51)

> plot(b1, b2, type = "n")
> lev <- c(1,2,5,7,10,15,20)
> contour(b1, b2, fstat, levels= lev, labex = 0.75, lty=2, add=T)
> contour(b1, b2, fstat, levels= qf(0.95, 2,21), add=T, labex = 0)
> text(31.6, 0.3, labels = "95% CR", adj= 0, cex= 0.75)
> points(bc[1], bc[2], pch=3, mkh=0.1)
> par(pty = "m")


Josep Perarnau i Codina,      josep.perarnau at upc.es
Universitat Politècnica de Catalunya
Facultat de Matemàtiques i Estadística
LIOS - Laboratori d'Investigació Operativa i Simulació
Pau Gargallo 5                     
08028 Barcelona
Tel:  +34 93 4016947          Fax:  +34 93 4015855 


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._