Proyecto Aplicación

Modelo Beta Prime

Introdución

La distribución beta prime (BP) es una distribución de dos parámetros en la línea real positiva, que puede ser interpretada como la distribución de la razón de posibilidades de una variable distribuida según la distribución beta, es decir, si $X \sim Beta(\alpha, \beta)$, entonces $Y = \displaystyle \frac{X}{1-X} \sim BP(\alpha, \beta)$ con $\alpha > 0$ y $\beta >0$ estamos adoptando la parametrización de la distribución BP en términos de la media y los parámetros de precisión que fue propuesta por Bourguignon et al. (2018).

Una ventaja de usar esta parametrización es que podemos introducir estructuras de regresión para cada media y parámetros de precisión y la interpretación de los coeficientes de regresión es sencilla en términos de ellos como en los modelos lineales generalizados.

Distribución

La variable aleatoria $Y$ se distribuye BP (Bourguignon et al. (2018)) con función de densidad de probabilidad dada por: \[ f(y;\mu , \phi)= \frac{y^{\mu(\phi+1)-1}(1+y)^{[-\mu(\phi+1)+\phi+2]}}{B(\mu(1+\phi),\phi+2)},y>0\]

donde $\mu > 0$ y $\phi > 0$ son parámetros de media y precisión, respectivamente, $B(\mu(1+\phi),\phi+2) = \Gamma (\mu(1+\phi))\Gamma(\phi+2)$ es la función beta y $\gamma(\mu(1+\phi)) = \int_{0}^{\infty} \omega^{\mu(1+\phi)-1}e^{-\omega} \cdot d\omega$ es la función gamma. De ahora en adelante, usamos la $Y \sim BP(\mu,\phi)$ para indicar que $Y$ es un variable aleatoria que sigue una distribución BP. La media y la varianza de $Y$ son: \[ E[Y]=\mu , Var[Y]=\frac{\mu(1+\mu)}{\phi}\]

Modelo de regresión Beta prime

Considere $n$ variables aleatorias independientes $Y_{1}, . . . , Y_n$ donde cada $Y_i$ , $i = 1, . . . , n$ tiene PB distribución con media $\mu_i$ y parámetro de precisión $\phi_i$.\ Bourguignon et al. (2018) propusieron el modelo de regresión BP que se define por dos relaciones funcionales \[g_1(\mu_i)= \eta_{1i}=x_i^T \beta ,g_2(\mu_i)= \eta_{2i}=z_i^T v\]

Donde $g_1:\mathbb{R}\longrightarrow \mathbb{R}^+$ y $g_2:\mathbb{R}\longrightarrow \mathbb{R}^+$ son funciones de enlace estrictamente monótonas, positivas y al menos dos veces diferenciables, $\eta_{1i}$ y $\eta_{2i}$ son los predictores lineales, $\beta=(\beta_1,...\beta_p)^T(\beta \in \mathbb{R}^p)$ y $V = (v_1, ..., v_q)^T (V \in \mathbb{R}^q,q<n-p)$ son vectores de parámetros desconocidos a estimar, y $x_i=(x_{i1}, … ,x_{ip})^T $ y $z_i=(z_{i1}, ... ,z_{iq})^T$ son observaciones sobre p y q regresores conocidos, para $i = 1, . . . , n$. Además, suponemos que las matrices de covariables $X = (x_1 ,...,x_n)^T$ y $z = (z_1 ,...,z_n)^T$ tienen rango p y q, respectivamente.

MLE

En la literatura existen pocos trabajos que traten de la distribución de la PB. McDonald (1987) discutió sus propiedades y obtuvo las estimaciones de máxima verosimilitud (ML) de los parámetros del modelo. Las versiones corregidas de sesgo de los MLE de los parámetros que indexan la distribución de BP fueron obtenido por $Stosíc$ y $Cordeiro (2009)$. Vale la pena mencionar que todos los trabajos relacionados anteriormente tienennconsiderado la parametrización habitual de la distribución de la BP. Teniendo en cuenta la parametrización adoptamos, Bourguignon et al. (2018) utilizaron el método ML para estimar los parámetros que indexan el modelo de regresión de BP. Sin embargo, en muestras de pequeño tamaño, los estimadores ML de estos parámetros (especialmente para estructura de precisión) puede ser extremadamente parcial. Por lo tanto, es importante considerar estimadores alternativos con sesgos más pequeños.

La función logarítmica de verosimilitud para $(\beta, V)$ dados los valores observados $y_1, . . . , y_n$ es \[\ell(\beta,V)=\sum_{i=1}^{n}\ell(\mu_i, \phi_i) \] siendo \[ \ell(\mu_i,\phi_i)= [\mu_i(1+\phi_i)-1]log(y_i)-[\mu_i(1+\phi_i)+\phi_i + 2]log(1+y_i)\] \[log[\Gamma(\mu_i(1+\phi_i))]-log[\Gamma(\phi_i+2)]+log[\Gamma(\mu_i(1+\phi_i)+\phi_i+2)]\] $\mu_i=g_1^{-1}(\eta_{1i})$ y $\mu_i=g_2^{-1}(\eta_{2i})$ son funciones de $\beta$ y $V$, respectivamente como se define en el modelo de regresión. Se describe un método para obtener las estimaciones de los parámetros del modelo BP definido por la distribución y en la regresión del modelo se describe en detalle en Bourguignon et al. (2018)

MLE en R

Funciones Preliminares

####################################################################
#
# Librerias
#
###################################################################@

library(Formula)
library(LaplacesDemon)
library(ggplot2)
library(extraDistr)
library(gamlss)

library(Matrix)
library(nleqslv)

############################################################################
#              Función de densidad de probabilidad
###########################################################################
dBP <- function(x,mu=1,sigma=1,log=FALSE)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivo", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(x <= 0)) 
    stop(paste("x debe ser positivo", "\n", ""))
  
  a <- mu*(1+sigma)
  b <- 2 + sigma 
  
  fy <- dbetapr(x, shape1 = a, shape2 = b, scale = 1, log = log)
  fy
  
}

#############################################################################
#               Función de distribución acumulativa 
############################################################################
pBP <-  function(q,mu=1,sigma=1, lower.tail = TRUE, log.p = FALSE)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivoe", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(q < 0)) 
    stop(paste("q debe ser positivo", "\n", ""))
  
  a <- mu*(1+sigma)
  b <- 2 + sigma
  
  cdf <- pbetapr(q, shape1 = a, shape2 = b, scale=1, lower.tail = lower.tail, 
               log.p = log.p)
  cdf
}

#############################################################################
#                              Números Aleatorios
#############################################################################
rBP <- function(n,mu=1,sigma=1)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivo", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(n <= 0)) 
    stop(paste("n debe ser entero positivo", "\n", ""))
  
  n <- ceiling(n)
  
  a <- mu*(1+sigma)
  b <- 2 + sigma  
  
  r <- rbetapr(n,shape1=a,shape2=b,scale=1)
  
  r
}

#############################################################################
#                           Función Cuantil
############################################################################
qBP <- function(p, mu = 1, sigma = 1, lower.tail = TRUE, log.p = FALSE)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivo", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(n <= 0)) 
    if (any(p <= 0) | any(p >= 1)) 
      stop(paste("p debe estar entre 0 y 1", "\n", ""))  
  
  a <- mu*(1+sigma)
  b <- 2 + sigma 
  
  q <- qbetapr(p, shape1 = a, shape2 = b,scale=1, lower.tail = lower.tail, log.p = log.p)
  
  q
}

#############################################################################
#                           familia
############################################################################
BP <- function (mu.link = "log", sigma.link = "log") 
{
  mstats <- checklink("mu.link", "Beta Prime", substitute(mu.link), c("log", "identity", "sqrt"))
  dstats <- checklink("sigma.link", "Beta Prime",substitute(sigma.link), c("log", "identity", "sqrt"))
  structure(list(family = c("BP", "Beta Prime"), 
                 parameters = list(mu = TRUE,sigma = TRUE), nopar = 2, type = "Continuous", 
                 mu.link = as.character(substitute(mu.link)), 
                 sigma.link = as.character(substitute(sigma.link)), 
                 mu.linkfun = mstats$linkfun, 
                 sigma.linkfun = dstats$linkfun, 
                 mu.linkinv = mstats$linkinv, 
                 sigma.linkinv = dstats$linkinv, 
                 mu.dr = mstats$mu.eta, 
                 sigma.dr = dstats$mu.eta, 
                 
                 
                 dldm = function(y, mu, sigma){
                 a <- mu*(1+sigma)
                 b <- mu*(1+sigma)+sigma+2
                 Phi <-  (1+sigma)  
                 yast <- log(y) - log(1+y) 
                 muast <- digamma(a) - digamma(b) 
                 dldm <- Phi*(yast - muast)
                 
                 dldm
                 }, 
                 d2ldm2 = function(mu, sigma){ 
                   Phi2 <- (1+sigma)^2
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2
                   d2dldm2 <- -Phi2*(trigamma(a) - trigamma(b))
                   
                   d2dldm2
                   },
                 dldd = function(y, mu, sigma){ 
                   Phi <-  (1+sigma)
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2 
                   ystar <- mu*log(y) - (1+mu)*log(1+y) 
                   mustar <- mu*digamma(a) - (1+mu)*digamma(b) + digamma(Phi+1)
                   
                   dldd <- ystar - mustar
                   
                   dldd
                  }, 
                 d2ldd2 = function(mu,sigma){
                   Phi <-  (1+sigma)
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2 
                   
                   d2ldd2 <- -(mu^2)*trigamma(a) + ((1+mu)^2)*trigamma(b) - trigamma(Phi+1)
                   
                   d2ldd2
                   
                }, 
                 d2ldmdd = function(mu,sigma){
                   
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2
                   Phi <-  (1+sigma)  
                   gammaast <- Phi*(trigamma(b) + mu*(trigamma(b)-trigamma(a)))
                   
                   d2ldmdd <- gammaast
                   
                   d2ldmdd
                
                   }, 
                 G.dev.incr = function(y, mu, sigma,...){-2*dBP(y, mu, sigma, log = TRUE)}, 
                 rqres = expression(rqres(pfun = "pBP", type = "Continuous", y = y, mu = mu, sigma = sigma)), 
                 mu.initial = expression({mu <- mean(y)}), 
                 sigma.initial = expression({sigma <-  mean(y)*(1+mean(y))/var(y) }), 
                 mu.valid = function(mu) all(mu > 0), 
                 sigma.valid = function(sigma) all(sigma > 0), 
                 y.valid = function(y) all(y > 0)),
                class = c("gamlss.family","family"))
}

#####################################################################
#          BIAS OF THE BETA AND NU - COX-SNELL
#####################################################################
bias <- function(x){ 
  beta1 = matrix(x[1:p], ncol=1)
  nu1 = matrix(x[-(1:p)], ncol=1)
  eta1.est = X%*%beta1
  eta2.est = Z%*%nu1
      
    mu.est <- as.vector(exp(eta1.est))
    phi.est <- as.vector(exp(eta2.est))

  dg1 <- psigamma(mu.est*(1+phi.est), deriv=1)
    dg11 <- psigamma(mu.est*(1+phi.est) + phi.est + 2, deriv=1)
  dg12 <- psigamma(phi.est + 2, deriv=1)
    dg2 <- psigamma(mu.est*(1+phi.est), deriv=2)
    dg21 <- psigamma(mu.est*(1+phi.est) + phi.est + 2, deriv=2)
  dg22 <- psigamma(phi.est + 2, deriv=2)

    a1 <- dg1 - dg11
  b1 <- mu.est^2*dg1 - (1+mu.est)^2*dg11 + dg12
    c1 <- dg2 - dg21
    d1 <- (1+mu.est)^2*dg21 - mu.est^2*dg2
    e1 <- (1+mu.est)^3*dg21 - mu.est^3*dg2 - dg22

  m1 <- -( (1+phi.est)^3*c1*(mu.est)^3 )/2 - ( (1+phi.est)^2*a1*mu.est*mu.est )/2 
    m2 <- ( (1+phi.est)*(-a1*mu.est + dg11)*mu.est*phi.est )/2 - ( (1+phi.est)^2*(mu.est*c1 - dg21)*mu.est^2*phi.est )/2
    m3 <- - ( (1+phi.est)*( 2*a1 + (1+phi.est)*mu.est*c1 - (1+phi.est)*dg21 )*mu.est^2*phi.est )/2  - ( (1+phi.est)*(dg11 - a1*mu.est)*mu.est*phi.est )/2
    m4 <- ( (1+phi.est)*d1 + 2*dg11 - 2*mu.est*a1 )*mu.est*phi.est^2/2 - ( (1+phi.est)*(dg11 - mu.est*a1)*mu.est*phi.est )/2
    m5 <- ( (1+phi.est)*d1/2 )*phi.est^2*mu.est + ( (1+phi.est)*(dg11 - mu.est*a1)*mu.est*phi.est )/2
    m6 <- (e1*phi.est^3 - b1*phi.est*phi.est)/2

  M1 <- diag(m1)
    M2 <- diag(m2)
    M3 <- diag(m3)
    M4 <- diag(m4)
    M5 <- diag(m5)
    M6 <- diag(m6)
    
  E1 <- diag((1+phi.est)^2*(a1)*mu.est^2)
  E2 <- diag((1+phi.est)*(mu.est*a1 - dg11)*mu.est*phi.est)
  E4 <- diag(b1*phi.est^2)

  K11 <- t(X)%*%E1%*%X
    K12 <- t(X)%*%E2%*%Z
    K21 <- t(K12)
    K22 <- t(Z)%*%E4%*%Z
  k1 <- cbind(K11, K12)
  k2 <- cbind(K21, K22)
  IE <- rbind(k1, k2)
  K <- solve(IE)

    Kbb <- K[1:p,1:p]
  Kbnu <- K[1:p,(p+1):(p+q)]
  Knub <- t(Kbnu)
  Knunu <- K[(p+1):(p+q),(p+1):(p+q)]
    Pbb <- matrix(diag(X%*%Kbb%*%t(X)))
    Pbnu <- matrix(diag(X%*%Kbnu%*%t(Z)))
    #Pnub <- t(Pbnu)
    Pnunu <- matrix(diag(Z%*%Knunu%*%t(Z)))
    
  biasbeta <- Kbb%*%t(X)%*%(M1%*%Pbb +(M2+M3)%*%Pbnu + M5%*%Pnunu) + Kbnu%*%t(Z)%*%(M2%*%Pbb +(M4+M5)%*%Pbnu + M6%*%Pnunu)

    biasnu <- Knub%*%t(X)%*%(M1%*%Pbb +(M2+M3)%*%Pbnu + M5%*%Pnunu) + Knunu%*%t(Z)%*%(M2%*%Pbb +(M4+M5)%*%Pbnu + M6%*%Pnunu)
      
  return(c(biasbeta, biasnu))
}

#####################################################################
#          VECTOR SCORE - DAVID FIRTH
#####################################################################
score_DF <- function(x, v){
  U           = numeric(p+q)
  beta.df     = matrix(x[1:p], ncol=1)
  nu.df       = matrix(x[-(1:p)], ncol=1)
  eta1.est.df = X%*%beta.df
  eta2.est.df = Z%*%nu.df
      
    mu.est.df   = as.vector(exp(eta1.est.df))
    phi.est.df  = as.vector(exp(eta2.est.df))
   
  dg1         = psigamma(mu.est.df*(1+phi.est.df), deriv=1)
    dg11        = psigamma(mu.est.df*(1+phi.est.df) + phi.est.df + 2, deriv=1)
  dg12        = psigamma(phi.est.df + 2, deriv=1)

  a1 = dg1 - dg11
  b1 = mu.est.df^2*dg1 - (1+mu.est.df)^2*dg11 + dg12

    E1 = diag((1+phi.est.df)^2*(a1)*mu.est.df^2)
  E2 = diag((1+phi.est.df)*(mu.est.df*a1 - dg11)*mu.est.df*phi.est.df)
  E4 = diag(b1*phi.est.df^2)

  K11 = t(X)%*%E1%*%X
    K12 = t(X)%*%E2%*%Z
    K21 = t(K12)
    K22 = t(Z)%*%E4%*%Z
  k1  = cbind(K11, K12)
  k2  = cbind(K21, K22)
  IE  = rbind(k1, k2)
   
  A        = matrix(x, ncol=p+q)
  res.bias = bias(A)
  vies     = matrix(c(res.bias), ncol=1) 
    
  M       = IE%*%vies
  PHI     = diag(1+phi.est.df)
  D1      = diag(mu.est.df)
  y.star  = matrix(log(v/(1+v)), ncol=1)
  mu.star = matrix(psigamma(mu.est.df*(1+phi.est.df), deriv=0) - psigamma(mu.est.df*(1+phi.est.df) + phi.est.df + 2, deriv=0), ncol=1)
  D2      = diag(phi.est.df)
  v.nu    = matrix(mu.est.df*(y.star-mu.star) + psigamma(mu.est.df*(1+phi.est.df) + phi.est.df + 2, deriv=0) - log(1+v) - psigamma(2+phi.est.df, deriv=0), ncol=1 )
 
  U.beta  = t(X)%*%PHI%*%D1%*%(y.star-mu.star)
  U.nu    = t(Z)%*%D2%*%v.nu
  U       = c(as.vector(U.beta) - M[1:p] , as.vector(U.nu) - M[(p+1):(p+q)] )
    # print(U)
  U
    }

library(COVID19)
library(xts)
library(stats)

###############################################################################
#             Datos 
###############################################################################

data_h <- covid19("Honduras", verbose = FALSE)  # download the data

## Registered S3 method overwritten by 'bit':
##   method   from  
##   print.ri gamlss

ts <- xts(data_h[,c("confirmed", "deaths", "recovered")], order.by = data_h$date)  # convertir a un objeto xts

ts_2 <- na.omit(ts)
ts_r <- c(ts_2$confirmed)
head(ts_r)

##            confirmed
## 2020-03-28        95
## 2020-04-05       268
## 2020-04-10       382
## 2020-04-15       419
## 2020-04-17       442
## 2020-04-19       472

rb = 5000              # Repeticiones para Bootstrap
p = 3                  # dimensiones de beta
q = 2                  # dimensiones de gamma
n  = length(ts_2$recovered)
X  =  matrix(c(rep(1,n), ts_2$confirmed, ts_2$deaths), ncol = p)
Z  = matrix(c(rep(1,n), ts_2$confirmed), ncol = q)
Dobs = ts_2$recovered
n

## [1] 416

#####################################################################
#          INFORMACIÓN DE FISHER
#####################################################################
InfFisher <- function(x){
  beta1 = matrix(x[1:p], ncol=1)
  nu1 = matrix(x[-(1:p)], ncol=1)
  eta1.est = X%*%beta1
  eta2.est = Z%*%nu1
  
  mu.est <- as.vector(exp(eta1.est))
  phi.est <- as.vector(exp(eta2.est))
  
  dg1 <- psigamma(mu.est*(1+phi.est), deriv=1)
  dg11 <- psigamma(mu.est*(1+phi.est) + phi.est + 2, deriv=1)
  dg12 <- psigamma(phi.est + 2, deriv=1)
  dg2 <- psigamma(mu.est*(1+phi.est), deriv=2)
  dg21 <- psigamma(mu.est*(1+phi.est) + phi.est + 2, deriv=2)
  dg22 <- psigamma(phi.est + 2, deriv=2)
  
  a1 <- dg1 - dg11
  b1 <- mu.est^2*dg1 - (1+mu.est)^2*dg11 + dg12
  c1 <- dg2 - dg21
  d1 <- (1+mu.est)^2*dg21 - mu.est^2*dg2
  e1 <- (1+mu.est)^3*dg21 - mu.est^3*dg2 - dg22
  
  E1 <- diag((1+phi.est)^2*(a1)*mu.est^2)
  E2 <- diag((1+phi.est)*(mu.est*a1 - dg11)*mu.est*phi.est)
  E4 <- diag(b1*phi.est^2)
  
  K11 <- t(X)%*%E1%*%X
  K12 <- t(X)%*%E2%*%Z
  K21 <- t(K12)
  K22 <- t(Z)%*%E4%*%Z
  k1 <- cbind(K11, K12)
  k2 <- cbind(K21, K22)
  IE <- rbind(k1, k2)
  K <- solve(IE)
  return(diag(K))
  
}

########################################################################################
#                         M.A.I.N   F.U.N.C.T.I.O.N
########################################################################################
conh = gamlss.control(trace = FALSE, autostep = FALSE, save = TRUE)
fit  = gamlss(ts_2$recovered ~ ts_2$confirmed + ts_2$deaths, ts_2$recovered ~ ts_2$confirmed, family = BP(mu.link="log", sigma.link = "log"), control=conh, method=RS())

beta.est.mv = fit$mu.coefficients
nu.est.mv   = fit$sigma.coefficients

A = c(beta.est.mv, nu.est.mv)

#####################################################################
#     DAVID FIRTH
#####################################################################
est.mv = c(beta.est.mv,nu.est.mv)
est_df = nleqslv(est.mv, score_DF, method="Broyden",global="hook",jacobian=TRUE, control=list(btol=1e-7), v=Dobs) #control=list(trace=0,maxit=500,allowSingular=TRUE, cndtol=1e-5)
beta.est.df = est_df$x[1:p]  
nu.est.df = est_df$x[-(1:p)]

#####################################################################
#                          RESULTS 
#####################################################################
result_est = c(beta.est.mv, nu.est.mv, beta.est.df, nu.est.df)

result_est

#####################################################################
#                       PARAMETROS ESTIMADOS
#####################################################################
nEst = matrix(round(c(result_est),2), ncol=2)
colnames(nEst) = c("MLE","Firth")
rownames(nEst) = c("beta0","beta1","beta2","nu0","nu1")

print.est <- function(par){
  cat("\n Beta's estimates: \n")
  print(par[1:p,])
  cat("\n Nu's estimates:\n") 
  print(par[(p+1):(p+q),])
}
print.est(nEst)

####################################################################
#                      ERRORES STANDARD
####################################################################
variance = data.frame(apply(nEst, 2, InfFisher))
erro.pd  = round(sqrt(variance), 4)
rownames(erro.pd) = c("b0","b1","b2","nu0", "nu1")
print.se <- function(par){
  cat("\n standard error - beta's: \n")
  print(par[1:p,])
  cat("\n standard error - nu's:\n") 
  print(par[(p+1):(p+q),])
}
print.se(erro.pd)

Aplicaciones

$\textbf{Análisis del riesgo-rendimiento con el modelo beta en las Cooperativas de Ahorro y Crédito (COAC's) (Andrés Felipe Sandoval Vera 2018)}$

La implementación del Código Orgánico y Financiera a partir del año 2015, generó diferentes tipos de reacciones en sector de la economía ecuatorial sobre todo en el sector de los Seguros y de las Cooperativas de ahorro y crédito, así que proponen analizar este comportamiento en el sector cooperativo como un indicador de Beta Comparable ayuda a visualizar el panorama futuro del riesgo-rendimiento de las cooperativas analizadas frente a su mercado donde se desarrollan, mediante recopilación de datos entregados por la SEPS.

$\textbf{La tasa social de descuento en la evaluación de proyectos de inversión: Aplicación para el Ecuador}$ $\textbf{( José Gabriel Castillo y Donald Zhangallimbay 2021)}$

En este trabajo presentan una alternativa de determinación de la tasa social de descuento basada en el modelo de estimación gamma (Weitzman, 2001), en el contexto de países en desarrollo.

Modelo Gamma

Introdución

La distribución gamma se puede utilizar con gran flexibilidad en el análisis de variables aleatorias positivas. Por lo tanto, los modelos de regresión gamma se aplican en una amplia gama de aplicaciones empíricas, como en el proceso de fijación de tasas en el marco de trabajo de carteras de seguros heterogéneas, que es la función más importante de las aseguradoras (Krishnamoorthy, 2006), y en los ingresos hospitalarios por enfermedades raras donde las series temporales son muy escasas debido a la infrecuencia de los eventos (Winklemann, 2008).

Distribución gamma

La distribución Gamma tiene la función de densidad de probabilidad:

Si $\displaystyle{X \sim \Gamma (\alpha ,\lambda )}$ entonces su función de densidad es

\[{\displaystyle f(x)={\displaystyle\frac {\lambda (\lambda x)^{\alpha -1}e^{-\lambda x}}{\Gamma (\alpha )}}}\] para $\displaystyle x >0$ $\beta > 0$ es el parámetro de escala y $\alpha > 0$ es el parámetro de forma donde: $\Gamma (\alpha ) = \displaystyle\int_{0}^{\infty }t^{\alpha -1}e^{-t}dt$

Modelo de regresión gamma

Sea $Yi \sim G(\mu_i, \alpha)$, $i = 1, . . . , n$, sea una muestra aleatoria de tamaño $n$. En los modelos de regresión gamma con parámetro de forma constante, la estructura de regresión media se define por: \[g(\mu_i) = \eta_i = x_i'\beta\] donde $g$ es la función de enlace, $\beta = (\beta_0, . . ., \beta_p)'$ es el vector de media de los parámetros de regresión, $x_i$ es el i-ésimo valor del vector de las variables explicativas y $\eta_i$ es un predictor lineal. Algunas funciones de vínculo usuales en la regresión gamma son:

la función logarítmica, $g(\mu)$ = $log(\mu)$;
la función identidad, $g(\mu) = \mu$,
la función inversa $g(\mu) = 1/\mu$.

En los modelos lineales generalizados, el vínculo canónico de la media es la función inversa.

MLE

Introducción y Función Log-Verosimilitud

\[L(\alpha, \lambda) = \displaystyle\prod_{i=0}^n\left[ \displaystyle\frac{\lambda^{\alpha}e^{-\alpha x_i}(x_i)^{\alpha-1}}{\Gamma(\alpha)}\right] \] \[log(L(\alpha, \lambda)) = \log \left( \displaystyle \displaystyle\frac{\lambda^{n\alpha}e^{-\alpha \sum x_i}( \displaystyle\prod_{i=1}^n x_i)^{\alpha-1}}{[\Gamma(\alpha)]^n} \right)\] \[n\cdot\alpha log(\lambda) - \lambda\sum x_i + (\alpha - 1)log(\displaystyle\prod_{i=1}^n x_i) - log([\Gamma(\alpha)]^n) \] \[n\cdot\alpha log(\lambda) - \lambda\sum x_i + (\alpha - 1)\displaystyle\sum_{i=1}^n log(x_i) - n\cdot log(\Gamma(\alpha)) \]

MLE de $\lambda$

\[\displaystyle\frac{\partial l(\alpha, \lambda)}{\partial \lambda} = \displaystyle\frac{n\cdot\alpha}{\lambda} - \lambda\sum x_i = 0\] \[\displaystyle\frac{n\cdot\alpha}{\lambda} = \lambda\sum x_i \Rightarrow \lambda = \displaystyle\frac{n\cdot\alpha}{\lambda\sum x_i} \]

\[\hat{\lambda} =\displaystyle\frac{n\cdot\hat{\alpha}}{\sum x_i} = \displaystyle\frac{n\cdot\hat{\alpha}}{n \cdot \bar{x}} = \displaystyle\frac{\hat{\alpha}}{\bar{x}}\]

MLE de $\alpha$

\[\displaystyle\frac{\partial l(\alpha, \lambda)}{\partial \alpha} = n \cdot log(\lambda) + \displaystyle\sum_{i=1}^nlog(x_i) - \displaystyle\frac{n \cdot\Gamma'(\alpha)}{\Gamma(\alpha)} \] \[n \cdot log(\lambda) + \displaystyle\sum_{i=1}^nlog(x_i) - n \cdot DigammaFuntion(\alpha)\] \[n \cdot log(\displaystyle\frac{\hat{\alpha}}{\bar{x}}) + \displaystyle\sum_{i=1}^nlog(x_i) - n \cdot DigammaFuntion(\alpha)\] \[\hat{\alpha} = n \cdot log(\hat{\alpha}) - n \cdot log(\bar{x}) + \displaystyle\sum_{i=1}^nlog(x_i) - n \cdot DigammaFuntion(\hat{\alpha})\] \[\hat{\alpha} = n [ log(\hat{\alpha}) - log(\bar{x}) + \displaystyle\frac{1}{n}\sum_{i=1}^nlog(x_i) - DigammaFuntion(\hat{\alpha}) ]\] \[log(\hat{\alpha}) - log(\bar{x}) + \displaystyle\frac{1}{n}\sum_{i=1}^nlog(x_i) - DigammaFuntion(\hat{\alpha}) = 0\]

MLE en R

library(COVID19)
library(xts)
library(stats)

data_h <- covid19("Honduras", verbose = FALSE)  # download the data

ts <- xts(data_h[,c("confirmed", "deaths", "recovered")], order.by = data_h$date)  # convertir a un objeto xts

Omitiremos los NA

ts_2 <- na.omit(ts)
head(ts_2)

##            confirmed deaths recovered
## 2020-03-28        95      1         3
## 2020-04-05       268     22         6
## 2020-04-10       382     23         7
## 2020-04-15       419     31         9
## 2020-04-17       442     41        10
## 2020-04-19       472     46        15

Estimaremos los $\alpha$ y $\lambda$ para los datosde personas recuperadas.

# Sea x_bar la media muestral de los datos.
x_bar  <- mean(ts_2$recovered) # alfa/lambda es la media


# Tome el logaritmo natural de cada uno de los valores
lnSample <- log(ts_2$recovered)


# Obtener la media de los datos de muestra de registro natural
meanlnSample <- mean(lnSample)


# Escriba una función que corresponda a la derivada parcial que derivamos de la función de verosimilitud logarítmica
# con respecto a alfa.

f <- function(x){
  # En este caso de funciones, tomamos el valor x, que será el valor de alfa que estamos resolviendo.
  log(x) - digamma(x) - log(x_bar) + meanlnSample
}

# Use la función uniroot para encontrar dónde la derivada es igual a cero para que podamos encontrar
# la raíz (valor x o sombrero alfa) en la que ocurre.
a <- uniroot(f,lower = 0.5, upper = 2.5)

# Establecer alpha_hat igual a la raíz de la función f en el intervalo (2.0,2.2)
alpha_hat <- a$root
alpha_hat

## [1] 0.6132038

# Recuerde que se encontró que la estimación de probabilidad máxima de Lambda era alpha_hat sobre x_bar.
lambda_hat <- alpha_hat/x_bar

lambda_hat

## [1] 1.430313e-05

GLM

summary(glm(ts_2$recovered ~ ts_2$confirmed + ts_2$deaths, family=Gamma))

## 
## Call:
## glm(formula = ts_2$recovered ~ ts_2$confirmed + ts_2$deaths, 
##     family = Gamma)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8931  -0.9076  -0.0058   0.4683   0.8597  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     6.282e-05  2.504e-06   25.09   <2e-16 ***
## ts_2$confirmed -1.098e-09  9.919e-11  -11.07   <2e-16 ***
## ts_2$deaths     3.407e-08  3.594e-09    9.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.4828182)
## 
##     Null deviance: 835.68  on 415  degrees of freedom
## Residual deviance: 577.22  on 413  degrees of freedom
## AIC: 9455
## 
## Number of Fisher Scoring iterations: 6

Aplicaciones

$\textbf{Modelación de la distribución gamma en matlab para aplicaciones de radar ( José Raul Machado-Fernandez 2016)}$

El clutter es una señal de interferencia que aparece en los sistemas de radar. El autor modeló la distribución Gamma, íntimamente relacionado con clutter, en MATLAB. La implementación permite un acceso fácil a la manipulación de las funciones de densidad de probabilidad y de distribución, la generación de muestras, el cálculo de los momentos, la aplicación de algoritmos de bondad de ajuste y la estimación de parámetros a partir de muestras con distribución gamma.

$\textbf{La tasa social de descuento en la evaluación de proyectos de inversión: Aplicación para el Ecuador}$ $\textbf{( José Gabriel Castillo y Donald Zhangallimbay 2021)}$

En este trabajo presentan una alternativa de determinación de la tasa social de descuento basada en el modelo de estimación gamma (Weitzman, 2001), en el contexto de países en desarrollo.

$\textbf{Reclamaciones de seguros}$

La distribución gamma se ha utilizado para modelar el tamaño de las reclamaciones de seguros$^{[1]}$ y las precipitaciones$^{[2]}$. Esto significa que las reclamaciones de seguros agregadas y la cantidad de lluvia acumulada en un embalse se modelan mediante un proceso gamma , al igual que la distribución exponencial genera un proceso de Poisson.

$\textbf{Comunicación Inalámbrica}$

La distribución gamma se utiliza para modelar el desvanecimiento de la potencia de la señal por múltiples rutas.

$\textbf{Oncología}$

La distribución por edades de la incidencia de cáncer a menudo sigue la distribución gamma, mientras que los parámetros de forma y escala predicen, respectivamente, el número de eventos impulsores y el intervalo de tiempo entre ellos$^{[3]}$.

$\textbf{Neurociencia}$

la distribución gamma se usa a menudo para describir la distribución de intervalos entre picos$^{[4][5]}$.

$\textbf{Expresión génica bacteriana}$

El número de copias de una proteína expresada constitutivamente a menudo sigue la distribución gamma, donde el parámetro de escala y forma son, respectivamente, el número medio de explosiones por ciclo celular y el número medio de moléculas de proteína producidas por un solo ARNm durante su vida$^{[6]}$.

$\textbf{En genómica}$

La distribución gamma se aplicó en el paso de llamada de pico (es decir, en el reconocimiento de la señal) en el análisis de datos de ChIP-chip$^{[7]}$ y ChIP-seq$^{[8]}$.

$\textbf{Estadística bayesiana}$

La distribución gamma se usa ampliamente como un previo conjugado en la estadística bayesiana. Es el conjugado previo para la precisión (es decir, inversa de la varianza) de una distribución normal. También es el conjugado previo de la distribución exponencial.

Notas

p. 43, Philip J. Boland, Métodos estadísticos y probabilísticos en ciencia actuarial, Chapman & Hall CRC 2007
Aksoy, H. (2000) “Uso de distribución gamma en análisis hidrológico” , Turk J. Engin Environ Sci , 24, 419 - 428.
Belikov, Aleksey V. (22 de septiembre de 2017). “El número de eventos cancerígenos clave se puede predecir a partir de la incidencia del cáncer” . Informes científicos. 7 (1): 12170. doi : 10.1038 / s41598-017-12448-7. PMC 5610194 . PMID 28939880 .
JG Robson y JB Troy, “Naturaleza de la descarga mantenida de las células ganglionares retinianas Q, X e Y del gato”, J. Opt. Soc. Soy. A 4, 2301-2307 (1987)
MCM Wright, IM Winter, JJ Forster, S. Bleeck “La respuesta a las ráfagas de tono de la mejor frecuencia en el núcleo coclear ventral se rige por las estadísticas ordenadas del intervalo entre picos”, Hearing Research 317 (2014)
N. Friedman, L. Cai y XS Xie (2006) “Vinculación de la dinámica estocástica con la distribución de la población: un marco analítico de expresión génica”, Phys. Rev. Lett. 97, 168302.
DJ Reiss, MT Facciotti y NS Baliga (2008) “Deconvolución basada en modelos de la unión del ADN en todo el genoma” , Bioinformática , 24, 396–403.
MA Mendoza-Parra, M Nowicka, W Van Gool, H Gronemeyer (2013) “Caracterización de patrones de unión de ChIP-seq por deconvolución de forma de pico basada en modelo” , BMC Genomics , 14: 834

Dase de Datos

Introducción

Podemos analizar los datos mundiales de casos reportados para la enfermedad por coronavirus (CoViD-19) del repositorio del Centro de Ciencias e Ingeniería de Sistemas de la Universidad Johns Hopkins (JHU CSSE) https://github.com/CSSEGISandData/COVID -19.

Los conjuntos de datos están disponibles en dos modalidades, como secuencias de series de tiempo y según el estado de la persona (confirmado, fallecido o recuperado).

Gracias a este paquete tenemos disponibles varias funciones de análisis, visualización y modelado que nos permitirán calcular y visualizar el número total de casos, el número total de cambios y la tasa de crecimiento a nivel mundial o para una ubicación geográfica específica.

Incluso tenemos para generar el modelo Susceptible-Infected-Recovered (SIR) para la propagación de la enfermedad. (beta)

Instalamos la librería en nuestro entorno de R

install.packages('covid19.analytics')

Cargamos la librería para poder usarla

library(covid19.analytics)

Explicación de la base de datos

La función covid19.data() permite a los usuarios obtener datos en tiempo real sobre los casos notificados por CoViD19 del repositorio CCSE de la JHU, en las siguientes modalidades:

aggregated para el último día, con una gran ‘granularidad’ de regiones geográficas (es decir, ciudades, provincias, estados, países).
time series para regiones geográficas acumuladas más grandes (provincias / países).
deprecated: también incluimos el estilo de datos original en el que se informaron inicialmente estos conjuntos de datos.

Podemos jugar con estas secuencias de tiempo y los diferentes estados.

$\textbf{Time Series data}$

ts-confirmed Datos confirmados

ts-deaths Datos de fallecidos

ts-recovered Datos de recuperados

ts-ALL Datos combinados

Vamos a obtener por secuencia de tiempo todos los datos confirmados por COVID19.

# Cargamos toda la base de datos.
data <- covid19.data()

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
## ================================================================================ 
## ================================================================================ 
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

data_c <- covid19.data(case="ts-confirmed")

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
## --------------------------------------------------------------------------------

No mostraremos la data devido a que es demasiado grande, pero mostraremos los casos confirmados para un país en específico, además nos mostrará un modelo de regresión lineal donde intenta predecir la variable $y$ y utiliza los valores de la variable $x$.

tots.per.location(data_c, geo.loc = "Honduras")

## HONDURAS  --  472250 
## ===============================   running models...=============================== 
##   Linear Regression (lm): 
## 
## Call:
## lm(formula = y.var ~ x.var)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -83171 -26344  -5143  28531  78550 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -26322.571   2226.421  -11.82   <2e-16 ***
## x.var          508.962      3.372  150.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37610 on 1141 degrees of freedom
## Multiple R-squared:  0.9523, Adjusted R-squared:  0.9523 
## F-statistic: 2.279e+04 on 1 and 1141 DF,  p-value: < 2.2e-16
## 
## -------------------------------------------------------------------------------- 
##   Linear Regression (lm): 
## 
## Call:
## lm(formula = y.var ~ x.var)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.9217 -0.8980  0.5388  1.6095  2.0237 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.5999522  0.1305452   58.22   <2e-16 ***
## x.var       0.0065658  0.0001977   33.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.205 on 1141 degrees of freedom
## Multiple R-squared:  0.4915, Adjusted R-squared:  0.4911 
## F-statistic:  1103 on 1 and 1141 DF,  p-value: < 2.2e-16
## 
## -------------------------------------------------------------------------------- 
##   GLM using Family [1] "poisson" : 
## 
## Call:
## glm(formula = y.var ~ x.var, family = family)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -384.46  -171.54    -4.13   139.70   257.87  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.106e+01  1.657e-04   66720   <2e-16 ***
## x.var       2.099e-03  1.989e-07   10552   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 170543659  on 1142  degrees of freedom
## Residual deviance:  43454804  on 1141  degrees of freedom
## AIC: 43469796
## 
## Number of Fisher Scoring iterations: 5
## 
## --------------------------------------------------------------------------------

Tambien podemos sacar los datos de fallecidos a nivel mundial.

data_d <- covid19.data(case="ts-deaths")

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
## --------------------------------------------------------------------------------

live.map(data_d)

Proyecto Aplicación

Terlyn Ulloa

2022-08-16

Modelo Beta Prime

Introdución

Distribución

Modelo de regresión Beta prime

MLE

MLE en R

Aplicaciones

Modelo Gamma

Introdución

Distribución gamma

Modelo de regresión gamma

MLE

Introducción y Función Log-Verosimilitud

MLE de \(\lambda\)

MLE de \(\alpha\)

MLE en R

GLM

Aplicaciones

Notas

Dase de Datos

Introducción

Instalamos la librería en nuestro entorno de R

Cargamos la librería para poder usarla

Explicación de la base de datos