Distribuição Logística e seus modelos de regressão
aplicações e exemplos
Marcos S. Medeiros
22 de junho de 2022
Introdução
Para variáveis respostas do tipo dicotômicas é comum a codificação 0 (ausência) e 1 (presença). Nesse caso, o usuário define qual das situações corresponde ao “sucesso” ou evento de interesse. Esse conceito nos remete à variável aleatória de Bernoulli, em que uma variável assume apenas os valores 0 ou 1. O exemplo clássico é o jogo da moeda, em que o sucesso corresponde a cara, com probabilidade \(p=50\%\), assumindo que a moeda não seja viciada.
Podemos, a partir daí, modelar genericamente a contagem do número de vezes em que a moeda resultou “cara”. Analogamente, em um grupo de pacientes submetidos a uma determinada intervenção, podemos observar quantos apresentaram determinado desfecho. Esse é o caso especial de uma variável aleatória binomial em que \(N=1\) e tem \(S=\{0,1\}\), mais conhecida como a variável de Bernoulli. O parâmetro \(p\) pode ser interpretado como a probabilidade de um indivíduo no grupo apresentar o evento de interesse, variando entre 0 e 1.
Variável aleatória binomial genérica:
\(Prob[X=x]=f(x)= \binom{N}{x}p^{x}(1-p)^{N-x}\) em que \(x=0,1,2,...,N\) e \(\binom{N}{x}=\frac{N!}{x!(N-x)!}\)
Exemplos
- Um voo numa determinada companhia aérea num Boeing 747 tem 410 assentos. Sabendo que a probabilidade de no-show é de 0.3 para aquele voo, calcule a probabilidade de overbooking nos casos em que a empresa tenha vendido 510 passagens e 560 passagens.
success <- 0:410
plot(success, dbinom(success, size=510, prob=0.7),type='h')
pbinom(410,510,0.7, lower.tail = FALSE)## [1] 3.739798e-08pbinom(410,560,0.7, lower.tail = FALSE)## [1] 0.04283985- Caso a probabilidade de no-show se mantenha constante, simule um experimento de ocupação de 30 voos dessa companhia, com venda de 560 passagens.
set.seed(91271)
voos <- rbinom(30, size=560, prob=0.7)
voos## [1] 374 401 410 403 389 417 390 396 385 409 383 373 378 401 384 388 361 383 393
## [20] 401 404 393 391 399 395 411 396 385 394 408mean(voos)## [1] 393.1667Como aponta Aljarrah (2020), a distribuição logística é a distribuição limitante da média dos maiores e menores valores de uma amostra aleatória de tamanho \(n\) de uma distribuição simétrica do tipo exponencial. Khan (2017) apresenta uma visualização das propriedades das variáveis de forma e escala:
# Khan, R. (2017)
pdf=function(x,mu,s){
k=(x-mu)/s
return(exp(-k)/(s*(1+exp(-k))^2))
}
cdf=function(x,mu,s){
k=(x-mu)/s
return(1/(1+exp(-k)))
}
x=seq(-10,10,0.01)
## PDF
layout(matrix(1:2,nrow=1))
plot(x,pdf(x,-1,1),type="l",lty=3,lwd=2,ylab="PDF")
lines(x,pdf(x,0,1),type="l",lty=1,lwd=2,col="gray50")
lines(x,pdf(x,1,1),type="l",lty=5,lwd=2,col="gray70")
abline(v=c(0,1,-1),col=2, lty=3)
legend("topright",c("\u03BC = -1, \u03C3 = 1","\u03BC = 0, \u03C3 = 1","\u03BC = 1, \u03C3 = 1"),
bty="n",lty=c(3,1,5),col=c("black","gray50","gray70"),lwd=2, cex=0.7)
## CDF
plot(x,cdf(x,-1,1),type="l",lty=3,lwd=2,ylab="CDF")
lines(x,cdf(x,0,1),type="l",lty=1,lwd=2,col="gray40")
lines(x,cdf(x,1,1),type="l",lty=5,lwd=2,col="gray70")
abline(h=0.5, v=c(0,1,-1),col=2, lty=3)
legend("topleft",c("\u03BC = -1, \u03C3 = 1","\u03BC = 0, \u03C3 = 1","\u03BC = 1, \u03C3 = 1"),
bty="n",lty=c(3,1,5),col=c("black","gray40","gray70"),lwd=2, cex=0.7)
Modelo de regressão logística
Consideremos o modelo de regressão linear tradicional:
\[y_{i}=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+... +\beta_{k}x_{ki}+\varepsilon_{i}\]
Nesse modelo, a variável dependente (VD) \(y_{i}\) é, via de regra, uma variável contínua com distribuição normal. Contudo, esses pressupostos não necessariamente precisam ser verdadeiros para que possamos utilizar um modelo de regressão linear.
Existem diversas situações em que é interessante o uso de um modelo de regressão adaptado, de acordo com a variável resposta. Uma dessas situações corresponde a uma variável binária, por exemplo, se um leito de hospital está ocupado ou livre, se um paciente submetido a uma intervenção apresentou ou não um determinado desfecho.
É possível adaptar o conceito de Bernoulli ao modelo de regressão. Consideremos uma base de dados de unidades observacionais em que cada unidade assume uma variável resposta de valor 0 ou 1 e um conjunto de variáveis potencialmente explicativas.
A ideia básica consiste em testar a probabilidade \(Prob[Y_{i}=1]=p_{i}=\beta_{0}+\beta_{1}x_{1i}+\beta_{2}x_{2i}+... +\beta_{k}x_{ki}\) para o intervalo \([0,1]\), o que implica: \[Prob[Y_{i}=1]=p_{i}=(\frac{\beta_{0}+\beta_{1}x_{1i}+\beta_{2}x_{2i}+... +\beta_{k}x_{ki}}{1+\beta_{0}+\beta_{1}x_{1i}+\beta_{2}x_{2i}+... +\beta_{k}x_{ki}})\]
Se \(\beta_{1}\) é positivo, quando \(x_{1}\) aumenta, a probabilidade de sucesso também aumenta.
O desfecho ou variável dependente (VD) ou ainda, variável resposta \(y_{i}\) da regressão logística pode ser binária ou categórica (ordenada ou não), porém a variável explicativa ou preditora \(x_{i}\) pode ser binária, categórica ou contínua.
Exemplos
Apresentou febre nos últimos 7 dias? (0 - não; 1 - sim)
Concorda com a afirmação seguinte? (1-não concordo; 2-neutro; 3-concordo)
Dentre os dois modelos seguintes, qual você prefere? (1-A; 2-B, 3-qualquer um)
Probability Density Function (PDF) e Cumulative Density Function (CDF).
A PDF de uma distribuição logística é definida por:
\[f_{X}(x)= \frac{e^{-\frac{x-\mu }{\sigma}}}{s(1+e^{-\frac{x-\mu}{\sigma}})^2}\] em que \(-\infty<x<\infty\).
A CDF de um distribuição logística é definida por:
\[f_{X}(x)= \int_{-\infty}^{x}xf_{X}(x)dx= \frac{1}{1+e^{-\frac{x-\mu}{\sigma}}}\]
ou \(f_{X}(x)= (1+exp(-\frac{x-\mu}{\sigma}))^{-1}\), em que \(-\infty<x<\infty\).
Exemplo
Dados originários de Hosmer e Lemeschow (2013), por meio do pacote “aplore3”:
Amostra de 100 pessoas, em que a variável dependente é a ocorrência ou não de doença coronária cardíaca (CHD), associando-se à idade (AGE).
library(aplore3)
library(ggplot2)
head(chdage)## id age agegrp chd
## 1 1 20 20-39 No
## 2 2 23 20-39 No
## 3 3 24 20-39 No
## 4 4 25 20-39 No
## 5 5 25 20-39 Yes
## 6 6 26 20-39 Nosummary(chdage)## id age agegrp chd
## Min. : 1.00 Min. :20.00 55-59 :17 No :57
## 1st Qu.: 25.75 1st Qu.:34.75 30-34 :15 Yes:43
## Median : 50.50 Median :44.00 40-44 :15
## Mean : 50.50 Mean :44.38 45-49 :13
## 3rd Qu.: 75.25 3rd Qu.:55.00 35-39 :12
## Max. :100.00 Max. :69.00 20-39 :10
## (Other):18ggplot(chdage, aes(x=age, y=chd)) + theme_bw() +
geom_point() +
stat_smooth(method="glm", method.args=list(family="binomial"), se=FALSE)## `geom_smooth()` using formula 'y ~ x'## Warning: Computation failed in `stat_smooth()`:
## y values must be 0 <= y <= 1
A partir desse conjunto de observações, pode-se modelar a regressão logística:
m1=glm(chd~age, family = binomial(link="logit"), data = chdage)
summary(m1)##
## Call:
## glm(formula = chd ~ age, family = binomial(link = "logit"), data = chdage)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9718 -0.8456 -0.4576 0.8253 2.2859
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.30945 1.13365 -4.683 2.82e-06 ***
## age 0.11092 0.02406 4.610 4.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 136.66 on 99 degrees of freedom
## Residual deviance: 107.35 on 98 degrees of freedom
## AIC: 111.35
##
## Number of Fisher Scoring iterations: 4Observa-se que o valor do intercepto é -5,309 para idade igual a zero. Embora não seja um resultado prático, o importante é a relação com o coeficiente idade, de 0,1109. Pelo fato de ser positiva, significa que há um incremento de chance de ocorrência de chd conforme a idade aumenta, sendo p<0.001, o que é estatisticamente significante. Assim, teremos:
\[ln(\frac{p_{chd}}{1-p_{chd}})=-5.309+0.1109age\]
Podemos construir a curva CDF:
IDADE <- chdage[,1]
# Criando a varável preditora
chdage$PRED=predict(m1, data=IDADE, type="response")
# Plotando o gráfico CDF
ggplot(chdage, aes(x=age, y=PRED)) + theme_bw() +
geom_point()
Propriedades e pesquisas recentes
A CDF de uma distribuição logística é \(F(y)=(1+e^{-y})^-1, -\infty<y<\infty\). A função de densidade com curtose 4.2 é simétrica em relação a zero e possuem caudas mais longas e pesadas que as da função de densidade normal, o que torna a distribuição logística uma escolha popular para dados não normais de ajuste aproximadamente simétrico.
O primeiro tipo de distribuição para valores extremos é conhecido como Distribuição de Gumbel devido ao trabalho desse autor (1958), que fez diversas contribuições para a análise de valores extremos, com aplicações na linha de tempo de vida, emissões radioativas e análise de enchentes.
As CDFs de valores máximos e mínimos de Gumbel são definidas, respectivamente, como sendo: \[f_{Gu-max}(x;\mu,\sigma)= exp\left \{ -exp\left ( -\frac{x-\mu}{\sigma} \right ) \right \}, -\infty<x<\infty, -\infty<\mu<\infty, \sigma>0\]
\[f_{Gu-min}(x;\mu,\sigma)= 1 - exp\left \{ -exp\left ( -\frac{x-\mu}{\sigma} \right ) \right \}, -\infty<x<\infty, -\infty<\mu<\infty, \sigma>0\]
É interessante notar que a distribuição de Gumbel é uma boa escolha para ajustar dados assimétricos enquanto que a distribuição logística se amolda a dados simétricos. Uma propriedade notável é a relação entre ambas as distribuições: Se \(X\sim Gumbel(\mu_{X},\sigma)\) e \(Y \sim Gumbel(\mu_{Y}, \sigma)\) , então \((X-Y) \sim logistic(\mu_{x}-\mu_{Y},\sigma)\).
Quando buscamos estabelecer uma relação de dependência ou intercorrelação entre variáveis aleatórias no intervalo [0,1], usamos o modelo de cópula, descrito por Abe Sklar (1959), cujo teorema afirma que “qualquer distribuição conjunta multivariada pode ser escrita em termos de funções de distribuição marginal univariadas e uma cópula, que descreve a estrutura de dependência entre as variáveis”.
Exemplos:
library(gumbel)
#independência
rand <- t(rgumbel(200,1))
plot(rand[1,], rand[2,], col="green", main="Gumbel copula")
#dependência positiva
rand <- t(rgumbel(200,2))
plot(rand[1,], rand[2,], col="red", main="Gumbel copula")
#Densidade de distribuição 3D
x <- seq(.05, .95, length = 30)
y <- x
z <- outer(x, y, dgumbel, alpha=2)
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 100, shade = 0.25, ticktype = "detailed",
xlab = "x", ylab = "y", zlab = "density")
#CDF 3D
z <- outer(x, y, pgumbel, alpha=2)
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 100, shade = 0.25, ticktype = "detailed",
xlab = "u", ylab = "v", zlab = "CDF")
O Modelo de E-L{GW}, descrito por Aljarrah (2020)
Aljarrah (2020) propõe um modelo da família das distribuições logísticas generalizadas, a distribuição Exponencial-Logística Generalizada de Weibull (E-L {GW}).
Seja \(R\) uma variável aleatória em uma distribição logística padrão. Usando um parâmetro de forma \(\xi>0\), parâmetro de localização \(-\infty<\mu<\infty\), parâmetro de reflexão de escala \(\sigma\neq0\), usando a técnica descrita pelo mesmo autor (2019) para combinar a função exponencial-normal generalizada de Weibull, a distribuição E-L{GW} é definida por:
\[F_{X}(x)=0.5+sng(\sigma)(0.5-exp\left \{ -((\bar{F}_{R}(\frac{x-\mu}{\sigma}))^{-\xi}-1) /\xi\right \})\] em que \(sng(\sigma)\) é o sinal (positivo ou negativo) do parâmetro \(\sigma\). Observe que a CDF reduz-se ao termo \(\bar{F}_{R}(\frac{x-\mu}{\left | \sigma\right |})\) quando \(\xi\rightarrow0\)
O modelo possui três parâmetros, sendo \(\xi>0\) o parâmetro de forma, \(-\infty<\mu<\infty\) o parâmetro de localização e \(\sigma\neq0\) o parâmetro de reflexão de escala.
A PDF correspondente é dada por:
\[f_{X}(x)= \frac{1}{\left | \sigma \right |}exp(\frac{x-\mu}{\sigma})(1+(\frac{x-\mu}{\sigma}))^{\xi-1}exp\left \{ (1-(1+ exp{(\frac{x-\mu}{\sigma})} )^{\xi} /\xi\right \}\]
em que \(-\infty<x,\mu<\infty\) , \(\sigma\neq0\) e \(\xi>0\).
Propriedades
Quando \(\xi\rightarrow0\) a distribuição logística é \(\bar{F}_{R}(\frac{x-\mu}{\left | \sigma\right |})\)
Quando \(\xi=1\) e \(\sigma<0\), a PDF assume a distribuição máxima de Gumbel.
Quando \(\xi=1\) e \(\sigma>0\), a PDF assume a distribuição mínima de Gumbel.
Comportamento
##### PDF
f_PDF<-function(par,x)
{
mu=par[1]
sg=par[2]
ep=par[3]
z=sign(sg)*(x-mu)/abs(sg)
n1=dlogis(z, location = 0, scale = 1)
n2=plogis(z, location = 0, scale = 1)
c1=((1-n2)^-ep-1)/ep
c2=1/(abs(sg)*(1-n2)^(ep+1))
f=n1*c2*exp(-c1)
return(f)
}
pelgw <- f_PDF
curve(f_PDF(c(0, 1, 0.2), x), col = "blue", lwd = 2, xlim = c(-5, 5), ylim = c(0, 1), xlab = "x",
ylab = "fx(x)", main = "Comportamento da Curva PDF sgn(\u03C3) positivo")
curve(f_PDF(c(0, 1, 1), x), col = "red", lwd = 2, add = TRUE)
curve(f_PDF(c(0, 1, 2), x), col = "green", lwd = 2, add = TRUE)
curve(f_PDF(c(0, 1,3), x), col = "orange", lwd = 2, add = TRUE)
curve(f_PDF(c(0, 1, 10), x), col = "violet", lwd = 2, add = TRUE)
legend("topright",
legend = c("\u03BE = 0.2 \u03C3 = 1",
"\u03BE = 1, \u03C3 = 1",
"\u03BE = 2, \u03C3 = 1",
"\u03BE = 3, \u03C3 = 1",
"\u03BE = 10, \u03C3 = 1"),
col = c("blue", "red", "green", "orange", "violet"), lty = 1, lwd = 2)
curve(f_PDF(c(0, -1, 0.2), x), col = "blue", lwd = 2, xlim = c(-5, 5), ylim = c(0, 1), xlab = "x",
ylab = "fx(x)", main = "Comportamento da Curva PDF sgn(\u03C3) negativo")
curve(f_PDF(c(0, -1, 1), x), col = "red", lwd = 2, add = TRUE)
curve(f_PDF(c(0, -1, 2), x), col = "green", lwd = 2, add = TRUE)
curve(f_PDF(c(0, -1,3), x), col = "orange", lwd = 2, add = TRUE)
curve(f_PDF(c(0, -1, 10), x), col = "violet", lwd = 2, add = TRUE)
legend("topleft",
legend = c("\u03BE = 0.2 \u03C3 = -1",
"\u03BE = 1, \u03C3 = -1",
"\u03BE = 2, \u03C3 = -1",
"\u03BE = 3, \u03C3 = -1",
"\u03BE = 10, \u03C3 = -1"),
col = c("blue", "red", "green", "orange", "violet"), lty = 1, lwd = 2)
##### CDF
F_CDF<-function(par,x)
{
mu=par[1]
sg=par[2]
ep=par[3]
z=(x-mu)/abs(sg)
n1=dlogis(z, location = 0, scale = 1)
n2=plogis(z, location = 0, scale = 1)
c1=((1-n2)^-ep-1)/ep
c2=((n2)^-ep-1)/ep
F=(sg <0)*exp(-c2)+(sg>0)*(1-exp(-c1))
return(F)
}
curve(F_CDF(c(0, 1, 0.2), x), col = "blue", lwd = 2, xlim = c(-5, 5), ylim = c(0, 1), xlab = "y",
ylab = "f(x|a, b)", main = "Comportamento da Curva CDF sgn(\u03C3) positivo")
curve(F_CDF(c(0, 1, 1), x), col = "red", lwd = 2, add = TRUE)
curve(F_CDF(c(0, 1, 2), x), col = "green", lwd = 2, add = TRUE)
curve(F_CDF(c(0, 1, 3), x), col = "orange", lwd = 2, add = TRUE)
curve(F_CDF(c(0, 1, 10), x), col = "violet", lwd = 2, add = TRUE)
legend("topleft",
legend = c("\u03BE = 0.2 \u03C3 = 1",
"\u03BE = 1, \u03C3 = 1",
"\u03BE = 2, \u03C3 = 1",
"\u03BE = 3, \u03C3 = 1",
"\u03BE = 10, \u03C3 = 1"),
col = c("blue", "red", "green", "orange", "violet"), lty = 1, lwd = 2)
curve(F_CDF(c(0, -1, 0.2), x), col = "blue", lwd = 2, xlim = c(-5, 5), ylim = c(0, 1), xlab = "y",
ylab = "f(x|a, b)", main = "Comportamento da Curva CDF sgn(\u03C3) negativo")
curve(F_CDF(c(0, -1, 1), x), col = "red", lwd = 2, add = TRUE)
curve(F_CDF(c(0, -1, 2), x), col = "green", lwd = 2, add = TRUE)
curve(F_CDF(c(0, -1, 3), x), col = "orange", lwd = 2, add = TRUE)
curve(F_CDF(c(0, -1, 10), x), col = "violet", lwd = 2, add = TRUE)
legend("topleft",
legend = c("\u03BE = 0.2 \u03C3 = -1",
"\u03BE = 1, \u03C3 = -1",
"\u03BE = 2, \u03C3 = -1",
"\u03BE = 3, \u03C3 = -1",
"\u03BE = 10, \u03C3 = -1"),
col = c("blue", "red", "green", "orange", "violet"), lty = 1, lwd = 2)
#### QELGW
qelgw <- function(par,u)
{
mu=par[1]
sg=par[2]
ep=par[3]
z=sign(sg)*(u-0.5)
c1=log(0.5-z)
c2=log(((1-ep*c1)**(1/ep))-1)
x <- mu + sg*c2
return(x)
}
qelgw(c(-3.36,-8.46,8.69), 1) ## [1] Infcurve(qelgw(c(-3,-12, 8), x), col = "blue", lwd = 2, xlim = c(0, 1), ylim = c(0, 40), xlab = "u",
ylab = "Qx(u|\u03BC, \u03C3, \u03BE)", main = "Comportamento da função quantílica (\u03BC = - 3, \u03BE = 8)")
curve(qelgw(c(-3,-8,8), x), col = "red", lwd = 2, add = TRUE)
curve(qelgw(c(-3,-6,8), x), col = "green", lwd = 2, add = TRUE)
curve(qelgw(c(-3,-4,8), x), col = "orange", lwd = 2, add = TRUE)
curve(qelgw(c(-3,-2,8), x), col = "violet", lwd = 2, add = TRUE)
legend("topleft",
legend = c("\u03C3 = -12",
"\u03C3 = -8",
"\u03C3 = -6",
"\u03C3 = -4",
"\u03C3 = -2"),
col = c("blue", "red", "green", "orange", "violet"), lty = 1, lwd = 2)
curve(qelgw(c(-3,-8, 10), x), col = "blue", lwd = 2, xlim = c(0, 1), ylim = c(0, 40), xlab = "u",
ylab = "Qx(u|\u03BC, \u03C3, \u03BE)", main = "Comportamento da função quantílica (\u03BC = -3, \u03C3 = -8)")
curve(qelgw(c(-3,-8,15), x), col = "red", lwd = 2, add = TRUE)
curve(qelgw(c(-3,-8,20), x), col = "green", lwd = 2, add = TRUE)
curve(qelgw(c(-3,-8,30), x), col = "orange", lwd = 2, add = TRUE)
curve(qelgw(c(-3,-8,50), x), col = "violet", lwd = 2, add = TRUE)
legend("topleft",
legend = c("\u03BE = 10",
"\u03BE = 15",
"\u03BE = 20",
"\u03BE = 30",
"\u03BE = 50"),
col = c("blue", "red", "green", "orange", "violet"), lty = 1, lwd = 2)
Estudo de Simulação
##### pelgw
pelgw<-function(q, mu=1, sg=2, ep=3)
{
z=sign(sg)*(q-mu)/abs(sg)
n1=dlogis(z, location = 0, scale = 1)
n2=plogis(z, location = 0, scale = 1)
c1=((1-n2)^-ep-1)/ep
c2=1/(abs(sg)*(1-n2)^(ep+1))
f=n1*c2*exp(-c1)
return(f)
}
##### delgw
delgw<-function(x, mu=1, sg=2, ep=3)
{
z=(x-mu)/abs(sg)
n1=dlogis(z, location = 0, scale = 1)
n2=plogis(z, location = 0, scale = 1)
c1=((1-n2)^-ep-1)/ep
c2=((n2)^-ep-1)/ep
F=(sg <0)*exp(-c2)+(sg>0)*(1-exp(-c1))
return(F)
}
#### qelgw
qelgw <- function(q, mu=1, sg=2, ep=3)
{
z=sign(sg)*(q-0.5)
c1=log(0.5-z)
c2=log(((1-ep*c1)**(1/ep))-1)
x <- mu + sg*c2
return(x)
}
### relgw
relgw <- function(n, mu=1, sg=2, ep=3)
{
u <- runif(n)
z=sign(sg)*(u-0.5)
c1=log(0.5-z)
c2=log(((1-ep*c1)**(1/ep))-1)
x <- mu + sg*c2
return(x)
}
library(fitdistrplus)
B <- 2000
n <- seq(20, 200, 30)
L <- length(n)
par <- c(0, -1, 3)
mle <- list()
mle <- matrix(ncol = 3, nrow = B)
vies <- rmse <- matrix(ncol = 3, nrow = L)
for(j in 1:L)
{
set.seed(912)
for(i in 1:B)
{
dados3 <- relgw(n[j], par[1], par[2], par[3])
mle[i,] <- fitdist(dados3, distr = 'elgw', method = "mle",
start = list(mu=par[1], sg=par[2], ep=par[3]))$estimate
}
aux <- mle - matrix(par, ncol = 3, nrow = B, byrow = T)
vies[j,] <- apply(aux, 2, mean, na.rm = T)
rmse[j,] <- apply(aux ^ 2, 2, mean, na.rm = T)
cat(i, n[j], "\n")
}## 2000 20
## 2000 50
## 2000 80
## 2000 110
## 2000 140
## 2000 170
## 2000 200par(mfrow = c(1, 2))
plot(n, vies[,1], type = "b", pch = 17, main = "Viés de \u03BC")
plot(n, rmse[,1], type = "b", pch = 17, main = "RMSE de \u03BC")
plot(n, vies[,2], type = "b", pch = 17, main = "Viés de \u03C3")
plot(n, rmse[,2], type = "b", pch = 17, main= "RMSE de \u03C3")
plot(n, vies[,3], type = "b", pch = 17, main = "Viés de \u03BE")
plot(n, rmse[,3], type = "b", pch = 17, main = "RMSE de \u03BE")
Podemos observar que, no estudo de simulação, o viés e erro médio quadrático é considerável para \(n<100\), o que nos impõe muita cautela quando pretendemos afirmar, em dados reais, que uma distribuição segue o modelo E-L{GW}. Por outro lado, quando \(n>100\), a qualidade do ajuste se aproxima consideravelmente ao modelo teórico, como é o caso da aplicação que demonstraremos a seguir.
Aplicações
O modelo E-L{GW} foi demonstrado originalmente em 4 aplicações diferentes, sendo as duas primeiras em R e as duas seguintes em SAS.
Reproduzimos a primeira aplicação, no estudo de Patrício, M., Pereira, J., Crisóstomo, J. et al. Using Resistin, glucose, age and BMI to predict the presence of breast cancer. BMC Cancer 18, 29 (2018).
Trata-se de um conjunto de dados consistente em 116 medições de Adiponectina (\(\mu g/mL\))
library(AdequacyModel)
library(AdequacyModel)
library(rmutil)
library(cubature)
library(QRM)
library(gdata)
library(moments)
library(Bolstad2)
# Conjunto de dados
x <- c(9.7024,5.429285,22.43204,7.16956,4.81924,13.67975,5.589865,13.25132,10.358725,11.57899,13.11,26.72,23.67,36.06,17.95,
20.32,38.04,7.780255,5.46262,5.1,8.6,7.64276,6.209635,9.048185,9.73138,4.617125,4.667645,10.35526,6.966895,5.065915,
7.53955,20.03,4.7942,8.300955,3.19209,2.19428,4.267105,6.796985,2.19428,8.12555,7.36996,9.34663,9.000805,11.78796,
8.13,7.652055,3.74122,10.567295,8.2863,7.9317,4.935635,9.75326,21.056625,8.237405,8.412175,21.823745,8.462915,8.574655,
5.80762,11.236235,4.77192,4.138025,3.886145,9.349775,4.230105,4.294705,3.70523,4.783985,18.55,20.37,16.67,16.1,12.76,
17.86,9.76,6.420295,11.018455,12.71896,5.357135,13.494865,16.44048,10.636525,21.57,5.4861,6.695585,7.28287,10.22231,
7.901685,4.104105,10.26266,6.160995,9.92365,6.26854,10.79394,8.40443,10.73174,5.47817,6.78387,1.65602,5.288025,7.65222,
8.42996,3.71009,2.78491,2.36495,3.335665,6.644245,9.16,11.9,10.06,8.01,12.1,21.42,22.54,33.75,14.11)
dados<-x
# Medida notáveis
In_mean1=mean(dados)
In_mean1 # Média## [1] 10.18087In_sd1=sd(dados)
In_sd1 # Desvio-Padrão## [1] 6.843341sign_sn=-skewness(dados)/abs(skewness(dados))
sign_sn # Sinal## [1] -1In_sd2=sign_sn*sd(dados)
In_sd2 # Sinal*Desvio-Padrão## [1] -6.843341kskewd=skewness(dados)
kskewd # Skewness## [1] 1.794174kurtd=kurtosis(dados)
kurtd # Kurtosis## [1] 6.710771Podemos comparar o modelo de Regressão Logística Padrão sobre o histograma da distribuição dos dados, a fim de verificar sua (in)adequação (“goodness”).
# Distribuição Logística
# PDF
L_PDF<-function(par,x)
{
mu=par[1]
sg=par[2]
z=(x-mu)/(sg)
f=(1/sg)*dlogis(z, location = 0, scale = 1)
return(f)
}
##### CDF
L_CDF<-function(par,x)
{
mu=par[1]
sg=par[2]
z=(x-mu)/(sg)
f=plogis(z, location = 0, scale = 1)
return(f)
}
x = seq(0,ceiling(max(dados)), length.out = 500)
frame()
L_PDF(c(In_mean1,In_sd1),dados)## [1] 0.036487251 0.032459873 0.017912690 0.034818976 0.031453032 0.034244712
## [7] 0.032711155 0.034753272 0.036525693 0.036153290 0.034908443 0.010987553
## [13] 0.015682117 0.003182827 0.026894371 0.022049408 0.002410087 0.035430631
## [19] 0.032512534 0.031926248 0.036048781 0.035303800 0.033621010 0.036282794
## [25] 0.036492487 0.031102744 0.031191025 0.036525931 0.034589195 0.031869664
## [31] 0.035204380 0.022639455 0.031410059 0.035851225 0.028447088 0.026454074
## [37] 0.030478617 0.034386712 0.026454074 0.035720263 0.035033291 0.036396471
## [43] 0.036261626 0.036032776 0.035723722 0.035312577 0.029503994 0.036502756
## [49] 0.035840706 0.035562799 0.031651159 0.036496225 0.020569362 0.035805053
## [55] 0.035928515 0.019067375 0.035962280 0.036033310 0.033041991 0.036315510
## [61] 0.031371720 0.030243260 0.029776617 0.036397488 0.030411427 0.030528592
## [67] 0.029435851 0.031392493 0.025674306 0.021948030 0.029410281 0.030468903
## [73] 0.035264724 0.027075791 0.036497338 0.033907014 0.036395389 0.035303827
## [79] 0.032345019 0.034471080 0.029842247 0.036491402 0.019558991 0.032549471
## [85] 0.034261697 0.034941752 0.036531526 0.035537238 0.030180975 0.036530557
## [91] 0.033553234 0.036518961 0.033702224 0.036458662 0.035923278 0.036472746
## [97] 0.032537011 0.034370716 0.025355463 0.032233898 0.035312733 0.035940457
## [103] 0.029445063 0.027642991 0.026799842 0.028726776 0.034197231 0.036329368
## [109] 0.035961513 0.036529012 0.035627997 0.035822911 0.019852122 0.017711540
## [115] 0.004382040 0.033679170hist(dados, breaks = seq(0, 40, 3),prob=TRUE,xlim=c(-3,42), ylim=c(0,0.12))
curve(L_PDF(c(In_mean1,In_sd1),x), add=T,lty=1,col=1,lwd=3)
box()
Como podemos observar, a curva de distribuição logística não se amolda ao conjunto de dados.
Vamos, agora, modelar o conjunto de dados usando a nova distribuição proposta pelo estudo mencionado:
# Modelo E-L{GW}
# PDF
f_PDF<-function(par,x)
{
mu=par[1]
sg=par[2]
ep=par[3]
z=sign(sg)*(x-mu)/abs(sg)
n1=dlogis(z, location = 0, scale = 1)
n2=plogis(z, location = 0, scale = 1)
c1=((1-n2)^-ep-1)/ep
c2=1/(abs(sg)*(1-n2)^(ep+1))
f=n1*c2*exp(-c1)
return(f)
}
# CDF
F_CDF<-function(par,x)
{
mu=par[1]
sg=par[2]
ep=par[3]
z=(x-mu)/abs(sg)
n1=dlogis(z, location = 0, scale = 1)
n2=plogis(z, location = 0, scale = 1)
c1=((1-n2)^-ep-1)/ep
c2=((n2)^-ep-1)/ep
F=(sg <0)*exp(-c2)+(sg>0)*(1-exp(-c1))
return(F)
}Vamos testar a qualidade do ajuste:
res_L=goodness.fit(pdf=L_PDF,cdf=L_CDF,
starts=c(8.46,8.69),
method="B", data=dados, domain=c(-Inf,Inf),mle=NULL)
res_E_L_GW=goodness.fit(pdf=f_PDF,cdf=F_CDF,
#starts=c(-3000,-1000,10),
#starts=c(-10,-5.803,15), lbest ks values
starts=c(-3.36,-8.46,8.69),
method="B", data=dados, domain=c(-Inf,Inf),mle=NULL)
#Convergence 0 indicates successful completion and 1 indicates that
#the iteration limit maxit had been reached.
print(c(res_L, res_E_L_GW))## $W
## [1] 0.6975544
##
## $A
## [1] 4.083836
##
## $KS
##
## One-sample Kolmogorov-Smirnov test
##
## data: data
## D = 0.10846, p-value = 0.1305
## alternative hypothesis: two-sided
##
##
## $mle
## [1] 9.108284 3.389387
##
## $AIC
## [1] 760.0408
##
## $`CAIC `
## [1] 760.147
##
## $BIC
## [1] 765.5479
##
## $HQIC
## [1] 762.2764
##
## $Erro
## [1] 0.5345197 0.2716438
##
## $Value
## [1] 378.0204
##
## $Convergence
## [1] 0
##
## $W
## [1] 0.04600238
##
## $A
## [1] 0.2907324
##
## $KS
##
## One-sample Kolmogorov-Smirnov test
##
## data: data
## D = 0.044703, p-value = 0.9746
## alternative hypothesis: two-sided
##
##
## $mle
## [1] -3.361406 -8.460868 8.688507
##
## $AIC
## [1] 714.5514
##
## $`CAIC `
## [1] 714.7657
##
## $BIC
## [1] 722.8122
##
## $HQIC
## [1] 717.9048
##
## $Erro
## [1] 9.667282 3.049097 8.918288
##
## $Value
## [1] 354.2757
##
## $Convergence
## [1] 0Pelo Teste de qualidade do ajuste de Cramer-von Mises (W), a qualidade de ajuste é muito superior na distribuição E-L{GW}. Graficamente, podemos visualizar novamente o histograma, agora com a curva de Distribuição E-L{GW}:
x = seq(0,ceiling(max(dados)), length.out = 500)
frame()
#hist(dados,breaks=seq(0,100,l=21), prob=T, main="",xlim=c(0,100), ylim=c(0,0.04))
hist(dados, breaks = seq(0, 40, 3),prob=TRUE,xlim=c(-3,42), ylim=c(0,0.12))
curve(f_PDF(res_E_L_GW$mle,x), add=T,lty=1,col=1,lwd=3)
box()
plot(ecdf(dados), lwd = 1)
curve(F_CDF(c(-3.361406, -8.460868, 8.688507), x),
add=T, lty=1,col="blue",lwd=2)
Conclusão
As novas pesquisas em modelos de distribuição estão longe de exaurimento. A cada dia, novos conjuntos de dados com variáveis resposta com forte grau de interação ainda carecem de modelagem para prever seu comportamento. Em Bioestatística, os desafios são ainda maiores tendo em vista o imenso campo de aplicações. O novo modelo proposto por Aljarrah(2020) demonstrou-se confiável, reproduzível e verificável na vida real, a partir de observações e testes estatísticos.
Referências:
Aljarrah, M.A., Famoye, F., Lee, C.: A new generalized normal distribution: properties and applications. Commun. Stat. Theory Methods. 48(18), 4474–4491 (2019)
Aljarrah, M.A., Famoye, F., Lee, C.: Generalized logistic distribution and its regression model. Journal of Statistical Distributions and Applications (2020) 7:7
Khan, R. (2017) Logistic Distribution Basics
Hosmer Jr, DW., Lemeshow, S., Sturdivant R. Applied Logistic Regression 3rd ed. Hoboken: Wiley, 2013.
Patrício, M., Pereira, J., Crisóstomo, J., Matafome, P., Gomes, M., Seiça, R., Caramelo, F. Using Resistin, glucose, age and BMI to predict the presence of breast cancer. BMC Cancer. 18, 29 (2018).
Smolski, F.M. da S. Regressão Logística.