Jacob Brogren
2011-Jun-28 10:51 UTC
[R] coxph() - unexpected result using Crawley's seedlings data (The R Book)
Hi, I ran the example on pp. 799-800 from Machael Crawley's "The R Book" using package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The model is a Cox's Proportional Hazards model. The result was quite different compared to the R Book. I have compared my code to the code in the book but can not find any differences in the function call. My results are attached as well as a link to the results presented in the book (link to Google Books). When running the examples on pp. 797-799 I can't detect any differences in results so I don't think there are errors in the data set or in the creation of the status variable. --------------------------------- Original from the R Book: http://books.google.com/books?id=8D4HVx0apZQC&lpg=PA799&ots=rQgd_8ofeS&dq=r%20coxph%20crawley&pg=PA799#v=onepage&q&f=false --------------------------------- My result:> summary(model1)Call: coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize, data = seedlings) n= 60, number of events= 60 coef exp(coef) se(coef) z Pr(>|z|) gapsize -0.001893 0.998109 0.593372 -0.003 0.997 gapsize:strata(cohort)cohort=September 0.717407 2.049112 0.860807 0.833 0.405 exp(coef) exp(-coef) lower .95 upper .95 gapsize 0.9981 1.002 0.3120 3.193 gapsize:strata(cohort)cohort=September 2.0491 0.488 0.3792 11.074 Rsquare= 0.022 (max possible= 0.993 ) Likelihood ratio test= 1.35 on 2 df, p=0.5097 Wald test = 1.32 on 2 df, p=0.5178 Score (logrank) test = 1.33 on 2 df, p=0.514 Anyone have an idea why this is occurring? Kind Regards Jacob
David Winsemius
2011-Jun-28 13:30 UTC
[R] coxph() - unexpected result using Crawley's seedlings data (The R Book)
On Jun 28, 2011, at 6:51 AM, Jacob Brogren wrote:> Hi, > > I ran the example on pp. 799-800 from Machael Crawley's "The R Book" > using package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The > model is a Cox's Proportional Hazards model. The result was quite > different compared to the R Book. I have compared my code to the > code in the book but can not find any differences in the function > call. My results are attached as well as a link to the results > presented in the book (link to Google Books).Shouldn't this instead go to Mr Crawley? Despite his unfortunately successful efforts to give the impression that his book is authoritative and somehow official, it is neither. I have never seen a contribution to the R effort from Mr. Crawley. -- David.> > When running the examples on pp. 797-799 I can't detect any > differences in results so I don't think there are errors in the data > set or in the creation of the status variable. > > --------------------------------- > Original from the R Book: > http://books.google.com/books?id=8D4HVx0apZQC&lpg=PA799&ots=rQgd_8ofeS&dq=r%20coxph%20crawley&pg=PA799#v > =onepage&q&f=false > > --------------------------------- > My result: >> summary(model1) > Call: > coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize, > data = seedlings) > > n= 60, number of events= 60 > > coef exp(coef) > se(coef) z Pr(>|z|) > gapsize -0.001893 0.998109 0.593372 > -0.003 0.997 > gapsize:strata(cohort)cohort=September 0.717407 2.049112 > 0.860807 0.833 0.405 > > exp(coef) exp(-coef) lower .95 > upper .95 > gapsize 0.9981 1.002 > 0.3120 3.193 > gapsize:strata(cohort)cohort=September 2.0491 0.488 > 0.3792 11.074 > > Rsquare= 0.022 (max possible= 0.993 ) > Likelihood ratio test= 1.35 on 2 df, p=0.5097 > Wald test = 1.32 on 2 df, p=0.5178 > Score (logrank) test = 1.33 on 2 df, p=0.514 > > Anyone have an idea why this is occurring? > > Kind Regards > > Jacob > ______________________________________________ > 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.David Winsemius, MD West Hartford, CT
Robert A LaBudde
2011-Jun-28 13:48 UTC
[R] coxph() - unexpected result using Crawley's seedlings data (The R Book)
Did you create the 'status' variable the way indicated on p. 797? Frequently with Surv() it pays to use syntax such as Surv(death, status==1) to make a clear logical statement of what is an event (status==1) vs. censored. PS. Next time include head(seedlings) and str(seedlings) to make clear what you are using as data. At 06:51 AM 6/28/2011, Jacob Brogren wrote:>Hi, > >I ran the example on pp. 799-800 from Machael Crawley's "The R Book" >using package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The >model is a Cox's Proportional Hazards model. The result was quite >different compared to the R Book. I have compared my code to the >code in the book but can not find any differences in the function >call. My results are attached as well as a link to the results >presented in the book (link to Google Books). > >When running the examples on pp. 797-799 I can't detect any >differences in results so I don't think there are errors in the data >set or in the creation of the status variable. > >--------------------------------- >Original from the R Book: >http://books.google.com/books?id=8D4HVx0apZQC&lpg=PA799&ots=rQgd_8ofeS&dq=r%20coxph%20crawley&pg=PA799#v=onepage&q&f=false > >--------------------------------- >My result: > > summary(model1) >Call: >coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize, > data = seedlings) > > n= 60, number of events= 60 > > coef > exp(coef) se(coef) z Pr(>|z|) >gapsize -0.001893 0.998109 0.593372 >-0.003 0.997 >gapsize:strata(cohort)cohort=September 0.717407 2.049112 0.860807 > 0.833 0.405 > > exp(coef) exp(-coef) lower > .95 upper .95 >gapsize 0.9981 1.002 >0.3120 3.193 >gapsize:strata(cohort)cohort=September 2.0491 0.488 >0.3792 11.074 > >Rsquare= 0.022 (max possible= 0.993 ) >Likelihood ratio test= 1.35 on 2 df, p=0.5097 >Wald test = 1.32 on 2 df, p=0.5178 >Score (logrank) test = 1.33 on 2 df, p=0.514 > >Anyone have an idea why this is occurring? > >Kind Regards > >Jacob >______________________________________________ >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.===============================================================Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: ral at lcfltd.com Least Cost Formulations, Ltd. URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239 Fax: 757-467-2947 "Vere scire est per causas scire"