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:tlumley at u.washington.edu]
Enviada em: Wednesday, March 08, 2006 9:15 AM
Para: Paulo Brando
Cc: r-help at 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, datawsuv)
>
>
> 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: pmbrando at ipam.org.br