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

Frank E Harrell Jr f.harrell at vanderbilt.edu
Wed Sep 9 15:02:26 CEST 2009


Petar Milin wrote:
> Thank you very much for the help!
> Now, I have an additional question regarding transcending from Design to 
> rms: Before, it was possible to plot interaction, and get lines with
>    plot(ols2, gender=NA, marital.status=NA, xlab='gender',
>    ylab='anxiety', conf.int=FALSE)
> Now, I am lost how to do that. If model is:
>    ols2 <- ols(anxiety ~ marital.status * gender)
> Then, I need:
>    p2 <- Predict(ols2, marital.status=., gender=.)
> And plot(p2) does not give me lines, but points with conf.int. Also, 
> axes are transposed as to what old plot(ols2) would give.

You're right, the plot is not optimal for the case of two factor 
predictors.  Here is a way to take control:

m <- sample(c('a','b'), 2*n, TRUE)
d <-  datadist(d, m)
f <- ols(anxiety ~ gender*m)
p <- Predict(f, m=., gender=.)
plot(p)  # default plot
Key(.5, .5)

im <- as.integer(p$m)
xYplot(Cbind(yhat,lower,upper) ~ im, groups=gender,
        type='l', data=p, ylim=c(5,30), xlab='m',
        scales=list(x=list(at=1:2, labels=levels(p$m))),
        label.curve=list(offset=unit(.3,"in")))

For the future I'll see if I can incorporate this automatically in 
plot.Predict with there are only 2 predictors being varied and they are 
both factors.


> In addition, if I try old function (plot(ols2), above), error message 
> returns:
>    Error in predictDesign(fit, adj, type = "x", non.slopes = non.slopes)
>    could not find function "Varcov"

That is because we need to update Design to be compatible with the 
latest Hmisc.

Frank

> 
> Best,
> PM
> 
> Frank E Harrell Jr wrote:
>> Petar Milin wrote:
>>> I would like to have a line on this plot, instead of two points:
>>>
>>> x1 = rnorm(100, 10, 2.5)
>>> x2 = rnorm(100, 26, 3.2)
>>> x1 = as.data.frame(x1)
>>> x2 = as.data.frame(x2)
>>> colnames(x1) = 'anxiety'
>>> colnames(x2) = 'anxiety'
>>> x1$gender = 'male'
>>> x2$gender = 'female'
>>> dat = rbind(x1, x2)
>>>
>>> require(Design)
>>>
>>> attach(dat)
>>> d=datadist(gender)
>>> options(datadist="d")
>>>
>>> ols1 <- ols(anxiety ~ gender, data=dat, x=T, y=T)
>>>
>>> plot(ols1, gender=NA, xlab="gender", ylab="anxiety",
>>> ylim=c(5,30), conf.int=FALSE)
>>>
>>> detach(dat)
>>
>> Your code has many problems and inefficiencies in it.  Here are some 
>> suggested fixes and the commands needed using the new rms package:
>>
>> require(rms)
>> n <- 100
>> anxiety <- c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2))
>> gender  <- c(rep('male', n), rep('female',n))
>> d <- datadist(gender); options(datadist='d')
>> f <- ols(anxiety ~ gender)
>> p <- Predict(f, gender=.)
>> # For this case could also do p <- Predict(f); plot(p) would give a
>> # vertical dot chart
>> p           # print estimates
>> plot(p)     # horizontal dot chart; preferred for categorical predictors
>> # To take control using lines:
>> with(p, plot(1:2, yhat, type='l', xlab='gender numeric'))
>>
>> Frank
>>
>>>
>>> Thanks!
>>> PM
>>>
>>> Frank E Harrell Jr wrote:
>>>> Petar Milin wrote:
>>>>> Hello ALL!
>>>>> I have a problem to plot factor (lets say gender) as a line, or at 
>>>>> least both line and point, from ols model:
>>>>> ols1 <- ols(Y ~ gender, data=dat, x=T, y=T)
>>>>> plot(ols1, gender=NA, xlab="gender", ylab="Y",
>>>>> ylim=c(5,30), conf.int=FALSE)
>>>>>
>>>>> If I convert gender into discrete numeric predictor, and use 
>>>>> forceLines=TRUE, plot is not nice and true, since it shows values 
>>>>> between 1 and 2.
>>>>>
>>>>> Thanks!
>>>>> PM
>>>>>
>>>>
>>>> Petar,
>>>>
>>>> forceLines seems to be doing what it was intended to do.  I'm not 
>>>> clear on why you need a line, though.  If you provide self-contained 
>>>> code and data that replicate your problem I may be able to help 
>>>> more, or you can try a new package I'm about to announce.
>>>>
>>>> Frank
>>>
>>
>>
> 


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




More information about the R-help mailing list