[R] contour plot for non-linear models

[Ko-Kang Kevin Wang]

Hi,

You should be able to if you go to $R_HOME\library\MASS\scripts3\ch08.R and
find the relevant codes.

But just in case, I simply copied it out from that file:
data(stormer)
attach(stormer)
fm0 <- lm(Wt*Time ~ Viscosity + Time - 1,  data=stormer)
b0 <- coef(fm0)
names(b0) <- c("b1", "b2")
b0
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)
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)
qf(0.95, 2, 21)
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)
detach()

------------------------------------------------
Ko-Kang Kevin Wang
Post Graduate PGDipSci Student
Department of Statistics
University of Auckland
New Zealand
www.stat.auckand.ac.nz/~kwan022

----- Original Message -----
From: "Josep Perarnau i Codina" <josep_perarnau at coeic.org>
To: <r-help at stat.math.ethz.ch>
Sent: Saturday, June 08, 2002 10:33 PM
Subject: [R] contour plot for non-linear models


> 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
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._
>


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._