[R] How to compare the effect of a variable across regression models?

[Greg Snow]

If you just want to visualize the effect on one variable on the response from some different models then you might try Predict.Plot from the TeachingDemos package.  It takes a little tweaking to get it to work with cph objects, but here is a basic example (partly stolen from the help page for cph):

library(rms)
library(TeachingDemos)


     set.seed(731)
     age <- 50 + 12*rnorm(n)
     label(age) <- "Age"
     sex <- factor(sample(c('Male','Female'), n, 
                   rep=TRUE, prob=c(.6, .4)))
     cens <- 15*runif(n)
     h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
     dt <- -log(runif(n))/h
     label(dt) <- 'Follow-up Time'
     e <- ifelse(dt <= cens,1,0)
     dt <- pmin(dt, cens)
     units(dt) <- "Year"
  

tmp.df <- data.frame(dt=dt, e=e, age=age, sex=sex)
f <- cph(Surv(dt,e) ~ rcs(age,4) + sex, data=tmp.df )
f$data <- tmp.df
Predict.Plot(f, 'age', age=50, sex='Male', type='lp', 
	plot.args=list(col='blue',ylim=c(-1,2)))
Predict.Plot(f, 'age', age=50, sex='Female', type='lp', add=TRUE,
	plot.args=list(col='red'))

What that all means depends on the things mentioned in the other replies.  Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Biau David
> Sent: Friday, August 13, 2010 7:42 AM
> To: r help list
> Subject: [R] How to compare the effect of a variable across regression
> models?
> 
> Hello,
> 
> I would like, if it is possible, to compare the effect of a variable
> across
> regression models. I have looked around but I haven't found anything.
> Maybe
> someone could help? Here is the problem:
> 
> I am studying the effect of a variable (age) on an outcome (local
> recurrence:
> lr). I have built 3 models:
> - model 1: lr ~ age      y = \beta_(a1).age
> - model 2: lr ~ age +  presentation variables (X_p)        y =
> \beta_(a2).age +
> \BETA_(p2).X_p
> - model 3: lr ~ age + presentation variables + treatment variables(
> X_t)
>        y = \beta_(a3).age  + \BETA_(p3).X_(p) + \BETA_(t3).X_t
> 
> Presentation variables include variables such as tumor grade, tumor
> size, etc...
> the physician cannot interfer with these variables.
> Treatment variables include variables such as chemotherapy, radiation,
> surgical
> margins (a surrogate for adequate surgery).
> 
> I have used cph for the models and restricted cubic splines (Design
> library) for
> age. I have noted that the effect of age decreases from model 1 to 3.
> 
> I would like to compare the effect of age on the outcome across the
> different
> models. A test of \beta_(a1) = \beta_(a2) = \beta_(a3) and then two by
> two
> comparisons or a global trend test maybe? Is that possible?
> 
> Thank you for your help,
> 
> 
> David Biau.
> 
> 
> 
> 
> 	[[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.