Régression bêta: simulations de Monte Carlo pour les coefficients

Question

J'utilise betaregbetaregresion pour le pourcentage de dégradation modèle au fil du temps: perc_degr~time. Je suit à demi-vie calculate ( en supposant fucntion lien logit),: th = -coef(model)[1]/coef(model)[2]. J'ai également utilisé la mvrnormcommande du MASSpackage pour faire des simulations de Monte Carlo pour avoir une idée de la distribution des valeurs et calculer plus tard les intervalles de confiance:

sims_coef = mvrnorm(n=10000, mu=coefs, Sigma=vcovmat, empirical=T)
th_sims = (-sims_coef[,1])/sims_coef[,2]

Maintenant, j'ai compris que cela betaregme donne le paramètre de précision phi. Comment incorporer cette valeur pour simuler les intervalles de distribution / confiance pour la demi-vie?
Merci beaucoup!

Answer 1

La stratégie consiste à simuler les coefficients à partir de la distribution normale conjointe de tous les paramètres, y compris le$\phi$paramètre. La raison en est que les coefficients pour$\mu$ et $\phi$ne sont pas orthogonales, c'est-à-dire dépendent les uns des autres. Après avoir obtenu le vecteur de paramètres simulé complet, vous pouvez calculer la quantité que vous voulez à partir du sous-ensemble de paramètres pertinent.

Je recommanderais cependant d'utiliser un lien de journal pour le $\phi$dans ce cas, car l'approximation normale sera généralement bien meilleure sur l'échelle logarithmique car les paramètres sont illimités. Le moyen le plus simple d'y parvenir est de spécifier une formule en deux parties: perc_degr ~ time | 1dans votre cas, déclarer que$\mu$dépend de timemais$\phi$est constante. De plus, je vous encourage à envisager un modèle avec des$\phi$cependant: perc_degr ~ time | timeoù les deux paramètres sont autorisés à changer time.

À titre d'exemple concret, considérons l'illustration suivante basée sur un modèle (trop) simple pour les ReadingSkillsdonnées disponibles dans le package:

library("betareg")
data("ReadingSkills", package = "betareg")
m <- betareg(accuracy ~ iq | 1, data = ReadingSkills)

Ensuite, nous pouvons simuler 10000 tirages à partir du vecteur de paramètres complet (en utilisant mvtnormplutôt que MASS):

library("mvtnorm")
set.seed(1)
cf <- rmvnorm(n = 10000, mean = coef(m), sigma = vcov(m))

Et puis nous pouvons calculer n'importe quelle quantité que nous voulons. Par exemple, nous pouvons calculer le rapport négatif des coefficients pour$\mu$ et comparez-la avec la quantité aux coefficients estimés:

hist(-cf[,1]/cf[,2], main = "")
abline(v = -coef(m)[1]/coef(m)[2], col = 2, lwd = 2)

Dans ce modèle particulier, cette intrigue semble plutôt mauvaise car nous devons extrapoler un peu. Par conséquent, je montre également une autre quantité qui semble un peu plus raisonnable sur cet ensemble de données.

Calculons également la valeur prédite $\mu$au régresseur iqdéfini sur -1, 0 et 1. Nous avons donc mis en place une matrice appropriée de 10 000 x 3 et calculons le$\mu$ avec le lien inverse de l'objet modèle ajusté:

x <- matrix(rep(c(-1, 0, 1), each = nrow(cf)), ncol = 3)
mu <- m$link$mean$linkinv(cf[,1] + cf[,2] * x)
mu0 <- predict(m, newdata = data.frame(iq = c(-1, 0, 1)), type = "response")

La quantité ajustée du modèle d'origine est également ajoutée à des fins de comparaison:

for(i in 1:3) {
  hist(mu[, i], main = "", xlab = paste("x =", x[1,i]), xlim = c(0.5, 1))
  abline(v = mu0[i], col = 2, lwd = 2)
}