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]]
______________________________________________
[email protected] 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
______________________________________________
[email protected] 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.
