graphics: 3D regression plane

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