Analisis Supervivencia PBC
Omar Escobar
26 de agosto de 2019
Curva de Kaplan Meier
Lo primero que vamos a hacer es cargar la libreria que vamos a utilizar para el analisis de supervivencia. y recomiendo para graficar con ggplot la libreria “survminer”, al igual de “muhaz” para utilizar una herramienta que permite calcular la funcion Hazard.
library(survival)
library(KMsurv)
library(survMisc)
library(survminer)
library(ggfortify)
library(flexsurv)
library(actuar)
library(dplyr)
library(survminer)
library(muhaz)## Warning: package 'survMisc' was built under R version 3.5.3## Warning: package 'survminer' was built under R version 3.5.3## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following object is masked from 'package:survMisc':
##
## autoplot## Loading required package: ggpubr## Loading required package: magrittr## Warning: package 'ggfortify' was built under R version 3.5.3## Warning: package 'flexsurv' was built under R version 3.5.3## Warning: package 'actuar' was built under R version 3.5.3##
## Attaching package: 'actuar'## The following objects are masked from 'package:flexsurv':
##
## dllogis, pllogis, qllogis, rllogis## The following object is masked from 'package:grDevices':
##
## cm## Warning: package 'dplyr' was built under R version 3.5.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: package 'muhaz' was built under R version 3.5.3Agregar la base de datos de Colangitis biliar primaria o PBC
antes de agregar la base de datos se puede modificar en excel o en SPSS, dado que hay algunas variables que se deben convertir como la albumina, protrombina y la bilirrubina en su escala logaritmica.
Luego de agregar la base de datos
library(haven)
PBC <- read_sav("~/Maestria Epidemiologia/Tercer Semestre/Bioestadistica III/Trabajo Final/PBC datos_2.sav")
View(PBC)AHORA SI A LO QUE VINIMOS
Se tiene que crear un Objet surv de la siguiente manera
## Add survival object. status == 1 is death
PBC.surv <- Surv(PBC$years, PBC$status)Luego creaos los estimadoes de Kaplan-Meier
PBC.km <- survfit(PBC.surv ~ 1, data = PBC, type = "kaplan-meier")
head(PBC.km)## Warning in if (j > 1) stop("subscript out of bounds"): la condición tiene
## longitud > 1 y sólo el primer elemento será usado## Call: survfit(formula = PBC.surv ~ 1, data = PBC, type = "kaplan-meier")
##
## n events median 0.95LCL 0.95UCL
## 312.00 125.00 9.30 8.45 10.55Y procedemos a realizar la grafica
ggsurvplot(fit = PBC.km, data = PBC, conf.int = T, title = "Curva de Supervivencia",
xlab = "Tiempo", ylab = "Probabilidad de supervivencia", legend.title = "Estimación",
legend.labs = "Kaplan-Meier")
o de la siguiente manera
plot(PBC.surv, xlab="Survival Time in years",
ylab="% Surviving", yscale=100,
main="Survival Distribution (Overall)") 
Estimacion de la funcion de riesgo acumulado
maxima verosimilitud
Hdet <- PBC.km %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))graficar el riesgo acumulado
plot(PBC.km, fun = "cumhaz", conf.int = F, main = "Riesgo Acumulado", col = 1:2,
xlab = "Tiempo", ylab = "Riesgo Acumulado")
Con covariables
La forma de hacerlo por estratos, yo lo hari ade la siguiente manera, ceando un vector para el objeto de supervivencia distribuido por la variable histologia, posterior agregar la columna riesgo acumulado
km.as.hist <- survfit(PBC.surv ~ histol, data = PBC)
Histol <- km.as.hist %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))graficar el riesgo acumulado
plot(km.as.hist, fun = "cumhaz", conf.int = F, main = "Riesgo Acumulado por Histologia", col = 1:2,
xlab = "Tiempo", ylab = "Riesgo Acumulado")
ggsurvplot(km.as.hist, fun = "cumhaz", xlab = "Tiempo", censor = T,
ylab = "Riesgo Acumulado", title = "Riesgo Acumulado por Histologia", legend.title = "Grado histologico")
Estimacion de Nelson Aalen
qplot(time, CumHaz, data = Hdet, geom = "step", xlab = "Tiempo (years)", ylab = "Riesgo Acumulado",
main = "Riesgo Acumulado")
qplot(time, CumHaz, col = strata, data = R, geom = "step", xlab = "Tiempo (Semanas)",
ylab = "Riesgo Acumulado", main = "Riesgo Acumulado")logRank
Despues vamos a crear la grafica de logRank para comparacion por dos variables, en el caso de este estudio es por el tratamiento (rx), se debe crear datos donde se separan por tratamiento.
al igual que podemos crear de una vez datos para comparar por sexo, que es otra de las variables dicotomicas.
km.by.sex <- survfit(PBC.surv ~ sex, data = PBC, conf.type = "log-log")
km.by.rx <- survfit(PBC.surv ~ rx, data = PBC, conf.type = "log-log")
## y el siguiente codigo permite crear un datos mas facil para posterior graficar por medio de Log-Log
km.as.one <- survfit(PBC.surv ~ 1, data = PBC, conf.type = "log-log")
## Show object
km.as.one## Call: survfit(formula = PBC.surv ~ 1, data = PBC, conf.type = "log-log")
##
## n events median 0.95LCL 0.95UCL
## 312.00 125.00 9.30 8.45 10.51km.by.sex## Call: survfit(formula = PBC.surv ~ sex, data = PBC, conf.type = "log-log")
##
## n events median 0.95LCL 0.95UCL
## sex=0 36 22 6.53 3.15 11.2
## sex=1 276 103 9.39 8.46 10.5km.by.rx## Call: survfit(formula = PBC.surv ~ rx, data = PBC, conf.type = "log-log")
##
## n events median 0.95LCL 0.95UCL
## rx=0 158 65 8.99 6.95 11.5
## rx=1 154 60 9.39 8.46 10.5Con este metodo podemos observar el numero de eventos de cada grupo su mediana y sus intervalos de confianza.
y con el siguiente se puede observar todos los datos.
## See survival estimates at given time (lots of outputs)
summary.data.frame (km.by.rx, maxsum = 10)## n time n.risk n.event
## Min. :154 Min. : 0.1123 Min. : 1.00 Min. :0.0000
## 1st Qu.:155 1st Qu.: 3.2752 1st Qu.: 39.25 1st Qu.:0.0000
## Median :156 Median : 5.1417 Median : 77.50 Median :0.0000
## Mean :156 Mean : 5.5165 Mean : 78.15 Mean :0.4085
## 3rd Qu.:157 3rd Qu.: 7.4134 3rd Qu.:117.00 3rd Qu.:1.0000
## Max. :158 Max. :12.4736 Max. :158.00 Max. :2.0000
## n.censor surv type strata
## Min. :0.0000 Min. :0.3186 Length:1 Min. :151
## 1st Qu.:0.0000 1st Qu.:0.5866 Class :character 1st Qu.:152
## Median :1.0000 Median :0.7052 Mode :character Median :153
## Mean :0.6111 Mean :0.6886 Mean :153
## 3rd Qu.:1.0000 3rd Qu.:0.7920 3rd Qu.:154
## Max. :2.0000 Max. :0.9937 Max. :155
## std.err lower upper conf.type
## Min. :0.006349 Min. :0.1735 Min. :0.4737 Length:1
## 1st Qu.:0.041406 1st Qu.:0.4900 1st Qu.:0.6779 Class :character
## Median :0.054304 Median :0.6227 Median :0.7730 Mode :character
## Mean :0.068981 Mean :0.6030 Mean :0.7602
## 3rd Qu.:0.077602 3rd Qu.:0.7187 3rd Qu.:0.8481
## Max. :0.248664 Max. :0.9559 Max. :0.9991
## conf.int call
## Min. :0.95 Length:4
## 1st Qu.:0.95 Class :call
## Median :0.95 Mode :call
## Mean :0.95
## 3rd Qu.:0.95
## Max. :0.95La forma de graficar mas sencillo es la siguiente
## Plotting without any specification
plot(km.as.one)
## Plotting without any specification
plot(km.by.sex)
## Plotting without any specification
plot(km.by.rx)
sin embargo tambien se pueden modificar los tirulos y marcas, hasta colores
plot(km.by.rx, xlab="Survival Time in Years",
ylab="% Surviving", yscale=100, col=c("red","blue"),
main="Survival Distributions by Treatment")
legend("topright", title="Tratamiento", c("Placebo", "DPCA"),
fill=c("red", "blue"))
Y para crear un grafico de riesgo acumlado por grupos
ggsurvplot(km.by.rx, fun = "event", censor = F, xlab = "Time (years)")
Estimacion de Mediana
Prueba de hipótesis para igualdad de dos o más funciones de supervivencia.
Y para el metodo estadistico para comparacion de curvas es el siguiente
Es necesario tener el objeto surv que habiamos creado al inicio
##### test for difference between placebo y DPCA
# survival curves (logrank test)
survdiff(PBC.surv~rx, data=PBC) ## Call:
## survdiff(formula = PBC.surv ~ rx, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=0 158 65 63.2 0.0502 0.102
## rx=1 154 60 61.8 0.0513 0.102
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7# survival curves (Peto&Peto)
survdiff(PBC.surv~rx, data=PBC, rho = 1) ## Call:
## survdiff(formula = PBC.surv ~ rx, data = PBC, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=0 158 49.7 49.0 0.00963 0.0243
## rx=1 154 46.8 47.5 0.00993 0.0243
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9Al igual se puede hacer con las demas variables para ver cuales on significativas y cuales se deben ingresar al modelo.
## Call:
## survdiff(formula = PBC.surv ~ sex, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=0 36 22 14.6 3.728 4.27
## sex=1 276 103 110.4 0.494 4.27
##
## Chisq= 4.3 on 1 degrees of freedom, p= 0.04## Call:
## survdiff(formula = PBC.surv ~ rx, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=0 158 65 63.2 0.0502 0.102
## rx=1 154 60 61.8 0.0513 0.102
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7## Call:
## survdiff(formula = PBC.surv ~ edema, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## edema=0 263 89 113.4 5.24 56.8
## edema=1 49 36 11.6 51.14 56.8
##
## Chisq= 56.8 on 1 degrees of freedom, p= 5e-14## Call:
## survdiff(formula = PBC.surv ~ asictes, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## asictes=0 288 102 121.28 3.06 104
## asictes=1 24 23 3.72 99.91 104
##
## Chisq= 104 on 1 degrees of freedom, p= <2e-16## Call:
## survdiff(formula = PBC.surv ~ spiders, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## spiders=0 222 73 98.1 6.43 30.3
## spiders=1 90 52 26.9 23.44 30.3
##
## Chisq= 30.3 on 1 degrees of freedom, p= 4e-08## Call:
## survdiff(formula = PBC.surv ~ histol, data = PBC)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## histol=1 16 1 9.9 8.00 8.78
## histol=2 67 16 32.3 8.22 11.17
## histol=3 120 43 51.2 1.30 2.22
## histol=4 109 65 31.6 35.17 47.84
##
## Chisq= 53.8 on 3 degrees of freedom, p= 1e-11Modelo Predictivo de COX
Se deben corroborar que las variables nominales sean un factor, de no ser asi toca convertilo con la funcion factor
PBC$sex <- factor(x = PBC$sex, levels = c(0, 1), labels = c("H", "M"))PBC$asictes <- factor(x = PBC$asictes, levels = c(0, 1), labels = c("Ausente", "Presente"))
PBC$hepatom <- factor(x = PBC$hepatom, levels = c(0, 1), labels = c("Ausente", "Presente"))
PBC$edema <- factor(x = PBC$edema, levels = c(0, 1), labels = c("Ausente", "Presente"))
PBC$spiders <- factor(x = PBC$spiders, levels = c(0, 1), labels = c("Ausente", "Presente"))Crear el modelo
Para crear el modelo se seleccionas la variables con significancia
ModCoxPBC <- coxph(PBC.surv ~ age + logbili + logalbu + logprot + copper , data = PBC)la variable edema la considere una posible variable confusora por lo que cree el siguiente modelo para realizar su comparacion.
ModCoxPBCedema <- coxph(PBC.surv ~ age + logbili + logalbu + logprot + edema , data = PBC)display results
ModCoxPBC## Call:
## coxph(formula = PBC.surv ~ age + logbili + logalbu + logprot +
## copper, data = PBC)
##
## coef exp(coef) se(coef) z p
## age 0.028833 1.029253 0.008567 3.366 0.000764
## logbili 0.843135 2.323641 0.108050 7.803 6.04e-15
## logalbu -3.407564 0.033122 0.687241 -4.958 7.11e-07
## logprot 3.576517 35.748802 0.997247 3.586 0.000335
## copper 0.002271 1.002273 0.001002 2.267 0.023390
##
## Likelihood ratio test=197.7 on 5 df, p=< 2.2e-16
## n= 310, number of events= 124
## (2 observations deleted due to missingness)ModCoxPBCedema## Call:
## coxph(formula = PBC.surv ~ age + logbili + logalbu + logprot +
## edema, data = PBC)
##
## coef exp(coef) se(coef) z p
## age 0.032602 1.033140 0.008626 3.779 0.000157
## logbili 0.894493 2.446095 0.098032 9.124 < 2e-16
## logalbu -3.228639 0.039611 0.707401 -4.564 5.02e-06
## logprot 3.235409 25.416777 0.998202 3.241 0.001190
## edema 0.401057 1.493402 0.221426 1.811 0.070103
##
## Likelihood ratio test=195.9 on 5 df, p=< 2.2e-16
## n= 312, number of events= 125validacion de supuestos
evaluar: proportional hazards assumption
cada covariable calcula el coeficiente rho, el estadístico de prueba chisq y su respectivo p-valor, además incluye una prueba global
y otro metodo mas grafico de evaluar la proporcionalidad de riesgos es la siguiente funcion que se encuentra en el paquete servminer
library(survminer)
COX.PBC<- cox.zph(ModCoxPBC)
COX.PBC.edema <- cox.zph(ModCoxPBCedema)
COX.PBC## rho chisq p
## age 0.00035 1.39e-05 0.9970
## logbili 0.13844 2.31e+00 0.1285
## logalbu 0.04341 2.37e-01 0.6262
## logprot -0.24080 5.75e+00 0.0165
## copper -0.08131 8.14e-01 0.3669
## GLOBAL NA 7.73e+00 0.1718COX.PBC.edema ## rho chisq p
## age -0.04078 0.17921 0.6721
## logbili 0.15785 2.80313 0.0941
## logalbu -0.00786 0.00797 0.9288
## logprot -0.20503 4.05904 0.0439
## edema -0.20243 4.53766 0.0332
## GLOBAL NA 12.34570 0.0303l valor-p no es menor a .05 por lo que no podemos rechazar la hipótesis de riesgos proporcionales.
la forma grafica de realizarlo es
ggcoxzph(COX.PBC)
otro comando para realizar esot seria
ggcoxzph(cox.zph(ModCoxPBC))Y obtenemos la misma grafica con el test de Shchoenfeld
Residuales
ggcoxdiagnostics(ModCoxPBC, type = "schoenfeld")
ggcoxdiagnostics(ModCoxPBC, type = "dfbetas", ox.scale = "observation.id", title = "Residuales dfbetas",
subtitle = "dfbetas vs id", caption = "Análisis residual")
linealidad
ggcoxdiagnostics(ModCoxPBC, type = "martingale")
Prediccion
Numero 1
datos1<- data.frame(age = c(52), logbili = c(log(0.5)), logalbu = c(log(4.5)), logprot = c(log(10.1)), copper= c(98))
Pred1 <- survfit(ModCoxPBC, newdata = datos1, conf.type = "log-log")
ggsurvplot(fit = Pred1, data = datos1, conf.int = T, title = "Curva de Supervivencia",
xlab = "Tiempo", ylab = "Probabilidad de supervivencia")
0.027*52 + 0.826*log(0.5) + 3.323*log(10.1) - 3.177*log(4.5) + 0.002*98## [1] 3.933562Numero 2
datos2 <- data.frame(age = c(52), logbili = c(log(13.9)), logalbu = c(log(2.8)), logprot = c(log(13.8)), copper= c(98))
Pred2 <- survfit(ModCoxPBC, data = datos2, conf.type = "log-log")
ggsurvplot(fit = Pred2, data = datos2, conf.int = T, title = "Curva de Supervivencia",
xlab = "Tiempo", ylab = "Probabilidad de supervivencia")
0.027*52 + 0.826*log(13.9) + 3.323*log(13.8) - 3.177*log(2.8) + 0.002*98## [1] 9.224613Estimacin de riesgo de Base
La estimación de riesgo base,predicción de la curva de supervivencia utilizando las covariables igual a 0
Rbase <- data.frame(age = c(0), logbili = c(log(0)), logalbu = c(log(0)), logprot = c(log(0)), copper= c(0))
Rbase <- survfit(ModCoxPBC, data = Rbase, conf.type = "log-log")
cbind(Tiempo = Rbase$time, RiesgoAcumulado = Rbase$cumhaz)## Tiempo RiesgoAcumulado
## [1,] 0.1122519 0.0008263517
## [2,] 0.1396304 0.0017807378
## [3,] 0.1943874 0.0027640869
## [4,] 0.2108145 0.0037585667
## [5,] 0.3011636 0.0047757565
## [6,] 0.3559206 0.0057941287
## [7,] 0.3586585 0.0068410770
## [8,] 0.3832991 0.0079570356
## [9,] 0.4900753 0.0091028593
## [10,] 0.5092403 0.0103001102
## [11,] 0.5229295 0.0115035537
## [12,] 0.5420945 0.0127127834
## [13,] 0.5667351 0.0139223069
## [14,] 0.5913758 0.0151813549
## [15,] 0.6105407 0.0165332335
## [16,] 0.7227926 0.0193504587
## [17,] 0.8323066 0.0208647755
## [18,] 0.8788501 0.0223915688
## [19,] 0.8925394 0.0239312822
## [20,] 0.9144422 0.0254863291
## [21,] 0.9527721 0.0271229393
## [22,] 1.0622861 0.0287809750
## [23,] 1.0951403 0.0304489169
## [24,] 1.2594113 0.0322293255
## [25,] 1.4099932 0.0340201108
## [26,] 1.4592744 0.0340201108
## [27,] 1.5030801 0.0358152255
## [28,] 1.5112936 0.0376766522
## [29,] 1.6344969 0.0395508763
## [30,] 1.6728269 0.0414428872
## [31,] 1.8425735 0.0433458848
## [32,] 1.9000684 0.0452693544
## [33,] 1.9383984 0.0471977901
## [34,] 2.0041068 0.0471977901
## [35,] 2.0068445 0.0491412256
## [36,] 2.0177960 0.0491412256
## [37,] 2.0533881 0.0511160818
## [38,] 2.0862422 0.0531093149
## [39,] 2.1054072 0.0551103149
## [40,] 2.1519508 0.0571206653
## [41,] 2.1574264 0.0571206653
## [42,] 2.1629021 0.0591434418
## [43,] 2.1820672 0.0611815282
## [44,] 2.1875427 0.0632624541
## [45,] 2.2915812 0.0632624541
## [46,] 2.2970569 0.0632624541
## [47,] 2.3271732 0.0653632929
## [48,] 2.3353868 0.0674780704
## [49,] 2.3518138 0.0696894631
## [50,] 2.4010952 0.0696894631
## [51,] 2.4366872 0.0719800995
## [52,] 2.4668036 0.0719800995
## [53,] 2.4750171 0.0743187585
## [54,] 2.5462012 0.0766717706
## [55,] 2.5708418 0.0766717706
## [56,] 2.5817933 0.0790923870
## [57,] 2.6584532 0.0815874049
## [58,] 2.6666667 0.0841122343
## [59,] 2.6830938 0.0866633220
## [60,] 2.7214236 0.0866633220
## [61,] 2.7351129 0.0892384338
## [62,] 2.7378509 0.0918238805
## [63,] 2.7707050 0.0944346952
## [64,] 2.8199863 0.0944346952
## [65,] 2.8391514 0.0970682567
## [66,] 2.9212868 0.0970682567
## [67,] 2.9486654 0.0998271937
## [68,] 2.9568789 0.1026065711
## [69,] 2.9650924 0.1053991054
## [70,] 2.9678302 0.1053991054
## [71,] 3.1457906 0.1053991054
## [72,] 3.1540041 0.1082475770
## [73,] 3.1567419 0.1082475770
## [74,] 3.1895962 0.1111434722
## [75,] 3.2032855 0.1141840981
## [76,] 3.2607803 0.1207502718
## [77,] 3.3182752 0.1243409527
## [78,] 3.3292265 0.1243409527
## [79,] 3.3319645 0.1279530405
## [80,] 3.3675566 0.1279530405
## [81,] 3.3785079 0.1279530405
## [82,] 3.3812456 0.1316500467
## [83,] 3.4223135 0.1316500467
## [84,] 3.4798083 0.1316500467
## [85,] 3.5400410 0.1316500467
## [86,] 3.5455167 0.1316500467
## [87,] 3.5509925 0.1354305969
## [88,] 3.5592060 0.1354305969
## [89,] 3.5619438 0.1354305969
## [90,] 3.5646818 0.1354305969
## [91,] 3.6139631 0.1354305969
## [92,] 3.6167009 0.1354305969
## [93,] 3.6386037 0.1354305969
## [94,] 3.6933608 0.1354305969
## [95,] 3.6960986 0.1397219084
## [96,] 3.7125256 0.1440761587
## [97,] 3.7234771 0.1485686846
## [98,] 3.7316906 0.1485686846
## [99,] 3.8357289 0.1485686846
## [100,] 3.8548939 0.1485686846
## [101,] 3.8658452 0.1485686846
## [102,] 3.8685832 0.1531605401
## [103,] 3.8822725 0.1531605401
## [104,] 3.8877480 0.1531605401
## [105,] 3.9069130 0.1580148804
## [106,] 3.9233401 0.1580148804
## [107,] 3.9260781 0.1629322870
## [108,] 3.9288158 0.1629322870
## [109,] 3.9534566 0.1678898292
## [110,] 3.9616702 0.1678898292
## [111,] 3.9835730 0.1678898292
## [112,] 3.9890485 0.1678898292
## [113,] 4.0547571 0.1678898292
## [114,] 4.0711842 0.1729343368
## [115,] 4.0848732 0.1780213718
## [116,] 4.1177278 0.1780213718
## [117,] 4.1752224 0.1780213718
## [118,] 4.2053390 0.1831772363
## [119,] 4.2217660 0.1831772363
## [120,] 4.2655716 0.1831772363
## [121,] 4.2929502 0.1831772363
## [122,] 4.2956877 0.1831772363
## [123,] 4.3148527 0.1884921581
## [124,] 4.3586583 0.1884921581
## [125,] 4.4188910 0.1884921581
## [126,] 4.4216290 0.1884921581
## [127,] 4.5338807 0.1884921581
## [128,] 4.5366187 0.1942118250
## [129,] 4.5612593 0.1942118250
## [130,] 4.5913758 0.1942118250
## [131,] 4.6050649 0.2001501887
## [132,] 4.6269679 0.2121734777
## [133,] 4.6570840 0.2121734777
## [134,] 4.6598220 0.2121734777
## [135,] 4.7501712 0.2121734777
## [136,] 4.7665982 0.2183761970
## [137,] 4.8323069 0.2183761970
## [138,] 4.8432579 0.2183761970
## [139,] 4.8459959 0.2183761970
## [140,] 4.8624229 0.2183761970
## [141,] 4.8815880 0.2183761970
## [142,] 4.8870635 0.2183761970
## [143,] 4.8898015 0.2260459551
## [144,] 4.9007530 0.2260459551
## [145,] 4.9555101 0.2260459551
## [146,] 5.0020533 0.2338930805
## [147,] 5.0130048 0.2338930805
## [148,] 5.0157428 0.2338930805
## [149,] 5.0568104 0.2420360642
## [150,] 5.1526351 0.2420360642
## [151,] 5.2238193 0.2420360642
## [152,] 5.2703629 0.2505180837
## [153,] 5.2895279 0.2505180837
## [154,] 5.3251200 0.2505180837
## [155,] 5.3415470 0.2505180837
## [156,] 5.3853526 0.2505180837
## [157,] 5.4154687 0.2505180837
## [158,] 5.4182067 0.2505180837
## [159,] 5.5359344 0.2505180837
## [160,] 5.5660505 0.2505180837
## [161,] 5.6125941 0.2505180837
## [162,] 5.6262832 0.2596597803
## [163,] 5.6974673 0.2688092272
## [164,] 5.7221084 0.2780515880
## [165,] 5.7631760 0.2873257188
## [166,] 5.7659140 0.2873257188
## [167,] 5.9055443 0.2873257188
## [168,] 5.9356604 0.2873257188
## [169,] 5.9411364 0.2873257188
## [170,] 5.9438739 0.2873257188
## [171,] 5.9575634 0.2873257188
## [172,] 5.9630389 0.2873257188
## [173,] 6.0095825 0.2873257188
## [174,] 6.0670772 0.2873257188
## [175,] 6.0807667 0.2873257188
## [176,] 6.0889802 0.2969998477
## [177,] 6.1355238 0.2969998477
## [178,] 6.1738534 0.2969998477
## [179,] 6.1765914 0.3067943219
## [180,] 6.2203970 0.3067943219
## [181,] 6.2642026 0.3177098150
## [182,] 6.2806296 0.3177098150
## [183,] 6.2888432 0.3288700257
## [184,] 6.2997947 0.3288700257
## [185,] 6.3463383 0.3288700257
## [186,] 6.3791924 0.3288700257
## [187,] 6.3846679 0.3288700257
## [188,] 6.4339495 0.3288700257
## [189,] 6.4531145 0.3288700257
## [190,] 6.4695415 0.3288700257
## [191,] 6.4750171 0.3288700257
## [192,] 6.5325122 0.3409355565
## [193,] 6.5708418 0.3533906733
## [194,] 6.6228609 0.3661897506
## [195,] 6.6885695 0.3661897506
## [196,] 6.7049966 0.3661897506
## [197,] 6.7132101 0.3661897506
## [198,] 6.7241616 0.3661897506
## [199,] 6.7515402 0.3792760929
## [200,] 6.7570157 0.3792760929
## [201,] 6.7761807 0.3792760929
## [202,] 6.8528404 0.3927180888
## [203,] 6.8555784 0.3927180888
## [204,] 6.9185491 0.3927180888
## [205,] 6.9541411 0.4064951457
## [206,] 6.9952087 0.4064951457
## [207,] 6.9979467 0.4064951457
## [208,] 7.0171118 0.4064951457
## [209,] 7.0444899 0.4064951457
## [210,] 7.0472279 0.4064951457
## [211,] 7.0527039 0.4064951457
## [212,] 7.0636549 0.4064951457
## [213,] 7.0718684 0.4228369901
## [214,] 7.1129365 0.4392342339
## [215,] 7.1430526 0.4392342339
## [216,] 7.1594796 0.4392342339
## [217,] 7.1841207 0.4392342339
## [218,] 7.2388773 0.4392342339
## [219,] 7.2744694 0.4392342339
## [220,] 7.2991104 0.4392342339
## [221,] 7.3620806 0.4564431645
## [222,] 7.3702941 0.4564431645
## [223,] 7.4277892 0.4564431645
## [224,] 7.4496918 0.4564431645
## [225,] 7.5811090 0.4743917211
## [226,] 7.5893226 0.4743917211
## [227,] 7.6550307 0.4925685673
## [228,] 7.6577687 0.4925685673
## [229,] 7.7618070 0.4925685673
## [230,] 7.7946610 0.5112049404
## [231,] 7.8384666 0.5112049404
## [232,] 7.8576317 0.5112049404
## [233,] 7.9151268 0.5112049404
## [234,] 8.0602331 0.5112049404
## [235,] 8.1478443 0.5112049404
## [236,] 8.1861734 0.5112049404
## [237,] 8.1998634 0.5112049404
## [238,] 8.3504448 0.5112049404
## [239,] 8.3750858 0.5112049404
## [240,] 8.4024639 0.5112049404
## [241,] 8.4490080 0.5322467712
## [242,] 8.4599590 0.5534919977
## [243,] 8.4654350 0.5534919977
## [244,] 8.4818621 0.5534919977
## [245,] 8.4845991 0.5534919977
## [246,] 8.6214924 0.5534919977
## [247,] 8.6242304 0.5534919977
## [248,] 8.6789865 0.5779303741
## [249,] 8.8213549 0.6030142389
## [250,] 8.8678989 0.6030142389
## [251,] 8.8815880 0.6291000125
## [252,] 8.9117041 0.6291000125
## [253,] 8.9856262 0.6561359330
## [254,] 9.0266943 0.6561359330
## [255,] 9.1334705 0.6561359330
## [256,] 9.1937027 0.6883523143
## [257,] 9.2758389 0.6883523143
## [258,] 9.2950039 0.7238004666
## [259,] 9.3689251 0.7238004666
## [260,] 9.3853521 0.7615571061
## [261,] 9.4318962 0.8018409351
## [262,] 9.4674883 0.8018409351
## [263,] 9.7850790 0.8459689158
## [264,] 9.7932920 0.8459689158
## [265,] 9.8042440 0.8459689158
## [266,] 9.8124571 0.8976938443
## [267,] 9.8863792 0.8976938443
## [268,] 10.0533876 0.8976938443
## [269,] 10.1492128 0.8976938443
## [270,] 10.2997942 0.9535661255
## [271,] 10.4585896 0.9535661255
## [272,] 10.4668036 0.9535661255
## [273,] 10.5106096 1.0162357007
## [274,] 10.5407257 1.0162357007
## [275,] 10.5489388 1.0883390064
## [276,] 10.7132101 1.0883390064
## [277,] 10.7679672 1.0883390064
## [278,] 10.9295006 1.0883390064
## [279,] 11.0198498 1.0883390064
## [280,] 11.0390148 1.0883390064
## [281,] 11.0581789 1.0883390064
## [282,] 11.0882959 1.0883390064
## [283,] 11.1676931 1.1878621915
## [284,] 11.2991104 1.1878621915
## [285,] 11.4551678 1.1878621915
## [286,] 11.4715948 1.1878621915
## [287,] 11.4743328 1.3003477453
## [288,] 11.4880219 1.3003477453
## [289,] 11.5865841 1.3003477453
## [290,] 11.6522932 1.3003477453
## [291,] 11.9507189 1.3003477453
## [292,] 12.1204653 1.3003477453
## [293,] 12.1916494 1.3003477453
## [294,] 12.2080765 1.3003477453
## [295,] 12.2299795 1.3003477453
## [296,] 12.3203287 1.3003477453
## [297,] 12.3449688 1.3003477453
## [298,] 12.3832989 1.3003477453
## [299,] 12.4736481 1.3003477453