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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
The current version is in MASS/scripts3/ch08.R. That runs under R. I will leave it to you to look for any differences. Without any idea what you mean by `S' I can't help you there, but the example in the book was generated using S-PLUS, and scripts are available .... On Sat, 8 Jun 2002, Josep Perarnau i Codina wrote:> 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 > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._