[R] contour plot for non-linear models

[Bill.Venables@cmis.csiro.au]

The reply from Kevin should work, but you need one more prior step:

	library(nls)

or, to be a bit more fussy,

	require(nls)

You don't need this in S-PLUS, of course.

Bill Venables.

>  -----Original Message-----
> From: 	Ko-Kang Kevin Wang [mailto:Ko-Kang at xtra.co.nz] 
> Sent:	Saturday, June 08, 2002 9:08 PM
> To:	Josep Perarnau i Codina; r-help at stat.math.ethz.ch
> Subject:	Re: [R] contour plot for non-linear models
> 
> 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
> _._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._