####################################################################################################################################
####################################################################################################################################
######################################################## tarea-examen 1 ########################################################
####################################################################################################################################
####################################################################################################################################
#### ejercicio 4
#library(OIsurv)
library(fitdistrplus)
## Warning: package 'fitdistrplus' was built under R version 3.5.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 3.5.3
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.5.3
## Loading required package: npsurv
## Loading required package: lsei
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 3.5.3
#library(xm12)
library(rvest)
## Warning: package 'rvest' was built under R version 3.5.3
## Loading required package: xml2
library(stats)
library(survminer)
## Warning: package 'survminer' was built under R version 3.5.3
## Loading required package: ggpubr
## Loading required package: magrittr
library(survival)
library(flexsurv)
## Warning: package 'flexsurv' was built under R version 3.5.3
library(rriskDistributions)
## Warning: package 'rriskDistributions' was built under R version 3.5.3
library(MASS)
#Supongamos que tenemos 21 pacientes con leucemia aguda con los siguientes tiempos de de mejora sustancial de su enfermedad
#La variable es medida en meses
PACIENTES=c(1,1,2,2,3,4,4,5,5,6,8,8,9,10,10,12,14,16,20,24,34)
P=data.frame(PACIENTES)
barplot(PACIENTES) # parece que si ajusta una exponencial
plotdist(PACIENTES,histo=TRUE, demp=TRUE)
descdist(PACIENTES)
## summary statistics
## ------
## min: 1 max: 34
## median: 8
## mean: 9.428571
## estimated sd: 8.387916
## estimated skewness: 1.562348
## estimated kurtosis: 5.562953
ajust = fitdist(PACIENTES, "exp") # Una manera de estimar los parametros de la distribucion ya con informacion adicional
summary(ajust)
## Fitting of the distribution ' exp ' by maximum likelihood
## Parameters :
## estimate Std. Error
## rate 0.1060606 0.02314226
## Loglikelihood: -68.11864 AIC: 138.2373 BIC: 139.2818
plot(ajust)
# si deseamos verlas por separado
denscomp(ajust)
cdfcomp(ajust)
qqcomp(ajust)
ppcomp(ajust)
#a) ?C?al es la tasa de recaida?
#procederemos a encontrar la tasa de recaida, por lo que buscaremos las h(t) pero suponemos una dstribucion exponencial por lo que esta coincidira con la lambda
#sabemos que la h(t) se define como lambda para este caso, entonces busdacaremos el ajuste de un modelo exponencial
#escojemos la exponencial con otra manera de estimar el parametro
Ajuste = fit.cont(data2fit=(PACIENTES))
##
## Begin fitting distributions ---------------------------------------
## * fitting normal distribution ... OK
## * fitting Cauchy distribution ... OK
## * fitting logistic distribution ... OK
## * fitting beta distribution ... failed
## * fitting exponential distribution ... OK
## * fitting chi-square distribution ... OK
## * fitting uniform distribution ... OK
## * fitting gamma distribution ... OK
## * fitting lognormal distribution ... OK
## * fitting Weibull distribution ... OK
## * fitting F-distribution ... OK
## * fitting Student's t-distribution ... OK
## * fitting Gompertz distribution ... OK
## * fitting triangular distribution ... failed
## End fitting distributions -----------------------------------------
## logL AIC BIC Chisq(value) Chisq(p) AD(value)
## Normal -73.95 151.9 153.99 3.84 0.15 0.98
## Cauchy -73.98 151.97 154.06 2.12 0.35 0.93
## Logistic -72.98 149.97 152.06 2.24 0.33 0.67
## Exponential -68.12 138.24 139.28 0.30 0.96 0.29
## Chi-square -78.01 158.02 159.07 25.40 0.00 3.38
## Uniform NULL NULL NULL 5.42 0.07 Inf
## Gamma -67.43 138.86 140.95 1.09 0.58 0.14
## Lognormal -67.59 139.18 141.27 1.75 0.42 0.22
## Weibull -67.57 139.14 141.23 0.82 0.66 0.15
## F -80.07 164.15 166.24 15.64 0.00 5.94
## Student -97.44 196.87 197.92 58.31 0.00 17.25
## Gompertz -67.88 139.77 141.85 0.41 0.81 0.20
## H(AD) KS(value) H(KS)
## Normal rejected 0.19 NULL
## Cauchy not rejected 0.20 NULL
## Logistic rejected 0.16 NULL
## Exponential not rejected 0.11 NULL
## Chi-square NULL 0.18 NULL
## Uniform NULL 0.20 NULL
## Gamma not rejected 0.08 NULL
## Lognormal not rejected 0.12 NULL
## Weibull not rejected 0.08 NULL
## F NULL 0.41 NULL
## Student NULL 0.69 NULL
## Gompertz NULL 0.09 NULL
##
## Chosen continuous distribution is: Normal (norm)
## Fitted parameters are:
## mean sd
## 9.428571 8.185768
# nos arroja que 1/lambda = 0.1060606
# con lo que obtenemos un lambda de 10 lo cual se traduce en una tasa de recaida de 10 meses
# b)
# ?Cu?l es el tiempo medio de remisi?n ?
# lo obtendremos con el comando summary
summary(PACIENTES)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 4.000 8.000 9.429 12.000 34.000
# Con lo que obtenemos un tiempo medio de 9.429 ~ 9.4
#El cu?l nos dice que el tiempo medio de remision es 9.429
# c)
# ?Se puede decir que el tiempo medio de remision es al menos de 5 meses con una conianza del 95% ?
# Si
#d)
# ? Cu?l es la probabilidad de permanecer en remision por lo menos 20 meses ?
#Aqu? utilizaremos la funcion de distribucion para calcular la probabilidad dado que ya conocemos el parametro
proba=pexp(20,10,lower.tail = TRUE)
proba
## [1] 1
# por lo que la probabilidad es de uno.
####################################################################################################################################
####################################################################################################################################
####### Ejercicio 7
# CREAREMOS NUESTRO VECTOR DE TIEMPO DE FALLA ASI COMO EL VECTOR DE CENSURAS
# Cabe mencinar que en este ejercicio en particular consideramos 0:tiempo de falla y 1:cesura
tiempos=c(2,19,19,25,30,35,40,45,45,48,60,62,69,89,90,110,145,160,9,10,20,40,50,110,130)
censura=c(rep(0,18),rep(1,7))
length(tiempos)==length(censura)
## [1] TRUE
censura
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
dat=data.frame(tiempos,censura)
dat
## tiempos censura
## 1 2 0
## 2 19 0
## 3 19 0
## 4 25 0
## 5 30 0
## 6 35 0
## 7 40 0
## 8 45 0
## 9 45 0
## 10 48 0
## 11 60 0
## 12 62 0
## 13 69 0
## 14 89 0
## 15 90 0
## 16 110 0
## 17 145 0
## 18 160 0
## 19 9 1
## 20 10 1
## 21 20 1
## 22 40 1
## 23 50 1
## 24 110 1
## 25 130 1
# nuestros tiempos
dat$tiempos
## [1] 2 19 19 25 30 35 40 45 45 48 60 62 69 89 90 110 145
## [18] 160 9 10 20 40 50 110 130
dat$censura
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
# a)
# Realice las gr?ficas p-p, q-q y de linearizaci?n de la funci?n de sueprvivencia para probar graficamente si los datos siguen una distribucion exponencial con lambda = 0.01
attach(dat)
## The following objects are masked _by_ .GlobalEnv:
##
## censura, tiempos
D1<-Surv(tiempos,censura)
D1
## [1] 2+ 19+ 19+ 25+ 30+ 35+ 40+ 45+ 45+ 48+ 60+ 62+ 69+ 89+
## [15] 90+ 110+ 145+ 160+ 9 10 20 40 50 110 130
f1<-survfit(D1~1)
f1
## Call: survfit(formula = D1 ~ 1)
##
## n events median 0.95LCL 0.95UCL
## 25 7 130 110 NA
summary(f1)
## Call: survfit(formula = D1 ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 9 24 1 0.958 0.0408 0.882 1.000
## 10 23 1 0.917 0.0564 0.813 1.000
## 20 20 1 0.871 0.0698 0.744 1.000
## 40 16 1 0.816 0.0840 0.667 0.999
## 50 11 1 0.742 0.1041 0.564 0.977
## 110 5 1 0.594 0.1567 0.354 0.996
## 130 3 1 0.396 0.1924 0.153 1.000
f1<-survfit(D1~1,conf.int=.95)
summary(f1)
## Call: survfit(formula = D1 ~ 1, conf.int = 0.95)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 9 24 1 0.958 0.0408 0.882 1.000
## 10 23 1 0.917 0.0564 0.813 1.000
## 20 20 1 0.871 0.0698 0.744 1.000
## 40 16 1 0.816 0.0840 0.667 0.999
## 50 11 1 0.742 0.1041 0.564 0.977
## 110 5 1 0.594 0.1567 0.354 0.996
## 130 3 1 0.396 0.1924 0.153 1.000
print(f1,print.rmean=FALSE)
## Call: survfit(formula = D1 ~ 1, conf.int = 0.95)
##
## n events median 0.95LCL 0.95UCL
## 25 7 130 110 NA
#detalles del objeto survfit
str(f1)
## List of 13
## $ n : int 25
## $ time : num [1:21] 2 9 10 19 20 25 30 35 40 45 ...
## $ n.risk : num [1:21] 25 24 23 22 20 19 18 17 16 14 ...
## $ n.event : num [1:21] 0 1 1 0 1 0 0 0 1 0 ...
## $ n.censor : num [1:21] 1 0 0 2 0 1 1 1 1 2 ...
## $ surv : num [1:21] 1 0.958 0.917 0.917 0.871 ...
## $ type : chr "right"
## $ std.err : num [1:21] 0 0.0426 0.0615 0.0615 0.0801 ...
## $ lower : num [1:21] 1 0.882 0.813 0.813 0.744 ...
## $ upper : num [1:21] 1 1 1 1 1 ...
## $ conf.type: chr "log"
## $ conf.int : num 0.95
## $ call : language survfit(formula = D1 ~ 1, conf.int = 0.95)
## - attr(*, "class")= chr "survfit"
#grafico
plot(f1,main="Estimador Kaplan Meier con intervalo del 95% de confianza",
xlab = "tiempo", ylab="Probabilidad de supervivencia")
####################################
#Funcion de riesgo acumulada
D1<-Surv(tiempos)
f1<-summary(survfit(D1~1))
#usando la estimacion de KM de la funcion de supervivencia
H1<- -log(f1$surv)
#usando Nelson-Aalen
h<-f1$n.event/f1$n.risk
H2<-cumsum(h)
#grafico
H1<- c(H1,tail(H1,1))
H2<-c(H2,tail(H2,1))
plot(c(f1$time,250),H1,xlab="tiempo",ylab="riesgo acumulado",type="s")
points(c(f1$time,250),H2,lty=2,type="s")
############
#BANDAS DE CONFIANZA
plot(survfit(D1 ~ 1), xlim=c(0, 400), xlab="tiempo",
ylab="Estimacion Funcion de Supervivencia",
main="Bandas de Confianza")
#######
#Prueba de modelos parametricos
d1<-survfit(D1 ~ 1)
#Exponencial
plot(d1$time,log(d1$surv))
#Weibull
plot(log(d1$time),log(-log(d1$surv)))
library(flexsurv)
fs1 <- flexsurvreg(D1~ 1, data = D1,dist = "weibull")
plot(fs1)
fs1 <- flexsurvreg(D1~ 1, data = D1,dist = "lognormal")
plot(fs1)
fs1 <- flexsurvreg(D1~ 1, data = D1,dist = "exp")
plot(fs1)
# a)
#Grafica P-P
Sw<- exp(-1/as.numeric(fs1$coefficients[2])*d1$time^abs((as.numeric(fs1$coefficients[1]))))
Sw<- exp(-1/100*d1$time)
plot(Sw,d1$surv)
abline(a=0,b=1)
#Grafica Q-Q
p = c(.01,.1,.2,.3,.4,.5,.6,.7,.8,.9)
Qw = qexp(p,.01)
#Qw = qweibull(p, shape=0.18, scale = 1/as.numeric(fs1$coefficients[2]), lower.tail = TRUE, log.p = FALSE)
Qe = quantile(d1$time,probs = p)
plot(Qw,Qe)
abline(a=0,b=1)
# Por lo que podemos concluir que no sigue una distribucion lambda = .01 debido a que en la linealizacion no parte de la identidad y no ajusta bien el modelo.
# b) se sugeririra utilizar una prueba log-rank o en su defecto una ji-cuadrada con una sola poblacion.