Thanks peter.


*Yao Zhu*


*Department of UrologyFudan University Shanghai Cancer CenterShanghai,
China*


2014/1/5 peter dalgaard <pdalgd@gmail.com>

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

	[[alternative HTML version deleted]]


