Thanks peter.
*Yao Zhu* *Department of UrologyFudan University Shanghai Cancer CenterShanghai, China* 2014/1/5 peter dalgaard <pda...@gmail.com> > > 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 > > > > > > > > > [[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.
