Dear Thomas, The head of my dataset
> head(wsuv) parcel sp time censo treatment species 1 S8 Poecilanthe effusa ( Hub. ) Ducke. 1 1 1 1 2 S8 Poecilanthe effusa ( Hub. ) Ducke. 1 1 1 1 3 S8 Poecilanthe effusa ( Hub. ) Ducke. 1 1 1 1 4 S8 Poecilanthe effusa ( Hub. ) Ducke. 1 1 1 1 5 S8 Poecilanthe effusa ( Hub. ) Ducke. 1 1 1 1 6 S8 Poecilanthe effusa ( Hub. ) Ducke. 1 1 1 1 ... 144361 > summary(model.fit) # just one species from one treatment shown below Call: survfit(formula = Surv(time, censo) ~ treatment + species, data = wsuv) treatment=0, species=1 time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 15440 386 0.975 0.00126 0.973 0.977 2 15054 336 0.953 0.00170 0.950 0.957 3 14668 302 0.934 0.00200 0.930 0.938 4 14282 296 0.914 0.00226 0.910 0.919 5 13896 281 0.896 0.00247 0.891 0.901 6 13510 264 0.878 0.00264 0.873 0.883 7 13124 251 0.861 0.00280 0.856 0.867 8 12738 232 0.846 0.00293 0.840 0.852 9 12352 216 0.831 0.00305 0.825 0.837 10 11966 206 0.817 0.00315 0.811 0.823 11 11580 190 0.803 0.00325 0.797 0.810 12 11194 179 0.790 0.00333 0.784 0.797 13 10808 167 0.778 0.00341 0.772 0.785 14 10422 167 0.766 0.00349 0.759 0.773 15 10036 145 0.755 0.00356 0.748 0.762 16 9650 142 0.744 0.00363 0.737 0.751 17 9264 135 0.733 0.00369 0.726 0.740 18 8878 122 0.723 0.00375 0.715 0.730 19 8492 99 0.714 0.00380 0.707 0.722 20 8106 84 0.707 0.00385 0.699 0.714 21 7720 68 0.701 0.00389 0.693 0.708 22 7334 66 0.694 0.00393 0.687 0.702 23 6948 51 0.689 0.00397 0.681 0.697 24 6562 40 0.685 0.00400 0.677 0.693 25 6176 38 0.681 0.00403 0.673 0.689 26 5790 37 0.676 0.00407 0.669 0.684 27 5404 33 0.672 0.00411 0.664 0.680 28 5018 31 0.668 0.00415 0.660 0.676 29 4632 26 0.664 0.00419 0.656 0.673 30 4246 22 0.661 0.00423 0.653 0.669 31 3860 15 0.658 0.00427 0.650 0.667 32 3474 14 0.656 0.00431 0.647 0.664 33 3088 14 0.653 0.00436 0.644 0.661 34 2702 13 0.650 0.00443 0.641 0.658 35 2316 12 0.646 0.00451 0.638 0.655 36 1930 11 0.643 0.00462 0.634 0.652 37 1544 12 0.638 0.00480 0.628 0.647 38 1158 10 0.632 0.00507 0.622 0.642 39 772 9 0.625 0.00557 0.614 0.636 40 386 8 0.612 0.00709 0.598 0.626 I don't get why with 8 leaves remaining (out of 384), the survival is about 0.6??? Call: survfit(formula = Surv(time, censo) ~ 1, data = wsuv) n events median 0.95LCL 0.95UCL 144361 58830 40 39 40 > survfit(Surv(timee,ind)~sp2,data=wsuv) Call: survfit(formula = Surv(timee, ind) ~ sp2, data = wsuv) n events median 0.95LCL 0.95UCL sp2=1 32226 10856 Inf Inf Inf sp2=2 23370 9824 38 37 39 sp2=3 31201 13275 40 39 41 sp2=4 28044 10401 41 40 41 sp2=5 29520 14474 31 30 31 > survfit(Surv(timee,ind)~parcel2,data=wsuv) Call: survfit(formula = Surv(timee, ind) ~ parcel2, data = wsuv) n events median 0.95LCL 0.95UCL parcel2=0 68183 28116 38 38 38 parcel2=1 76178 30714 41 41 41 > survfit(Surv(timee,ind)~interaction(parcel2,sp2),data=wsuv) Call: survfit(formula = Surv(timee, ind) ~ interaction(parcel2, sp2), data = wsuv) n events median 0.95LCL 0.95UCL interaction(parcel2, sp2)=0.1 15826 5070 Inf Inf Inf interaction(parcel2, sp2)=1.1 16400 5786 Inf Inf Inf interaction(parcel2, sp2)=0.2 9430 3935 38 37 39 interaction(parcel2, sp2)=1.2 13940 5889 38 37 39 interaction(parcel2, sp2)=0.3 14678 6021 40 39 41 interaction(parcel2, sp2)=1.3 16523 7254 39 37 41 interaction(parcel2, sp2)=0.4 14473 5758 38 37 39 interaction(parcel2, sp2)=1.4 13571 4643 Inf Inf Inf interaction(parcel2, sp2)=0.5 13776 7332 29 28 30 interaction(parcel2, sp2)=1.5 15744 7142 33 32 34 Sorry, but I still cannot find what is going wrong. Regards, Paulo -----Mensagem original----- De: Thomas Lumley [mailto:[EMAIL PROTECTED] Enviada em: Wednesday, March 08, 2006 9:15 AM Para: Paulo Brando Cc: r-help@stat.math.ethz.ch Assunto: Re: [R] survival On Wed, 8 Mar 2006, Paulo Brando wrote: > Dear R-helpers, > > We marked 6000 leaves from 5 SPECIES - 10 individuals/species - in two > different TREATMENTs: a control and a dry-plot from which 50% of > incoming precipitation was excluded. We followed those leaves for 42 > months and noted the presence and absence at each visit. I then carried > out a Cox Harzard model to see differences in leaf mortality between > parcels and among species over time: > > leaves.cox <- coxph(Surv(time, censo) ~ treatment + species, data= wsuv) > > > When I plot 'survfitt(leaves.cox)', I come up with a survivor curve that > starts at 1 ends at 0.4. The problem is that at time 42 almost all > leaves are dead. I wander if surfit plot at time 42 should also be close > to zero? Yes, it probably should. It is a bit hard to tell what is going wrong, though. If you do plot(survfit(Surv(time,censo)~1,data=wsuv) plot(survfit(Surv(time,censo)~species,data=wsuv) plot(survfit(Surv(time,censo)~treatment,data=wsuv) plot(survfit(Surv(time,censo)~interaction(treatment,species),data=wsuv) you will get Kaplan-Meier estimates of survival curves. Looking at these may tell you what is happening. -thomas > > I followed examples from Venables and Ripley' book. (These analysis are > quite new for me). > >> summary(leaves.cox) > > Call: > coxph(formula = Surv(time, censo) ~ (treatment) + species, data = wsuv) > > n= 140840 > coef exp(coef) se(coef) z p > treatment -0.0209 0.98 0.00847 -2.47 0.014 > species 0.0712 1.07 0.00296 24.07 0.000 > > exp(coef) exp(-coef) lower .95 upper .95 > treatment 0.98 1.021 0.963 0.996 > species 1.07 0.931 1.068 1.080 > > Rsquare= 0.004 (max possible= 1 ) > Likelihood ratio test= 590 on 2 df, p=0 > Wald test = 587 on 2 df, p=0 > Score (logrank) test = 588 on 2 df, p=0 > > > My best regards and thanks in advance! > > Paulo > ________________________________________ > Paulo M. Brando > Instituto de Pesquisa Ambiental da Amazonia (IPAM) > Santarem, PA, Brasil. > Av. Rui Barbosa, 136. > Fone: + 55 93 3522 55 38 > www.ipam.org.br > E-mail: [EMAIL PROTECTED] ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html