[R] Very basic question regarding plot.Design...

[Frank E Harrell Jr]

I just added a new option nlines to plot.Predict to make it easy to get 
interaction line plots for categorical predictors.  To get the new 
version type the following after running require(rms):

source('http://biostat.mc.vanderbilt.edu/tmp/plot.Predict.s')

Here are some examples:

require(rms)
n <- 100
set.seed(1)
gender <- c(rep('male', n), rep('female',n))
m <- sample(c('a','b'), 2*n, TRUE)
d <-  datadist(gender, m); options(datadist='d')
anxiety <- runif(2*n) + .2*(gender=='female') + .4*(gender=='female' & 
m=='b')
tapply(anxiety, llist(gender,m), mean)
f <- ols(anxiety ~ gender*m)
p <- Predict(f, gender=., m=.)
plot(p)     # horizontal dot chart; usually preferred for categorical 
predictors
Key(.5, .5)
plot(p, ~gender, groups='m', nlines=TRUE)
plot(p, ~m, groups='gender', nlines=TRUE)
plot(p, ~gender|m, nlines=TRUE)


Frank

Petar Milin wrote:
> Sorry, I hope this will be the last (for now, at least).
> 
> Following your advices, I did:
>  n <- 100
>  anxiety <- c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2),
>   rnorm(100, 25, 3.1), rnorm(100, 10, 2.6))
>  m.status <- c(rep('married', n*2), rep('not married', n*2))
>  gender  <- c(rep('male', n), rep('female', n),
>   rep('male', n), rep('female',n))
> 
>  im <- as.integer(p2$m.status)
>  m.status <- as.factor(m.status)
>  gender <- as.factor(gender)
> 
>  require(rms)
> 
>  d <- datadist(gender, m.status); options(datadist='d')
> 
>  ols1 <- ols(anxiety ~ m.status)
>  p1 <- Predict(ols1, m.status=.)
>  ols2 <- ols(anxiety ~ m.status * gender)
>  p2 <- Predict(ols2, m.status=., gender=.)
> 
>  pdf('figs/anova.pdf', height=6, width=8)
>  par(mfrow=c(1,2), mar=c(4,4,1,1), cex=1.2)
>  with(p1, plot(1:2, yhat, type='l', xlab='marital status',
>   ylab='anxiety', ylim=c(5,30), lwd=1.2))
>  xyplot(yhat ~ im, groups=gender,
>   type='l', data=p2, ylim=c(5,30), xlab='marital status',
>   ylab='anxiety', scales=list(x=list(at=1:2,
>   labels=levels(p2$m.status))), col=c('black','black'),
>   lwd=1.2, label.curve=list(offset=unit(.3,"in")))
>  par(mfrow=c(1,1))
>  dev.off()
> 
> Graphs are great! However, now, they does not appear in one panel.
> 
> Sorry, again, for asking so many questions. I am trying to wrap this.
> 
> Best,
> PM
> 


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University