On 04 Jan 2014, at 13:56 , zhu yao <mailzhu...@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 <pda...@gmail.com> > > On 04 Jan 2014, at 12:53 , Uwe Ligges <lig...@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@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@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....@cbs.dk Priv: pda...@gmail.com > > > > > > > > > -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@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.
