[R] p value of trends for odds ratios (or hazard ratios)

Sat Jan 4 19:56:05 CET 2014

On 04 Jan 2014, at 13:56 , zhu yao <mailzhuyao at gmail.com> wrote:

> Thanks for the suggestion.
> The results is presented in following table. The authors calculated p value for linearity and trend and they stated in the methods:
> " Linear and nonlinear trends of BMI associated with each mortality outcome were obtained by modeling BMI as a continuous variable and using the partial likelihood test for linearity. To test whether the associations between each adiposity measure and mortality were modified by race/ethnicity, we constructed a likelihood ratio test for heterogeneity of trends comparing 2 multivariate Cox proportional hazard models."
> 

Unfortunately, there is not enough information in the table to reconstruct the p-values. 

Apparently, those come from a Cox regression analysis, and given the discrepancy between the rates per 1000 years and the fully adjusted HR, there is little hope to get anywhere close using approximate methods. Notice, for instance that the rates of 14.61 and 22.04 have nearly the same HR (1.41 and 1.42), which presumably is due to the correction for race/ethnicity - one might conjecture that raw results are skewed by an excess of (say) underweight Asians and overweight Latin-Americans.  

If all the information was the table of person-years and incidence counts, you might assume a constant-hazard model and do a Poisson regression, like this

> count <- c(34,602,429,238,81,61)
> pyr <- c(2328, 59094, 37687, 16832, 5789, 2768)
> obese <- 1:6
> obese_f <- relevel(factor(obese),2)

> ### Fit full model with obese as a factor
> ### (Drop the "0 +" to get rate ratios with cat.2 as reference)
> fit1 <- glm(count ~ 0 + obese_f + offset(log(pyr)),family=poisson)

> ### This reproduces the rates:
> exp(coef(fit1))*1000
obese_f2 obese_f1 obese_f3 obese_f4 obese_f5 obese_f6 
10.18716 14.60481 11.38324 14.13973 13.99205 22.03757

> ### fit model with obese as numerical variable
> fit2 <- glm(count~obese +offset(log(pyr)),family=poisson)

> ### fit model without obese
> fit3 <- glm(count~offset(log(pyr)),family=poisson)

> ### trend test
> anova(fit3,fit2, test="Rao")
Analysis of Deviance Table

Model 1: count ~ offset(log(pyr))
Model 2: count ~ obese + offset(log(pyr))
  Resid. Df Resid. Dev Df Deviance    Rao  Pr(>Chi)    
1         5     44.365                                 
2         4     11.427  1   32.939 34.776 3.699e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> ### test for linearity
> anova(fit2,fit1, test="Rao")
Analysis of Deviance Table

Model 1: count ~ obese + offset(log(pyr))
Model 2: count ~ 0 + obese_f + offset(log(pyr))
  Resid. Df Resid. Dev Df Deviance    Rao Pr(>Chi)  
1         4     11.427                              
2         0      0.000  4   11.427 12.908  0.01173 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

but as you see, the results are markedly different from what you get in the adjusted analysis.

This is not Stata "tabodds" territory, more like "strate", but I don't see a trend test option for that (or any test option for that matter), so I expect that Stata-ists would go directly for Poisson regression too. (Jim Holtman's tag line is very appropriate here)

To get the results in the table, I suppose that you would use the full data and do similar model comparisons using 

fit1 <- coxph(Surv(time, event) ~ ethnicity + obese_f)
fit2 <- coxph(Surv(time, event) ~ ethnicity + obese)
fit3 <- coxph(Surv(time, event) ~ ethnicity)

-pd

> <image.png>
> 
> Yao Zhu
> Department of Urology
> Fudan University Shanghai Cancer Center
> Shanghai, China
> 
> 
> 2014/1/4 peter dalgaard <pdalgd at gmail.com>
> 
> On 04 Jan 2014, at 12:53 , Uwe Ligges <ligges at statistik.tu-dortmund.de> wrote:
> 
> >
> >
> > On 04.01.2014 06:39, zhu yao wrote:
> >> Dear Sir
> >> Many papers calculated the p value of trends for odds ratios of ordered
> >> category variables. I have found the tabodds command in Stata. But how to
> >> do it in R?
> >
> > Depends on the method you want to use ... and most ladies and gents on this list won't know Stata well enough to know what tabodds does.
> 
> We can google its help files, though.
> 
> I would expect that prop.trend.test is pretty much equivalent, or in general score tests in the appropriate logistic regression model. However, a worked example in Stata and some reproducible data would allow us to say with more confidence.
> 
> >
> > Best,
> > Uwe Ligges
> >
> >
> >
> >
> >> Thanks
> >>
> >> *Yao Zhu*
> >>
> >>
> >> *Department of UrologyFudan University Shanghai Cancer CenterShanghai,
> >> China*
> >>
> >>      [[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.
> >>
> >
> > ______________________________________________
> > 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.
> 
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
> 
> 
> 
> 
> 
> 
> 
> 
> 

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com