Contrasts in coxph

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Contrasts in coxph

Sorkin, John
I would like to define contrasts on the output of a coxph function. It appears that the contrast function from the contrast library does not have a method defined that will allow computation of contrasts on a coxph object.

How does one define and evaluate contrasts for a cox model?

Thank you,
John

John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric Medicine
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)


______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: Contrasts in coxph

Göran Broström-3
Dear John,

I googled "contrasts Cox regression" and found

https://stats.stackexchange.com/questions/15603/performing-contrasts-among-treatment-levels-in-survival-analysis

where Frank Harrell reports how to do what you want with his rms
package. At least I think so, haven't tried it myself.

Best, Göran

On 2021-04-06 05:28, Sorkin, John wrote:

> I would like to define contrasts on the output of a coxph function.
> It appears that the contrast function from the contrast library does
> not have a method defined that will allow computation of contrasts on
> a coxph object.
>
> How does one define and evaluate contrasts for a cox model?
>
> Thank you, John
>
> John David Sorkin M.D., Ph.D. Professor of Medicine Chief,
> Biostatistics and Informatics University of Maryland School of
> Medicine Division of Gerontology and Geriatric Medicine Baltimore VA
> Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD
> 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone
> number above prior to faxing)
>
>
> ______________________________________________ [hidden email]
> mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: Contrasts in coxph

John Fox
In reply to this post by Sorkin, John
Dear John,

It's not clear to me exactly what you have in mind, but
car::linearHypothesis(), multcomp::glht(), and the emmeans package work
with Cox models. I expect there are functions in other packages that
will work too.

Here's an example, surely simpler than what you have in mind, but you
can probably adapt it:

------------------ snip -----------------

 > library("survival")
 > library("car")
Loading required package: carData
 > mod.allison <- coxph(Surv(week, arrest) ~
+                        fin + age + race + wexp + mar + paro + prio,
+                      data=Rossi)
 > mod.allison
Call:
coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp +
     mar + paro + prio, data = Rossi)

                    coef exp(coef) se(coef)      z       p
finyes         -0.37942   0.68426  0.19138 -1.983 0.04742
age            -0.05744   0.94418  0.02200 -2.611 0.00903
raceother      -0.31390   0.73059  0.30799 -1.019 0.30812
wexpyes        -0.14980   0.86088  0.21222 -0.706 0.48029
marnot married  0.43370   1.54296  0.38187  1.136 0.25606
paroyes        -0.08487   0.91863  0.19576 -0.434 0.66461
prio            0.09150   1.09581  0.02865  3.194 0.00140

Likelihood ratio test=33.27  on 7 df, p=2.362e-05
n= 432, number of events= 114
 >
 > linearHypothesis(mod.allison, "finyes")
Linear hypothesis test

Hypothesis:
finyes = 0

Model 1: restricted model
Model 2: Surv(week, arrest) ~ fin + age + race + wexp + mar + paro + prio

   Res.Df Df  Chisq Pr(>Chisq)
1    426
2    425  1 3.9306    0.04742 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 >
 > library("multcomp")
Loading required package: mvtnorm
Loading required package: TH.data
Loading required package: MASS

Attaching package: ‘TH.data’

The following object is masked from ‘package:MASS’:

     geyser

 > summary(glht(mod.allison, "finyes=0"))

         Simultaneous Tests for General Linear Hypotheses

Fit: coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp +
     mar + paro + prio, data = Rossi)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)
finyes == 0  -0.3794     0.1914  -1.983   0.0474 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

 >
 > library(emmeans)
 > pairs(emmeans(mod.allison, ~ fin))
  contrast estimate    SE  df z.ratio p.value
  no - yes    0.379 0.191 Inf 1.983   0.0474

------------------ snip -----------------

Results are averaged over the levels of: race, wexp, mar, paro
Results are given on the log (not the response) scale.

John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
web: https://socialsciences.mcmaster.ca/jfox/

On 2021-04-05 11:28 p.m., Sorkin, John wrote:

> I would like to define contrasts on the output of a coxph function. It appears that the contrast function from the contrast library does not have a method defined that will allow computation of contrasts on a coxph object.
>
> How does one define and evaluate contrasts for a cox model?
>
> Thank you,
> John
>
> John David Sorkin M.D., Ph.D.
> Professor of Medicine
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology and Geriatric Medicine
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
>
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: Contrasts in coxph

David Winsemius
In reply to this post by Sorkin, John

On 4/5/21 8:28 PM, Sorkin, John wrote:
> I would like to define contrasts on the output of a coxph function. It appears that the contrast function from the contrast library does not have a method defined that will allow computation of contrasts on a coxph object.
>
> How does one define and evaluate contrasts for a cox model?


I suspect that you are hoping for post-hoc pairwise contrasts. If that's
the case, then the multcomp package's glht is generally very capable. It
has three examples in its vignette the involve survivla models. The last
one constructs a categorical variable for age. That's not a procedure
that I endorse, but it may have been done for illustration of a method
and would therefore be free of sin.


library(multcomp)  # console output follows

 > if (require("survival") && require("MASS")) {
+     ### construct 4 classes of age
+     Melanoma$Cage <- factor(sapply(Melanoma$age, function(x){if( x <=
25 ) return(1)
+         if( x > 25 & x <= 50 ) return(2)
+         if( x > 50 & x <= 75 ) return(3)
+         if( x > 75 & x <= 100) return(4) }
+     ))
+ }
 > summary(Melanoma)
       time          status          sex age             year       
thickness
  Min.   :  10   Min.   :1.00   Min.   :0.0000   Min.   : 4.00 Min.  
:1962   Min.   : 0.10
  1st Qu.:1525   1st Qu.:1.00   1st Qu.:0.0000   1st Qu.:42.00 1st
Qu.:1968   1st Qu.: 0.97
  Median :2005   Median :2.00   Median :0.0000   Median :54.00 Median
:1970   Median : 1.94
  Mean   :2153   Mean   :1.79   Mean   :0.3854   Mean   :52.46 Mean  
:1970   Mean   : 2.92
  3rd Qu.:3042   3rd Qu.:2.00   3rd Qu.:1.0000   3rd Qu.:65.00 3rd
Qu.:1972   3rd Qu.: 3.56
  Max.   :5565   Max.   :3.00   Max.   :1.0000   Max.   :95.00 Max.  
:1977   Max.   :17.42
      ulcer       Cage
  Min.   :0.000   1: 14
  1st Qu.:0.000   2: 73
  Median :0.000   3:104
  Mean   :0.439   4: 14
  3rd Qu.:1.000
  Max.   :1.000
 > cm <- coxph(Surv(time, status == 1) ~ Cage, data = Melanoma)
 > ### specify all pair-wise comparisons among levels of "Cage"
 > cm.glht <- glht(cm, mcp(Cage = "Tukey"))
 > cm.glht

      General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Linear Hypotheses:
            Estimate
2 - 1 == 0  -0.2889
3 - 1 == 0   0.1910
4 - 1 == 0   1.0742
3 - 2 == 0   0.4800
4 - 2 == 0   1.3631
4 - 3 == 0   0.8832


--

David.

>
> Thank you,
> John
>
> John David Sorkin M.D., Ph.D.
> Professor of Medicine
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology and Geriatric Medicine
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
>
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> 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.

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.