S_exp <- function(alpha, t, delta){
log_L <- -sum(delta * log(alpha)) - (1/alpha) * sum(t)
return(-log_L)
}
ajuste <- optimize(f = S_exp, c(0.01, 100), tempos, censura)
ajuste$minimum
[1] 20.41178
$objective
[1] 68.273895 Funções no optim
5.0.1 Exponencial (Criado)
5.0.2 Exponencial (Survival)
modelo_exp <- survreg(Surv(tempos, censura) ~ 1, dist = "exponential")
summary(modelo_exp)
Call:
survreg(formula = Surv(tempos, censura) ~ 1, dist = "exponential")
Value Std. Error z p
(Intercept) 3.016 0.243 12.4 <2e-16
Scale fixed at 1
Exponential distribution
Loglik(model)= -68.3 Loglik(intercept only)= -68.3
Number of Newton-Raphson Iterations: 4
n= 20 5.1 Weibull (Criado)
S_weib <- function(theta, t, delta){
alpha <- theta[1]
gamma <- theta[2]
r <- sum(delta) #soma de todos as falhas
Log_L <- r*log(gamma)-(r*gamma*log(alpha))+(gamma -1)*sum(delta*log(t))-(alpha^(-gamma))*sum(t^gamma)
if(!is.finite(Log_L)) return(1e10) #garante o tempo positivo apenas
return(-Log_L)
}
#?optim()
resultado <- optim(par = c(1, 1), S_weib, t = tempos, delta = censura, method = "BFGS")
resultado $par
[1] 21.334902 1.543351
$value
[1] 66.13336
$counts
function gradient
44 23
$convergence
[1] 0
$message
NULL5.2 Weibull (Survival)
modelo_weib <- survreg(Surv(tempos, censura) ~ 1, dist = "weibull")
summary(modelo_weib)
Call:
survreg(formula = Surv(tempos, censura) ~ 1, dist = "weibull")
Value Std. Error z p
(Intercept) 3.061 0.160 19.1 <2e-16
Log(scale) -0.434 0.189 -2.3 0.022
Scale= 0.648
Weibull distribution
Loglik(model)= -66.1 Loglik(intercept only)= -66.1
Number of Newton-Raphson Iterations: 6
n= 20 6 Log-Normal (Criado)
S_lognorm <- function(theta, t, delta) {
mu <- theta[1]
sigma <- theta[2]
log_S <- plnorm(t, meanlog = mu, sdlog = sigma, lower.tail = FALSE, log.p = TRUE)
log_f <- -0.5*log(2*pi) - log(t) - log(sigma) - (1/(2*sigma^2)) * (log(t) - mu)^2
L_LOG <- sum(delta * log_f + (1 - delta) * log_S)
return(-L_LOG)
}o_mu <- mean(tempos) #mean(log(tempos))
o_sigma <- sd(tempos) #sd(log(tempos))
fit_ln <- optim(par = c(o_mu, o_sigma), fn = S_lognorm,t = tempos,delta = censura, method = "BFGS"
)Warning in plnorm(t, meanlog = mu, sdlog = sigma, lower.tail = FALSE, log.p =
TRUE): NaNs produzidosWarning in log(sigma): NaNs produzidosfit_ln$par
[1] 2.7171762 0.7648184
$value
[1] 65.7399
$counts
function gradient
37 17
$convergence
[1] 0
$message
NULL6.1 Log-Normal (Survival)
modelo_logN <- survreg(Surv(tempos, censura) ~ 1, dist = "lognormal")
summary(modelo_logN)
Call:
survreg(formula = Surv(tempos, censura) ~ 1, dist = "lognormal")
Value Std. Error z p
(Intercept) 2.717 0.176 15.42 <2e-16
Log(scale) -0.268 0.174 -1.54 0.12
Scale= 0.765
Log Normal distribution
Loglik(model)= -65.7 Loglik(intercept only)= -65.7
Number of Newton-Raphson Iterations: 6
n= 20 6.2 Log-Logística (Criada)
S_loglogistica <- function(theta, t, delta) {
alpha <- exp(theta[1])
gamma <- exp(theta[2])
log_f <- log(gamma) - (gamma * log(alpha)) + (gamma - 1) * log(t) - 2 * log(1 + (t/alpha)^gamma)
log_S <- -log(1 + (t/alpha)^gamma)
LOG_L <- sum(delta * log_f + (1 - delta) * log_S)
return(-LOG_L)
}
# Chute inicial (log da média dos tempos e log de 1)
inicio <- c(log(mean(tempos)), 0)
fit_ll <- optim(par = inicio, fn = S_loglogistica, t = tempos, delta = censura, method = "BFGS")
fit_ll$par
[1] 2.7373688 0.7993457
$value
[1] 66.03053
$counts
function gradient
21 7
$convergence
[1] 0
$message
NULL6.3 Log-Logística (Survival)
modelo_LL <- survreg(Surv(tempos, censura) ~ 1, dist = "loglogistic")
modelo_LLCall:
survreg(formula = Surv(tempos, censura) ~ 1, dist = "loglogistic")
Coefficients:
(Intercept)
2.737415
Scale= 0.4496653
Loglik(model)= -66 Loglik(intercept only)= -66
n= 20 6.4 Gama (Criada)
S_gama <- function(theta, t, delta) {
shape <- exp(theta[1])
scale <- exp(theta[2])
log_f <- dgamma(t, shape = shape, scale = scale, log = TRUE)
log_S <- pgamma(t, shape = shape, scale = scale, lower.tail = FALSE, log.p = TRUE)
logl <- sum(delta * log_f + (1 - delta) * log_S)
return(-logl)
}
p <- log(mean(tempos))
fit_g <- optim(par = c(0, p), fn = S_gama, t = tempos, delta = censura, method = "BFGS")
fit_g$par
[1] 0.7626296 2.1937980
$value
[1] 65.84754
$counts
function gradient
30 9
$convergence
[1] 0
$message
NULL6.5 Gama (Flexsurvreg)
modelo_gama <- flexsurvreg(Surv(tempos, censura) ~ 1, dist = "gamma")
modelo_gamaCall:
flexsurvreg(formula = Surv(tempos, censura) ~ 1, dist = "gamma")
Estimates:
est L95% U95% se
shape 2.1439 1.1580 3.9692 0.6737
rate 0.1115 0.0540 0.2303 0.0413
N = 20, Events: 17, Censored: 3
Total time at risk: 347
Log-likelihood = -65.84754, df = 2
AIC = 135.69516.6 Gama Generalizada (Criada)
S_gamaG <- function(theta, t, delta) {
alpha <- exp(theta[1])
gamma <- exp(theta[2])
s <- exp(theta[3])
log_f <- log(gamma)-lgamma(s)-(gamma*s*log(alpha))+(gamma * s - 1)*log(t)-(t / alpha)^gamma
log_S <- pgamma((t / alpha)^gamma, shape = s, lower.tail = FALSE, log.p = TRUE)
logl <- sum(delta * log_f + (1 - delta) * log_S)
if(!is.finite(logl)) return(1e10)
return(-logl)
}
i <- c(log(mean(tempos)), 0, 0)
resultado_gg <- optim(par = i, fn = S_gamaG, t = tempos, delta = censura, method = "BFGS")
resultado_gg$par
[1] -2.8700735 -0.8786525 2.3651648
$value
[1] 65.6936
$counts
function gradient
202 100
$convergence
[1] 1
$message
NULL6.7 Gama Generalizada (Flexsurvreg)
modelo_GG <- flexsurvreg(Surv(tempos, censura) ~ 1, dist = "gengamma")
modelo_GGCall:
flexsurvreg(formula = Surv(tempos, censura) ~ 1, dist = "gengamma")
Estimates:
est L95% U95% se
mu 2.805 2.168 3.442 0.325
sigma 0.743 0.498 1.110 0.152
Q 0.247 -1.291 1.786 0.785
N = 20, Events: 17, Censored: 3
Total time at risk: 347
Log-likelihood = -65.69074, df = 3
AIC = 137.3815