[R] graphics: 3D regression plane
David Winsemius
dwinsemius at comcast.net
Wed Jan 12 16:41:25 CET 2011
On Jan 12, 2011, at 6:10 AM, Federico Bonofiglio wrote:
> Hello Masters,
> wishing you all a great 2011 I was also going to ask if anyone knows
> a quick
> and efficient way to plot a regression plane (z~x*y).
There are many. There are limitations to using the ?? operator in that
it only brings up functions that are installed on your machine but
when I enter:
??"3D"
... on my machine it nominates a variety of functions from these
packages:
ca
car
emdbook
grDevices
HH
igraph
lattice
locfit
misc3d
plotrix
raster
rgl
rpanel
scatterplot3d
sm
sna
spatstat
spancs
TeachingDemos
vcdExtra
If you installed the sos package you would have search access to all
of the functions in CRAN packages (and maybe more).
There are also a variety of graphic galleries:
http://research.stowers-institute.org/efg/R/
http://addictedtor.free.fr/graphiques/allgraph.php
http://rgm2.lab.nig.ac.jp/RGM2/images.php?show=all&pageID=1108
>
> I have tried the regr2.plot{HH} function but it is only an
> educational tool
> and has poor graphical properties.
Ah, a critic. And a very non-specific one at that.
>
> I also tried to run the following script on a fictitious longitudinal
> problem, with poor results
Because of poor programming and failure to read the manuals.
>
> set.seed(1234)
>
> id<-c(rep(1,3),rep(2,4),rep(3,2)) # subjects
> y<-rchisq(9,df=20) #response
> k<-rnorm(9,4,2) # x
> time<-as.Date(c("03/07/1981","15/11/1981","03/04/1983","08/12/1979",
> "30/12/1979","08/03/1980","12/08/1980","12/08/1973","28/03/1975"),
> format="%d/%m/%Y")
> fac<-c("m","m","m","f","f","f","f","m","m")# sex
>
>
> d1<-as.vector(by(time,factor(id),min))
> t0<-as.Date(d1,origin=as.Date("1970-01-01"));t0
>
> A<-data.frame(id=c("1","2","3"),t=t0)
> B<-data.frame(id=id,tempo=time)
> C<-merge(A,B);C
>
> rd<-as.vector(C$tempo-C$t);rd #time centered on sbj specific first
> occurrence
>
>
> mod<-lm(y~rd*k)
> newax<- expand.grid(
> days = giorni<-seq(min(rd),max(rd), length=100),
> expl= esplic<- seq(min(k), max(k), length=100)
> )
> fit <- predict(mod,data.frame(rd=giorni,k=esplic))
> graph <- persp(x=giorni, y=esplic,fit,
> expand=0.5, ticktype="detailed", theta=-45) #error : z argument not
> valid
>
> I would be grateful if someone would give me some suggestions.
First suggestion would be to re-read the predict help page:
You are throwing together symbols in a manner not expected by predict.
The argument to newdata is invalid because you did not construct your
newax dataframe correctly, resulting in only 100 predicted points (at
the original data).
"newax" should have had column names that match the variables in the
model. This is what you got:
str(newax)
'data.frame': 10000 obs. of 2 variables:
$ days: num 0 6.45 12.91 19.36 25.82 ...
$ expl: num 0.499 0.499 0.499 0.499 0.499 ...
- attr(*, "out.attrs")=List of 2
..$ dim : Named int 100 100
.. ..- attr(*, "names")= chr "days" "expl"
..$ dimnames:List of 2
.. ..$ days: chr "days= 0.000000" "days= 6.454545" "days=
12.909091" "days= 19.363636" ...
.. ..$ expl: chr "expl=0.4985331" "expl=0.5615784"
"expl=0.6246238" "expl=0.6876691" ...
Generally is is a bad idea to use "<-" inside data.frame(). I'm not
sure if it's illegal, but it certainly is confusing. And the predict
result might have had the correct length had you had used the "newax"
dataframe, but it needed to be passed to persp as a properly
dimensioned matrix:
?persp
At the end of your constructed example try this instead:
mod<-lm(y~rd*k)
newax<- expand.grid(
rd = seq(min(rd), max(rd), length=100),
k = seq(min(k), max(k), length=100)
)
fit <- predict(mod,newax)
graph <- persp(x=seq(min(rd), max(rd), length=100),
y=seq(min(k), max(k), length=100),
z= matrix(fit, 100, 100),
expand=0.5, ticktype="detailed", theta=-45)
persp is not a lattice plotting function, so it does its plotting by
side-effects. It does return a value but it is only a transformation
matrix and I do not see that you have intentions to use i the "graph"
object, but who knows.
> Thank u again and happy new year
>
> Federico Bonofiglio
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
West Hartford, CT
More information about the R-help
mailing list