[R] graphics: 3D regression plane

[David Winsemius]

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