Modelo Beta Prime
2/8/2022
Integrantes
Terlyn Ulloa Galeano 20171003625.
Bryan Daniel Colindres 20181000761.
Lizzy Abigail Rodas 20181006261.
Yitxi Sevilla 20182500027.
Nathalye Nicol Deras 20161002923.
Gabriela Eunice Romero Sierra 20201004426.
Base de Datos y Modelo
Ahora presentamos un estudio de simulación de Monte Carlo para investigar y comparar el rendimiento de los MLE junto con sus versiones corregidas propuestas en este proyecto en muestras de tamaño pequeño y moderado.
Usamos modelos de regresión de BP con covariables de dispersión y un enlace logarÃtmico. Consideramos el modelo
\[ log(\mu_i) = \beta_0 + \sum_{j=1}^p\,\beta_j\, x_{ij} \,\,\,y\,\,\, log(\phi_i) = v_0 + \sum_{j=1}^p\,\beta_j\, v_{ij}\] Para \(i=1, 2, .., n\)
Explicacion de los datos
Consideramos el conjunto de datos real utilizado en \(Bonnail\, et\, al.\, (2016)\). El estudio contiene 27 observaciones (muestra de tamaño pequeño), que midieron, entre otras caracterÃsticas
El peso seco del tejido de las almejas (seco, en g).
El peso húmedo del tejido (húmedo, en g).
las concentraciones de cesio (cs ) en el tejido blando.
Dichos minerales se consideraron en 100 microgramos por litro (100 µg/L).
Código
Adoptamos un modelo de BP para ajustar el tejido de peso seco de las almejas.
####################################################################
#
# Librerias
#
###################################################################@
library(Formula)
library(LaplacesDemon)
library(ggplot2)
library(extraDistr)
library(gamlss)############################################################################
# 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"))
}
Funciones de Enlace
MLE
Una implementación de R para la obtención de MLEs junto con sus versiones corregidas propuestas en este
#######################################################################
# Libreria
#######################################################################
library(Matrix)
library(nleqslv)###############################################################################
# Datos
###############################################################################
dry = c(0.2013, 0.1987, 0.2015, 0.2014, 0.2042, 0.2003, 0.2055, 0.2013, 0.2031,
0.2014, 0.1958, 0.2000, 0.2017, 0.1866, 0.1783, 0.2006, 0.2031, 0.2025,
0.2013, 0.2016, 0.2005, 0.2067, 0.1998, 0.2002, 0.1566, 0.1332, 0.2075)
wet = c(1.15, 0.83, 0.96, 1.39, 1.74, 1.41, 1.92, 2.17, 0.98, 1.75, 0.90,
1.22, 1.65, 1.01, 0.90, 1.67, 1.94, 1.79, 1.84, 1.29, 1.73, 1.82,
2.23, 2.05, 0.71, 0.69, 1.81)
cs = c(0.11, 0.09, 0.08, 0.01, 0.0, 0.09, 0.09, 0.10, 0.09, 0.08, 0.09,
0.08, 0.10, 0.11, 0.11, 0.12, 0.10, 0.10, 0.10, 0.11, 0.11, 0.11,
0.10, 0.10, 0.15, 0.18, 0.12)
rb = 1000 # Repetticiones para Bootstrap
p = 3 # dimensiones de beta
q = 2 # dimensiones de gamma
n = length(dry)
X = matrix(c(rep(1,n), wet, cs), ncol = p)
Z = matrix(c(rep(1,n), wet), ncol = q)
Dobs = dry#####################################################################
# 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))
}#####################################################################
# BIAS DE LA BETA Y 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))
}#####################################################################
# BIAS DE LA BETA Y NU - BOOTSTRAP
#####################################################################
bp.fun <- function(data, est.b){
eta1.est.boot <- X%*%est.b[1:p]
eta2.est.boot <- Z%*%est.b[-(1:p)]
mu.est.boot <- as.vector(exp(eta1.est.boot))
nu.est.boot <- as.vector(exp(eta2.est.boot))
est.bootp <- matrix(NA,rb,p+q)
set.seed(2020)
Dobs_b <- matrix(mapply(rBP, rb, mu.est.boot, nu.est.boot), nr = rb)
calc_boot <- function(dadosb){
loglik_boot <- function(vp){
beta_boot <- vp[1:p]
nu_boot <- vp[-(1:p)]
eta1_boot <- as.vector(X%*%beta_boot)
eta2_boot <- as.vector(Z%*%nu_boot)
mu_boot <- exp(eta1_boot)
phi_boot <- exp(eta2_boot)
a_boot <- mu_boot*(1+phi_boot)
b_boot <- 2 + phi_boot
fy_boot <- (a_boot-1)*log(dadosb) - (a_boot+b_boot)*log(1+dadosb) - lbeta(a_boot,b_boot)
return(sum(fy_boot))
}
v0 <- est.b
est_bootp <- try(optim(v0, loglik_boot, method = "BFGS", control=list(fnscale=-1)),T)
if(class(est_bootp)=="try-error"){
estb = rep(NA, p+q)
}else{
estb = est_bootp$par
}
return(estb)
}
calc_bootb = match.fun(calc_boot)
resultb = t(apply(Dobs_b, 1, calc_boot))
return(resultb)
}########################################################################################
# M.A.I.N F.U.N.C.T.I.O.N
########################################################################################
conh = gamlss.control(trace = FALSE, autostep = FALSE, save = TRUE)
fit = gamlss(dry~wet+cs,dry~wet,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)
res.bias = bias(A)#####################################################################
# COX-SNELL
#####################################################################
vies.beta.cs <- res.bias[1:p]
vies.nu.cs <- res.bias[-(1:p)]
beta.est.cs = beta.est.mv - vies.beta.cs
nu.est.cs = nu.est.mv - vies.nu.cs#####################################################################
# 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)]#####################################################################
# PARAMETRIC BOOTSTRAPPING
#####################################################################
results_bootp = data.frame(bp.fun(Dobs, est.mv))
est.bp = apply(results_bootp, 2, mean)
beta.est.boot = 2*beta.est.mv - est.bp[1:p]
nu.est.boot = 2*nu.est.mv - est.bp[-(1:p)]
##########################################################################################################################################
# RESULTS
#####################################################################
result_est = c(beta.est.mv, nu.est.mv, beta.est.cs, nu.est.cs,
beta.est.df, nu.est.df, beta.est.boot, nu.est.boot)
result_est#####################################################################
# PARAMETROS ESTIMADOS
#####################################################################
nEst = matrix(round(c(result_est),4), ncol=4)
colnames(nEst) = c("MLE", "Cox-Snell", "Firth", "p-boot")
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)
Parametros Estimados
####################################################################
# 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)
Errores Estandars