I'm new to R (and to S) and am wondering about code from pages 72 and
83 of MASS (Venables+Ripley, 3rd edition), to draw an envelope on a QQ
plot. Copying from the book, I've got:
#... code whose gist is "a.fit <- nls(..."
num.points <- length(resid(a.fit))
qqnorm(residuals(a.fit)) # illustrate data-model residuals
qqline(residuals(a.fit))
samp <- cbind(residuals(a.fit), matrix(rnorm(num.points*19),num.points,19))
samp <- apply(scale(samp), 2, sort)
rs <- samp[,1]
xs <- qqnorm(rs, plot=FALSE)$x
env <- t(apply(samp[,-1], 1, range)) ###########################3
xyul <- par("usr")
smidge <- min(diff(c(xyul[1], xs, xyul[2])))/2
segments(xs-smidge,env[,1], xs+smidge, env[,1])
segments(xs-smidge,env[,2], xs+smidge, env[,2])
where I've marked a confusing line with ########. From what I gather
(from section 2.1, or page 5, of "R complements to MASS" by VR),
indexing is different in R than in S. It's not clear to me whether
this R-to-S difference is giving me problems. I'm such a newbie that
I'm not entirely sure what the 'samp[,-1]' is supposed to be doing.
(Confession: I'm such a newbie that I've spent much of the day typing
"?apply" and then "?par" etc, working through the program to
try to
figure it out!)
In any case, my worry is visible in the QQ plot attached. Why do I
have get so few "ledges" from the segments() calls? It's not
clear to
me how the scale of the data-model deviations come into this problem.
Q: should I first be standardizing the residuals from the data-model
comparison? I would have guessed that something in this 'QQ envelope'
would be doing that.
I'm guessing any replies will fit into one of three categories ...
1) You idiot! QQ plots are supposed to have standardized residuals of
the data-model misfit. Do that and the QQ envelop is ok.
2) Don't copy S stuff into R code. The S code "samp[,-1]" from VR
should be written "..." for use in R.
3) Don't be so rude as to ask such questions (with length attachments,
no less) on a group like this. Give us a break!
... but I'm hoping that the mode won't be in the final category!
Many thanks in advance.
Dan E. Kelley internet: mailto:Dan.Kelley at Dal.CA
Oceanography Department phone: (902)494-1694
Dalhousie University fax: (902)494-2885
Halifax, NS, CANADA, B3H 4J1 http://www.phys.ocean.dal.ca/~kelley
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplots.ps.gz
Type: application/octet-stream
Size: 2860 bytes
Desc:
Url :
https://stat.ethz.ch/pipermail/r-help/attachments/20000224/7d190b34/Rplots.ps.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fit.R.gz
Type: application/octet-stream
Size: 1657 bytes
Desc:
Url :
https://stat.ethz.ch/pipermail/r-help/attachments/20000224/7d190b34/fit.R.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fbox.dat.gz
Type: application/octet-stream
Size: 380 bytes
Desc:
Url :
https://stat.ethz.ch/pipermail/r-help/attachments/20000224/7d190b34/fbox.dat.obj
On Thu, 24 Feb 2000, Dan E. Kelley wrote:> I'm new to R (and to S) and am wondering about code from pages 72 and > 83 of MASS (Venables+Ripley, 3rd edition), to draw an envelope on a QQ > plot. Copying from the book, I've got: > > #... code whose gist is "a.fit <- nls(..." > num.points <- length(resid(a.fit)) > qqnorm(residuals(a.fit)) # illustrate data-model residuals > qqline(residuals(a.fit)) > samp <- cbind(residuals(a.fit), matrix(rnorm(num.points*19),num.points,19)) > samp <- apply(scale(samp), 2, sort) > rs <- samp[,1] > xs <- qqnorm(rs, plot=FALSE)$x > env <- t(apply(samp[,-1], 1, range)) ###########################3 > xyul <- par("usr") > smidge <- min(diff(c(xyul[1], xs, xyul[2])))/2 > segments(xs-smidge,env[,1], xs+smidge, env[,1]) > segments(xs-smidge,env[,2], xs+smidge, env[,2]) > > where I've marked a confusing line with ########. From what I gather > (from section 2.1, or page 5, of "R complements to MASS" by VR), > indexing is different in R than in S. It's not clear to me whether > this R-to-S difference is giving me problems. I'm such a newbie that > I'm not entirely sure what the 'samp[,-1]' is supposed to be doing. > (Confession: I'm such a newbie that I've spent much of the day typing > "?apply" and then "?par" etc, working through the program to try to > figure it out!)You can look in the R scripts (in the MASS package in the VR bundle) to see fully tested R versions of the code. That works, although you left out a line (starting matplot). samp[, -1] drops the first column (pp.39-40), so the marked line applies range to the remaining columns. That code applies to a variable, not residuals, which may or may not matter depending on your model complexity.> In any case, my worry is visible in the QQ plot attached. Why do I > have get so few "ledges" from the segments() calls? It's not clear to > me how the scale of the data-model deviations come into this problem. > > Q: should I first be standardizing the residuals from the data-model > comparison? I would have guessed that something in this 'QQ envelope' > would be doing that.Yes, the `scale' function is already normalizing for you. -- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
"Dan E. Kelley" <kelley at Phys.Ocean.Dal.CA> writes:> Content-Type: TEXT/PLAIN; charset=US-ASCII > > I'm new to R (and to S) and am wondering about code from pages 72 and > 83 of MASS (Venables+Ripley, 3rd edition), to draw an envelope on a QQ > plot. Copying from the book, I've got: > > #... code whose gist is "a.fit <- nls(..." > num.points <- length(resid(a.fit)) > qqnorm(residuals(a.fit)) # illustrate data-model residuals > qqline(residuals(a.fit)) > samp <- cbind(residuals(a.fit), matrix(rnorm(num.points*19),num.points,19)) > samp <- apply(scale(samp), 2, sort) > rs <- samp[,1] > xs <- qqnorm(rs, plot=FALSE)$x > env <- t(apply(samp[,-1], 1, range)) ###########################3 > xyul <- par("usr") > smidge <- min(diff(c(xyul[1], xs, xyul[2])))/2 > segments(xs-smidge,env[,1], xs+smidge, env[,1]) > segments(xs-smidge,env[,2], xs+smidge, env[,2]) >....> 1) You idiot! QQ plots are supposed to have standardized residuals of > the data-model misfit. Do that and the QQ envelop is ok. > > 2) Don't copy S stuff into R code. The S code "samp[,-1]" from VR > should be written "..." for use in R. > > 3) Don't be so rude as to ask such questions (with length attachments, > no less) on a group like this. Give us a break! > > ... but I'm hoping that the mode won't be in the final category!Well, 2) is *not* the answer... Try removing the plot=F from the 2nd qqnorm and I think you'll see the light. Notice the scale on the y axis. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._