Transplante de medula óssea – clássico e mudança no tempo (Inca)
Ester Rosa,Gabriel Thompson,Mateus Elias,Paulo Henrique
13 dezembro, 2020
1 ESCOLHA DO BANCO DE DADOS
1.1 Transplante de medula óssea – clássico e mudança no tempo (Inca)
Estes dados provêm de uma coorte de 96 pacientes submetidos a transplante de medula óssea para tratamento de leucemia mieloide crônica, no período de junho de 1986 a junho de 1998, no Centro de Transplante de Medula Óssea do Instituto Nacional do Câncer (Cemo – Inca). O acompanhamento dessa coorte possibilitou o estudo do efeito de fatores prognósticos para a ocorrência de doença do enxerto contra hospedeiro aguda e crônica, da sobrevivência livre de doença e da sobrevivência global. As covariáveis registradas para cada paciente estão descritas na tabela a seguir. Para mais detalhes acerca dessa coorte e do ajuste de modelo de sobrevivência de Cox clássico e Cox estendido, consulte Byington (1999).
Estes dados estão organizados em dois arquivos. O arquivo tmoclas.dat tem o formato clássico, com uma linha para cada paciente.
O arquivo tmopc.csv está preparado para uma análise de sobrevivência com covariáveis tempo-dependentes em que a variável que muda no tempo é a recuperação das plaquetas (reqplaq), doença do enxerto crônica (decr) e aguda (deag).
1.1.1 Variável: Descrição
id: identificador do paciente
sexo: 1 = masculino, 2 = feminino
idade: idade na data do transplante (5 a 53 anos)
status: 0 = censura, 1 = óbito
deag: doença enxerto aguda: 0 = não, 1 = sim
decr: doença enxerto crônica: 0 = não, 1 = sim
recplaq: 0 = plaquetas não recuperadas, 1 = recuperadas
fasegr: fase da doença na data do transplante agrupada: CP1 = 1a crônica e Other = outras
inicio: data do transplante ou da mudança de covariável
fim: data de mudança de covariável ou do fim do estudo
1.1.2 Dados
## Loading required package: rio##
## Attaching package: 'rio'## The following object is masked from 'package:plotly':
##
## export2 ESTIMADOR KAPLAN-MEIER
## Call: survfit(formula = Surv(t, censur) ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 31 194 1 0.995 0.00514 0.985 1.000
## 32 190 2 0.984 0.00895 0.967 1.000
## 39 174 1 0.979 0.01054 0.958 1.000
## 40 173 1 0.973 0.01190 0.950 0.997
## 48 164 1 0.967 0.01322 0.942 0.993
## 51 160 1 0.961 0.01446 0.933 0.990
## 54 159 3 0.943 0.01757 0.909 0.978
## 63 149 1 0.937 0.01856 0.901 0.974
## 65 145 1 0.930 0.01952 0.893 0.969
## 69 144 1 0.924 0.02043 0.885 0.965
## 70 143 1 0.917 0.02128 0.876 0.960
## 71 142 1 0.911 0.02209 0.869 0.955
## 74 140 1 0.904 0.02287 0.861 0.950
## 76 138 1 0.898 0.02362 0.853 0.945
## 78 135 1 0.891 0.02437 0.845 0.940
## 79 133 1 0.884 0.02509 0.837 0.935
## 83 130 1 0.878 0.02580 0.828 0.930
## 84 127 1 0.871 0.02651 0.820 0.924
## 98 123 1 0.864 0.02722 0.812 0.919
## 100 121 1 0.856 0.02791 0.803 0.913
## 101 120 1 0.849 0.02858 0.795 0.907
## 104 118 1 0.842 0.02923 0.787 0.901
## 120 110 1 0.834 0.02995 0.778 0.895
## 128 107 1 0.827 0.03067 0.769 0.889
## 139 105 1 0.819 0.03137 0.760 0.883
## 149 97 1 0.810 0.03216 0.750 0.876
## 162 92 1 0.802 0.03300 0.739 0.869
## 185 86 1 0.792 0.03390 0.728 0.862
## 200 82 1 0.783 0.03484 0.717 0.854
## 210 81 1 0.773 0.03572 0.706 0.846
## 211 80 1 0.763 0.03656 0.695 0.838
## 214 77 1 0.753 0.03741 0.683 0.830
## 216 76 1 0.743 0.03820 0.672 0.822
## 244 67 1 0.732 0.03921 0.659 0.813
## 261 64 1 0.721 0.04023 0.646 0.804
## 263 63 1 0.709 0.04119 0.633 0.795
## 281 59 1 0.697 0.04221 0.619 0.785
## 313 58 1 0.685 0.04316 0.606 0.775
## 347 55 1 0.673 0.04414 0.592 0.765
## 370 53 1 0.660 0.04510 0.577 0.755
## 371 52 1 0.648 0.04598 0.563 0.744
## 415 49 1 0.634 0.04690 0.549 0.733
## 425 48 1 0.621 0.04775 0.534 0.722
## 427 47 1 0.608 0.04853 0.520 0.711
## 434 46 1 0.595 0.04924 0.506 0.699
## 453 42 1 0.581 0.05006 0.490 0.687
## 475 40 1 0.566 0.05087 0.475 0.675
## 488 38 1 0.551 0.05167 0.459 0.662
## 522 36 1 0.536 0.05245 0.442 0.649
## 672 30 1 0.518 0.05366 0.423 0.6352.1 Gráfico KAPLAN-MEIER
3 ESTIMADOR DE NELSON-AALEN
## Call: survfit(formula = coxph(Surv(t, censur) ~ 1, method = "breslow"))
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 31 194 1 0.995 0.00513 0.985 1.000
## 32 190 2 0.984 0.00891 0.967 1.000
## 39 174 1 0.979 0.01050 0.958 1.000
## 40 173 1 0.973 0.01186 0.950 0.997
## 48 164 1 0.967 0.01318 0.942 0.993
## 51 160 1 0.961 0.01441 0.933 0.990
## 54 159 3 0.943 0.01748 0.910 0.978
## 63 149 1 0.937 0.01846 0.901 0.974
## 65 145 1 0.931 0.01943 0.893 0.969
## 69 144 1 0.924 0.02033 0.885 0.965
## 70 143 1 0.918 0.02119 0.877 0.960
## 71 142 1 0.911 0.02199 0.869 0.955
## 74 140 1 0.905 0.02277 0.861 0.950
## 76 138 1 0.898 0.02353 0.853 0.945
## 78 135 1 0.892 0.02427 0.845 0.940
## 79 133 1 0.885 0.02499 0.837 0.935
## 83 130 1 0.878 0.02570 0.829 0.930
## 84 127 1 0.871 0.02641 0.821 0.925
## 98 123 1 0.864 0.02712 0.813 0.919
## 100 121 1 0.857 0.02781 0.804 0.913
## 101 120 1 0.850 0.02848 0.796 0.908
## 104 118 1 0.843 0.02913 0.788 0.902
## 120 110 1 0.835 0.02984 0.779 0.896
## 128 107 1 0.827 0.03056 0.770 0.889
## 139 105 1 0.820 0.03126 0.760 0.883
## 149 97 1 0.811 0.03205 0.751 0.876
## 162 92 1 0.802 0.03288 0.740 0.869
## 185 86 1 0.793 0.03378 0.730 0.862
## 200 82 1 0.783 0.03472 0.718 0.855
## 210 81 1 0.774 0.03560 0.707 0.847
## 211 80 1 0.764 0.03643 0.696 0.839
## 214 77 1 0.754 0.03727 0.685 0.831
## 216 76 1 0.744 0.03806 0.674 0.823
## 244 67 1 0.733 0.03906 0.661 0.814
## 261 64 1 0.722 0.04008 0.648 0.805
## 263 63 1 0.711 0.04103 0.635 0.796
## 281 59 1 0.699 0.04204 0.621 0.786
## 313 58 1 0.687 0.04299 0.608 0.776
## 347 55 1 0.674 0.04396 0.594 0.766
## 370 53 1 0.662 0.04491 0.579 0.756
## 371 52 1 0.649 0.04579 0.565 0.745
## 415 49 1 0.636 0.04670 0.551 0.735
## 425 48 1 0.623 0.04755 0.536 0.724
## 427 47 1 0.610 0.04832 0.522 0.712
## 434 46 1 0.597 0.04903 0.508 0.701
## 453 42 1 0.583 0.04985 0.493 0.689
## 475 40 1 0.568 0.05065 0.477 0.677
## 488 38 1 0.554 0.05144 0.461 0.664
## 522 36 1 0.538 0.05222 0.445 0.651
## 672 30 1 0.521 0.05341 0.426 0.6373.1 Gráfico NELSON-AALEN
4 TEMPO MEDIANO
## Call: survfit(formula = Surv(t, censur) ~ 1)
##
## n events median 0.95LCL 0.95UCL
## 259 53 NA 453 NA5 TEMPO MÉDIO
## Call: survfit(formula = Surv(t, censur) ~ 1, conf.type = "plain")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 31 194 1 0.995 0.00514 0.985 1.000
## 32 190 2 0.984 0.00895 0.967 1.000
## 39 174 1 0.979 0.01054 0.958 0.999
## 40 173 1 0.973 0.01190 0.950 0.996
## 48 164 1 0.967 0.01322 0.941 0.993
## 51 160 1 0.961 0.01446 0.933 0.989
## 54 159 3 0.943 0.01757 0.909 0.977
## 63 149 1 0.937 0.01856 0.900 0.973
## 65 145 1 0.930 0.01952 0.892 0.968
## 69 144 1 0.924 0.02043 0.884 0.964
## 70 143 1 0.917 0.02128 0.876 0.959
## 71 142 1 0.911 0.02209 0.867 0.954
## 74 140 1 0.904 0.02287 0.859 0.949
## 76 138 1 0.898 0.02362 0.851 0.944
## 78 135 1 0.891 0.02437 0.843 0.939
## 79 133 1 0.884 0.02509 0.835 0.934
## 83 130 1 0.878 0.02580 0.827 0.928
## 84 127 1 0.871 0.02651 0.819 0.923
## 98 123 1 0.864 0.02722 0.810 0.917
## 100 121 1 0.856 0.02791 0.802 0.911
## 101 120 1 0.849 0.02858 0.793 0.905
## 104 118 1 0.842 0.02923 0.785 0.899
## 120 110 1 0.834 0.02995 0.776 0.893
## 128 107 1 0.827 0.03067 0.767 0.887
## 139 105 1 0.819 0.03137 0.757 0.880
## 149 97 1 0.810 0.03216 0.747 0.873
## 162 92 1 0.802 0.03300 0.737 0.866
## 185 86 1 0.792 0.03390 0.726 0.859
## 200 82 1 0.783 0.03484 0.714 0.851
## 210 81 1 0.773 0.03572 0.703 0.843
## 211 80 1 0.763 0.03656 0.692 0.835
## 214 77 1 0.753 0.03741 0.680 0.827
## 216 76 1 0.743 0.03820 0.669 0.818
## 244 67 1 0.732 0.03921 0.655 0.809
## 261 64 1 0.721 0.04023 0.642 0.800
## 263 63 1 0.709 0.04119 0.629 0.790
## 281 59 1 0.697 0.04221 0.615 0.780
## 313 58 1 0.685 0.04316 0.601 0.770
## 347 55 1 0.673 0.04414 0.586 0.759
## 370 53 1 0.660 0.04510 0.572 0.749
## 371 52 1 0.648 0.04598 0.557 0.738
## 415 49 1 0.634 0.04690 0.542 0.726
## 425 48 1 0.621 0.04775 0.528 0.715
## 427 47 1 0.608 0.04853 0.513 0.703
## 434 46 1 0.595 0.04924 0.498 0.691
## 453 42 1 0.581 0.05006 0.482 0.679
## 475 40 1 0.566 0.05087 0.466 0.666
## 488 38 1 0.551 0.05167 0.450 0.652
## 522 36 1 0.536 0.05245 0.433 0.639
## 672 30 1 0.518 0.05366 0.413 0.6235.1 Gráfico TEMPO MÉDIO
5.1.1 Resultado Tempo médio
## [1] 476.79866 CURVA TTT
##
## -- Column specification --------------------------------------------------------
## cols(
## id = col_double(),
## sexo = col_double(),
## idade = col_double(),
## status = col_double(),
## inicio = col_double(),
## fim = col_double(),
## deag = col_double(),
## decr = col_double(),
## recplaq = col_double(),
## fasegr = col_character()
## )## time censur Cum.Total.Time TTT i/n
## 68 31 1 7412 0.1501104 0.01886792
## 71 32 1 7602 0.1539583 0.03773585
## 73 32 1 7602 0.1539583 0.05660377
## 86 39 1 8865 0.1795370 0.07547170
## 87 40 1 9038 0.1830407 0.09433962
## 97 48 1 10373 0.2100776 0.11320755
## 100 51 1 10856 0.2198594 0.13207547
## 102 54 1 11333 0.2295198 0.15094340
## 103 54 1 11333 0.2295198 0.16981132
## 104 54 1 11333 0.2295198 0.18867925
## 111 63 1 12700 0.2572048 0.20754717
## 115 65 1 12992 0.2631185 0.22641509
## 116 69 1 13568 0.2747838 0.24528302
## 117 70 1 13711 0.2776799 0.26415094
## 118 71 1 13853 0.2805557 0.28301887
## 121 74 1 14275 0.2891022 0.30188679
## 122 76 1 14551 0.2946919 0.32075472
## 125 78 1 14823 0.3002005 0.33962264
## 127 79 1 14956 0.3028941 0.35849057
## 132 83 1 15477 0.3134455 0.37735849
## 133 84 1 15604 0.3160176 0.39622642
## 137 98 1 17338 0.3511351 0.41509434
## 139 100 1 17581 0.3560565 0.43396226
## 140 101 1 17701 0.3584867 0.45283019
## 143 104 1 18057 0.3656966 0.47169811
## 150 120 1 19845 0.4019078 0.49056604
## 153 128 1 20711 0.4194463 0.50943396
## 155 139 1 21872 0.4429593 0.52830189
## 163 149 1 22884 0.4634546 0.54716981
## 168 162 1 24095 0.4879802 0.56603774
## 174 185 1 26136 0.5293153 0.58490566
## 178 200 1 27388 0.5546712 0.60377358
## 179 210 1 28198 0.5710756 0.62264151
## 181 211 1 28278 0.5726958 0.64150943
## 183 214 1 28511 0.5774146 0.66037736
## 184 216 1 28663 0.5804929 0.67924528
## 194 244 1 30687 0.6214837 0.69811321
## 196 261 1 31779 0.6435992 0.71698113
## 197 263 1 31905 0.6461510 0.73584906
## 201 281 1 32991 0.6681451 0.75471698
## 202 313 1 34847 0.7057334 0.77358491
## 205 347 1 36754 0.7443547 0.79245283
## 207 370 1 37988 0.7693461 0.81132075
## 208 371 1 38040 0.7703992 0.83018868
## 211 415 1 40242 0.8149948 0.84905660
## 212 425 1 40722 0.8247160 0.86792453
## 213 427 1 40816 0.8266197 0.88679245
## 214 434 1 41138 0.8331409 0.90566038
## 218 453 1 41977 0.8501327 0.92452830
## 220 475 1 42872 0.8682585 0.94339623
## 222 488 1 43378 0.8785062 0.96226415
## 224 522 1 44603 0.9033153 0.98113208
## 230 672 1 49377 1.0000000 1.00000000
## Barlow-Proschan's test
## W=24.5105 k=53
## Z=-0.7155 p.value=0.2371
7 OBTENDO AJUSTES DOS MODELOS EXPONENCIAL
## $par
## [1] 964.381
##
## $value
## [1] 417.2415
##
## $counts
## function gradient
## 11 10
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1]
## [1,] 5.712764e-05## l AIC CAIC BIC
## [1,] 834.4829 836.4829 836.577 840.0398## parametros EP tvalue valorp LI LS
## [1,] 964.381 132.3052 7.289064 1.889554e-12 703.8502 1224.9128 OBTENDO AJUSTES DOS MODELOS WEIBULL
## $par
## [1] 929.0164139 0.7295653
##
## $value
## [1] 371.7855
##
## $counts
## function gradient
## 160 100
##
## $convergence
## [1] 1
##
## $message
## NULL
##
## $hessian
## [,1] [,2]
## [1,] -1.831779e-05 0.05475863
## [2,] 5.475863e-02 137.15171549## l AIC CAIC BIC
## [1,] 743.5709 747.5709 747.7284 754.6846## parametros EP tvalue valorp LI LS
## [1,] 929.0164139 NaN NaN NaN NaN NaN
## [2,] 0.7295653 0.05765385 12.65424 3.375429e-29 0.6160354 0.84309539 SOBREVIVÊNCIA ESTIMADA VERSUS KAPLEN-MEIER
##
## -- Column specification --------------------------------------------------------
## cols(
## id = col_double(),
## sexo = col_double(),
## idade = col_double(),
## status = col_double(),
## inicio = col_double(),
## fim = col_double(),
## deag = col_double(),
## decr = col_double(),
## recplaq = col_double(),
## fasegr = col_character()
## )## The following objects are masked from tmopc (pos = 3):
##
## deag, decr, fasegr, fim, id, idade, inicio, recplaq, sexo, status9.1 Gráfico SOBREVIVÊNCIA VERSUS KAPLEN-MEIER
#AJUSTE DE REGRESSÃO (PARAMÉTRICO)
options(warn=-1)
require(survival)
library(readr)
tmopc <- read_delim("tmopc.csv", ";", escape_double = FALSE,
trim_ws = TRUE)##
## -- Column specification --------------------------------------------------------
## cols(
## id = col_double(),
## sexo = col_double(),
## idade = col_double(),
## status = col_double(),
## inicio = col_double(),
## fim = col_double(),
## deag = col_double(),
## decr = col_double(),
## recplaq = col_double(),
## fasegr = col_character()
## )## The following objects are masked from tmopc (pos = 3):
##
## deag, decr, fasegr, fim, id, idade, inicio, recplaq, sexo, status## The following objects are masked from tmopc (pos = 4):
##
## deag, decr, fasegr, fim, id, idade, inicio, recplaq, sexo, statust=tmopc$fim
censur=tmopc$status
RE=tmopc$recplaq
CR=tmopc$decr
AG=tmopc$deag
y=log(t)
##função para obter os criterios
criterios=function(value,d,n){
#value -valor do optim
#d - numero de parametros
#n - numero da dados
l=2*value
AIC=l+2*d
BIC=l+d*log(n)
CAIC=AIC+(2*(d+2)*(d+3))/(n-d-3)
resul=cbind(l,AIC,CAIC,BIC)
return(resul)
}
IC=function(parametros,hessiana,n){
Var=solve(hessiana)
EP=diag(sqrt(Var))
tvalue=parametros/EP
LI=parametros-qt(0.975,n)*EP
LS=parametros+qt(0.975,n)*EP
valorp=pt(tvalue,n,lower.tail = FALSE)/2
resul=cbind(parametros,EP,tvalue,valorp,LI,LS)
return(resul)
}
##exponencial
verosi<-function(param){
b0<-param[1]
b1<-param[2]
b2<-param[3]
b3<-param[4]
mu=b0+b1*CR++b2*AG+b3*RE
z=(y-mu)
f=exp(z)*exp(-exp(z))
sb=exp(-exp(z))
lv<-censur*(log(f))+(1-censur)*log(sb);
sum(-lv)
}
valorin<-c(1,1,1,1)
estima1<-optim(valorin,verosi,method="BFGS",hessian=T);
estima1## $par
## [1] 6.124641 0.283440 -1.223466 1.437739
##
## $value
## [1] 138.7313
##
## $counts
## function gradient
## 100 34
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1] [,2] [,3] [,4]
## [1,] 53.00000 17.000038 27.000028 34.99997
## [2,] 17.00004 17.000038 7.236813 16.80226
## [3,] 27.00003 7.236813 27.000028 16.23747
## [4,] 34.99997 16.802256 16.237470 34.99997## l AIC CAIC BIC
## [1,] 277.4625 287.4625 287.9087 305.2467## parametros EP tvalue valorp LI LS
## [1,] 6.124641 0.2884220 21.2349965 3.024913e-59 5.5566901 6.6925916
## [2,] 0.283440 0.3352652 0.8454201 9.966417e-02 -0.3767528 0.9436327
## [3,] -1.223466 0.2775569 -4.4079828 4.999962e-01 -1.7700214 -0.6769105
## [4,] 1.437739 0.3309213 4.3446555 5.011429e-06 0.7861002 2.0893778t=tmopc$fim
censur=tmopc$status
RE=tmopc$recplaq
CR=tmopc$decr
AG=tmopc$deag
##função para obter os criterios
criterios=function(value,d,n){
#value -valor do optim
#d - numero de parametros
#n - numero da dados
l=2*value
AIC=l+2*d
BIC=l+d*log(n)
CAIC=AIC+(2*(d+2)*(d+3))/(n-d-3)
resul=cbind(l,AIC,CAIC,BIC)
return(resul)
}
IC=function(parametros,hessiana,n){
Var=solve(hessiana)
EP=diag(sqrt(Var))
tvalue=parametros/EP
LI=parametros-qt(0.975,n)*EP
LS=parametros+qt(0.975,n)*EP
valorp=pt(tvalue,n,lower.tail = FALSE)/2
resul=cbind(parametros,EP,tvalue,valorp,LI,LS)
return(resul)
}
##Weibull
verosi<-function(param){
sigma<-param[1]
b0<-param[2]
b1<-param[3]
b2<-param[4]
b3<-param[5]
if(any(sigma < 1e-20)) return(.Machine$double.xmax^.5)
mu=b0+b1*CR+b2*AG+b3*RE
z=(y-mu)/sigma
f=(1/sigma)*exp(z)*exp(-exp(z))
sb=exp(-exp(z))
lv<-censur*(log(f))+(1-censur)*log(sb);
sum(-lv)
}
valorin<-c(1,6.124641, 0.283440, -1.223466, 1.437739)
valorin## [1] 1.000000 6.124641 0.283440 -1.223466 1.437739## $par
## [1] 0.6719950 5.4754078 0.3458521 -0.8683463 1.6008134
##
## $value
## [1] 132.6471
##
## $counts
## function gradient
## 40 11
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1] [,2] [,3] [,4] [,5]
## [1,] 339.20492 -105.66050 -30.26532 -26.61026 -59.41548
## [2,] -105.66050 117.37004 37.64756 59.79576 77.51054
## [3,] -30.26532 37.64756 37.64756 14.67174 37.17204
## [4,] -26.61026 59.79576 14.67174 59.79576 33.53393
## [5,] -59.41548 77.51054 37.17204 33.53393 77.51054## l AIC CAIC BIC
## [1,] 265.2941 275.2941 275.7404 293.0783## parametros EP tvalue valorp LI LS
## [1,] 0.6719950 0.06921188 9.709244 8.275721e-20 0.53570538 0.8082847
## [2,] 5.4754078 0.23931967 22.879055 9.974494e-65 5.00414777 5.9466679
## [3,] 0.3458521 0.22526146 1.535336 3.148048e-02 -0.09772504 0.7894292
## [4,] -0.8683463 0.20320253 -4.273304 4.999932e-01 -1.26848570 -0.4682069
## [5,] 1.6008134 0.22824826 7.013475 5.051930e-12 1.15135476 2.050272010 LOG NORMAL
options(warn=-1)
##log-normal
verosi<-function(param){
sigma<-param[1]
b0<-param[2]
b1<-param[3]
b2<-param[4]
b3<-param[5]
if(any(sigma < 1e-20)) return(.Machine$double.xmax^.5)
mu=b0+b1*CR+b2*AG+b3*RE
z=(y-mu)/sigma
f=((1)/(sqrt(2*pi)*sigma))*(exp(-0.5*((z)^(2))))
sb=1-pnorm(z)
lv<-censur*(log(f))+(1-censur)*log(sb);
sum(-lv)
}
n=nrow(tmopc)
valorin<-c(0.6719950, 5.4754078, 0.3458521, -0.8683463, 1.6008134)
estima3<-optim(valorin,verosi,method="BFGS",hessian=T);
estima3## $par
## [1] 0.9891433 5.3193416 0.4812332 -0.8285952 1.3303994
##
## $value
## [1] 127.6028
##
## $counts
## function gradient
## 36 11
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1] [,2] [,3] [,4] [,5]
## [1,] 196.95978 -95.99080 -24.77648 -25.47159 -59.92532
## [2,] -95.99080 112.00697 34.67778 45.40716 75.21735
## [3,] -24.77648 34.67778 34.67778 10.79456 33.65570
## [4,] -25.47159 45.40716 10.79456 45.40716 27.78629
## [5,] -59.92532 75.21735 33.65570 27.78629 75.21735## l AIC CAIC BIC
## [1,] 255.2055 265.2055 265.6517 282.9897## parametros EP tvalue valorp LI LS
## [1,] 0.9891433 0.09726706 10.169356 2.924466e-21 0.79760840 1.1806783
## [2,] 5.3193416 0.22505117 23.636143 3.299173e-67 4.87617861 5.7625046
## [3,] 0.4812332 0.22786696 2.111904 8.913090e-03 0.03252547 0.9299410
## [4,] -0.8285952 0.20130920 -4.116032 4.999870e-01 -1.22500634 -0.4321841
## [5,] 1.3303994 0.22325440 5.959118 2.062731e-09 0.89077456 1.770024311 LOG LOGÍSTICA
options(warn=-1)
##log-logística
verosi<-function(param){
sigma<-param[1]
b0<-param[2]
b1<-param[3]
b2<-param[4]
b3<-param[5]
if(any(sigma < 1e-20)) return(.Machine$double.xmax^.5)
mu=b0+b1*CR+b2*AG+b3*RE
z=(y-mu)/sigma
f=(1/sigma)*exp(z)*(exp(1+exp(z)))^(-2)
sb=(1)/(1+exp(z))
lv<-censur*(log(f))+(1-censur)*log(sb);
sum(-lv)
}
n=nrow(tmopc)
valorin<-c(0.9891433, 5.3193416, 0.4812332, -0.8285952, 1.3303994)
estima4<-optim(valorin,verosi,method="BFGS",hessian=T);
estima4## $par
## [1] 0.5791330 5.3969886 0.5234551 -0.6807722 1.2732405
##
## $value
## [1] 241.5104
##
## $counts
## function gradient
## 31 12
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1] [,2] [,3] [,4] [,5]
## [1,] 455.65810 -145.59387 -37.65247 -56.05451 -84.81516
## [2,] -145.59387 129.04409 40.10253 69.86588 78.81756
## [3,] -37.65247 40.10253 40.10253 18.23727 39.25945
## [4,] -56.05451 69.86588 18.23727 69.86588 33.13828
## [5,] -84.81516 78.81756 39.25945 33.13828 78.81756## l AIC CAIC BIC
## [1,] 483.0207 493.0207 493.467 510.8049## parametros EP tvalue valorp LI LS
## [1,] 0.5791330 0.06087213 9.513928 3.350202e-19 0.45926570 0.6990003
## [2,] 5.3969886 0.22466927 24.021926 1.843410e-68 4.95457759 5.8393996
## [3,] 0.5234551 0.22135468 2.364780 4.694654e-03 0.08757105 0.9593391
## [4,] -0.6807722 0.19221107 -3.541795 4.998822e-01 -1.05926762 -0.3022768
## [5,] 1.2732405 0.21926871 5.806759 4.647011e-09 0.84146412 1.7050169## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "exponential")
##
## Coefficients:
## (Intercept) CR AG RE
## 6.1246408 0.2834457 -1.2234639 1.4377342
##
## Scale fixed at 1
##
## Loglik(model)= -399.7 Loglik(intercept only)= -423.4
## Chisq= 47.36 on 3 degrees of freedom, p= 2.91e-10
## n= 259##
## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "exponential")
## Value Std. Error z p
## (Intercept) 6.125 0.288 21.23 < 2e-16
## CR 0.283 0.335 0.85 0.4
## AG -1.223 0.278 -4.41 1.0e-05
## RE 1.438 0.331 4.34 1.4e-05
##
## Scale fixed at 1
##
## Exponential distribution
## Loglik(model)= -399.7 Loglik(intercept only)= -423.4
## Chisq= 47.36 on 3 degrees of freedom, p= 2.9e-10
## Number of Newton-Raphson Iterations: 6
## n= 259## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "weibull")
##
## Coefficients:
## (Intercept) CR AG RE
## 5.4752915 0.3458481 -0.8682290 1.6009021
##
## Scale= 0.6719733
##
## Loglik(model)= -393.6 Loglik(intercept only)= -423
## Chisq= 58.84 on 3 degrees of freedom, p= 1.04e-12
## n= 259##
## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "weibull")
## Value Std. Error z p
## (Intercept) 5.475 0.239 22.88 < 2e-16
## CR 0.346 0.225 1.54 0.12470
## AG -0.868 0.203 -4.27 1.9e-05
## RE 1.601 0.228 7.01 2.3e-12
## Log(scale) -0.398 0.103 -3.86 0.00011
##
## Scale= 0.672
##
## Weibull distribution
## Loglik(model)= -393.6 Loglik(intercept only)= -423
## Chisq= 58.84 on 3 degrees of freedom, p= 1e-12
## Number of Newton-Raphson Iterations: 6
## n= 259## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "lognorm")
##
## Coefficients:
## (Intercept) CR AG RE
## 5.3193412 0.4812343 -0.8285947 1.3303993
##
## Scale= 0.9891421
##
## Loglik(model)= -388.6 Loglik(intercept only)= -415.9
## Chisq= 54.58 on 3 degrees of freedom, p= 8.44e-12
## n= 259##
## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "lognorm")
## Value Std. Error z p
## (Intercept) 5.3193 0.2251 23.64 < 2e-16
## CR 0.4812 0.2279 2.11 0.035
## AG -0.8286 0.2013 -4.12 3.9e-05
## RE 1.3304 0.2233 5.96 2.5e-09
## Log(scale) -0.0109 0.0983 -0.11 0.912
##
## Scale= 0.989
##
## Log Normal distribution
## Loglik(model)= -388.6 Loglik(intercept only)= -415.9
## Chisq= 54.58 on 3 degrees of freedom, p= 8.4e-12
## Number of Newton-Raphson Iterations: 6
## n= 259## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "loglogistic")
##
## Coefficients:
## (Intercept) CR AG RE
## 5.2755960 0.4864628 -0.8453347 1.3643114
##
## Scale= 0.5571912
##
## Loglik(model)= -391.1 Loglik(intercept only)= -419.6
## Chisq= 57.1 on 3 degrees of freedom, p= 2.44e-12
## n= 259##
## Call:
## survreg(formula = Surv(t, censur) ~ CR + AG + RE, dist = "loglogistic")
## Value Std. Error z p
## (Intercept) 5.276 0.228 23.13 < 2e-16
## CR 0.486 0.234 2.08 0.038
## AG -0.845 0.204 -4.13 3.6e-05
## RE 1.364 0.231 5.90 3.6e-09
## Log(scale) -0.585 0.106 -5.52 3.3e-08
##
## Scale= 0.557
##
## Log logistic distribution
## Loglik(model)= -391.1 Loglik(intercept only)= -419.6
## Chisq= 57.1 on 3 degrees of freedom, p= 2.4e-12
## Number of Newton-Raphson Iterations: 5
## n= 259## [1] 1.61807## [1] 0.4366625## [1] 3.78255412 AJUSTE DE REGRESSÃO (COX)
##
## -- Column specification --------------------------------------------------------
## cols(
## id = col_double(),
## sexo = col_double(),
## idade = col_double(),
## status = col_double(),
## inicio = col_double(),
## fim = col_double(),
## deag = col_double(),
## decr = col_double(),
## recplaq = col_double(),
## fasegr = col_character()
## )## The following objects are masked from tmopc (pos = 3):
##
## deag, decr, fasegr, fim, id, idade, inicio, recplaq, sexo, status## The following objects are masked from tmopc (pos = 4):
##
## deag, decr, fasegr, fim, id, idade, inicio, recplaq, sexo, status## The following objects are masked from tmopc (pos = 5):
##
## deag, decr, fasegr, fim, id, idade, inicio, recplaq, sexo, status## Call:
## coxph(formula = Surv(t, censur) ~ CR + AG + RE, data = tmopc,
## x = T, method = "breslow")
##
## n= 259, number of events= 53
##
## coef exp(coef) se(coef) z Pr(>|z|)
## CR -0.43153 0.64951 0.33928 -1.272 0.203411
## AG 0.97292 2.64566 0.28563 3.406 0.000659 ***
## RE -2.59944 0.07432 0.46353 -5.608 2.05e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## CR 0.64951 1.540 0.33403 1.2629
## AG 2.64566 0.378 1.51150 4.6308
## RE 0.07432 13.456 0.02996 0.1843
##
## Concordance= 0.787 (se = 0.027 )
## Likelihood ratio test= 56.88 on 3 df, p=3e-12
## Wald test = 52.06 on 3 df, p=3e-11
## Score (logrank) test = 80.92 on 3 df, p=<2e-16## [1] -240.3676 -211.9269
13 TRV
#HO= EXPONENCIAL MODELO MAIS ADEQUADO
#H1= EXPONENCIAL NÃO É O MAIS ADEQUADO
#TRV
trv=834.4829-743.5709
trv## [1] 90.912## [1] 014 CONCLUSÃO
HO= EXPONENCIAL MODELO MAIS ADEQUADO
H1= EXPONENCIAL NÃO É O MAIS ADEQUADO
Portanto o como foi aceito H0 o modelo exponencial é mais adequado para os dados
15 CRÉDITOS
Criado por :
Ester Rosa
Gabriel Thompson
Mateus Elias
Paulo Henrique