Comment obtenir la prédiction d'une variable spécifique dans WinBUGS?

10

Je suis un nouvel utilisateur de WinBUGS et j'ai une question pour votre aide. Après avoir exécuté le code suivant, j'ai obtenu les paramètres de beta0through beta4(statistiques, densité), mais je ne sais pas comment obtenir la prédiction de la dernière valeur de h, que j'ai définie NApour modéliser dans le code.

Quelqu'un peut-il me donner un indice? Tout avis serait grandement apprécié.

model {
for(i in 1: N) {
CF01[i] ~ dnorm(0, 20)
CF02[i]  ~ dnorm(0, 1)
h[i] ~ dpois (lambda [i])
log(lambda [i]) <- beta0 + beta1*CF03[i] + beta2*CF02[i] + beta3*CF01[i] + beta4*IND[i]
}
beta0 ~ dnorm(0.0, 1.0E-6)
beta1 ~ dnorm(0.0, 1.0E-6)
beta2 ~ dnorm(0.0, 1.0E-6)
beta3 ~ dnorm(0.0, 1.0E-6)
beta4  <- log(p)
p ~ dunif(lower, upper)
}

INITS
list(beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0, p = 0.9)

DATA(LIST)
list(N = 154, lower = 0.80, upper = 0.95,

h = c(1, 4, 1, 2, 1, 2, 1, 1, 1, 3, 3, 0, 0, 0, 2, 0, 1, 0, 4, 2,
3, 0, 2, 1, 1, 2, 2, 2, 3, 4, 2, 3, 1, 0, 1, 3, 3, 3, 1, 0, 1,
0, 5, 2, 1, 2, 1, 3, 3, 1, 1, 0, 2, 2, 0, 3, 0, 0, 3, 2, 2, 2,
1, 0, 3, 3, 1, 1, 1, 2, 1, 0, 1, 2, 1, 2, 0, 2, 1, 0, 0, 2, 5,
0, 2, 1, 0, 2, 1, 2, 2, 2, 0, 3, 2, 1, 3, 3, 3, 3, 0, 1, 3, 3,
3, 1, 0, 0, 1, 2, 1, 0, 1, 4, 1, 1, 1, 1, 2, 1, 3, 0, 0, 1, 1,
1, 1, 0, 2, 1, 0, 0, 1, 1, 5, 1, 1, 1, 3, 0, 1, 1, 1, 0, 2, 1,
0, 3, 3, 0, 0, 1, 2, 6, NA),

CF03 = c(-1.575, 0.170, -1.040, -0.010, -0.750,
0.665, -0.250, 0.145, -0.345, -1.915, -1.515,
0.215, -1.040, -0.035, 0.805, -0.860, -1.775,
1.725, -1.345, 1.055, -1.935, -0.160, -0.075,
-1.305, 1.175, 0.130, -1.025, -0.630, 0.065,
-0.665, 0.415, -0.660, -1.145, 0.165, 0.955,
-0.920, 0.250, -0.365, 0.750, 0.045, -2.760,
-0.520, -0.095, 0.700, 0.155, -0.580, -0.970,
-0.685, -0.640, -0.900, -0.250, -1.355, -1.330,
0.440, -1.505, -1.715, -0.330, 1.375, -1.135,
-1.285, 0.605, 0.360, 0.705, 1.380, -2.385, -1.875,
-0.390, 0.770, 1.605, -0.430, -1.120, 1.575, 0.440,
-1.320, -0.540, -1.490, -1.815, -2.395, 0.305,
0.735, -0.790, -1.070, -1.085, -0.540, -0.935,
-0.790, 1.400, 0.310, -1.150, -0.725, -0.150,
-0.640, 2.040, -1.180, -0.235, -0.070, -0.500,
-0.750, -1.450, -0.235, -1.635, -0.460, -1.855,
-0.925, 0.075, 2.900, -0.820, -0.170, -0.355,
-0.170, 0.595, 0.655, 0.070, 0.330, 0.395, 1.165,
0.750, -0.275, -0.700, 0.880, -0.970, 1.155, 0.600,
-0.075, -1.120, 1.480, -1.255, 0.255, 0.725,
-1.230, -0.760, -0.380, -0.015, -1.005, -1.605,
0.435, -0.695, -1.995, 0.315, -0.385, -0.175,
-0.470, -1.215, 0.780, -1.860, -0.035, -2.700,
-1.055, 1.210, 0.600, -0.710, 0.425, 0.155, -0.525,
-0.565),

CF02 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, 0.38, 0.06, -0.94,
-0.02, -0.28, -0.78, -0.95, 2.33, 1.43, 1.24, 1.26,
-0.75, -1.5, -2.09, 1.01, -0.05, 2.48, 2.48, 0.46,
0.46, -0.2, -1.11, 0.52, -0.37, 0.58, 0.86, 0.59,
-0.12, -1.33, 1.4, -1.84, -1.4, -0.76, -0.23,
-1.78, -1.43, 1.2, 0.32, 1.87, 0.43, -1.71, -0.54,
-1.25, -1.01, -1.98, 0.52, -1.07, -0.44, -0.24,
-1.31, -2.14, -0.43, 2.47, -0.09, -1.32, -0.3,
-0.99, 1.1, 0.41, 1.01, -0.19, 0.45, -0.07, -1.41,
0.87, 0.68, 1.61, 0.36, -1.06, -0.44, -0.16, 0.72,
-0.69, -0.94, 0.11, 1.25, 0.33, -0.05, 0.87, -0.37,
-0.2, -2.22, 0.26, -0.53, -1.59, 0.04, 0.16, -2.66,
-0.21, -0.92, 0.25, -1.36, -1.62, 0.61, -0.2, 0,
1.14, 0.27, -0.64, 2.29, -0.56, -0.59, 0.44, -0.05,
0.56, 0.71, 0.32, -0.38, 0.01, -1.62, 1.74, 0.27, 0.97,
1.22, -0.21, -0.05, 1.15, 1.49, -0.15, 0.05, -0.87,
-0.3, -0.08, 0.5, 0.84, -1.67, 0.69, 0.47, 0.44,
-1.35, -0.24, -1.5, -1.32, -0.08, 0.76, -0.57,
-0.84, -1.11, 1.94, -0.68),

CF01 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, -0.117, -0.211, -0.333, -0.229, -0.272,
-0.243, -0.148, 0.191, -0.263, -0.239, -0.168,
-0.381, -0.512, -0.338, -0.296, 0.067, 0.104,
-0.254, -0.167, -0.526, -0.096, -0.43, 0.013,
-0.438, -0.297, -0.131, -0.098, -0.046, -0.063,
-0.194, -0.155, -0.645, -0.603, -0.374, -0.214,
-0.165, -0.509, -0.171, -0.442, -0.468, -0.289,
-0.427, -0.519, -0.454, 0.046, -0.275, -0.401,
-0.542, -0.488, -0.52, -0.018, -0.551, -0.444,
-0.254, -0.286, 0.048, -0.03, -0.015, -0.219,
-0.029, 0.059, 0.007, 0.157, 0.141, -0.035, 0.136,
0.526, 0.113, 0.22, -0.022, -0.173, 0.021, -0.027,
0.261, 0.082, -0.266, -0.284, -0.097, 0.097, -0.06,
0.397, 0.315, 0.302, -0.026, 0.268, -0.111, 0.084,
0.14, -0.073, 0.287, 0.061, 0.035, -0.022, -0.091,
-0.22, -0.021, -0.17, -0.184, 0.121, -0.192,
-0.24, -0.283, -0.003, -0.45, -0.138, -0.143,
0.017, -0.245, 0.003, 0.108, 0.015, -0.219, 0.09,
-0.22, -0.004, -0.178, 0.396, 0.204, 0.342, 0.079,
-0.034, -0.122, -0.24, -0.125, 0.382, 0.072, 0.294,
0.577, 0.4, 0.213, 0.359, 0.074, 0.388, 0.253, 0.167),

IND = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
bayesian bugs prediction winbugs Bo Yu
la source
Ne demandez-vous pas simplement la valeur de lambda [N]?
whuber
@whuber oui, je pense que c'est correct, mais plus fondamentalement, vous devez avoir des choses que vous voulez prédire (c'est-à-dire, avoir une distribution postérieure pour) être distinctes des choses que vous avez déjà observées. Vous pouvez soit faire la prédiction explicitement dans les winbugs ou dans le post-traitement en utilisant les échantillons des bêtas.
atiretoo - réintègre monica le
@atiretoo Pour autant que je sache, les lambdas sont exactement ce que l'on veut prédire: il s'agit d'un modèle linéaire généralisé pour une distribution de Poisson avec loglog et les lambdas sont les paramètres de Poisson prédits. Ils n'ont pas été observés. Je crois que tout ce qu'il faut faire ici est de régler un moniteur sur lambda [N].
whuber
@whuber, je dirais plutôt surveiller h[N]plutôt que lambda[N]... et vous obtenez la distribution postérieure de la valeur prédite.
Curieux
@tomek mais h[N]n'est pas la valeur prédite: ce sera une collection de tirages à partir d'un ensemble de distributions de Poisson prédites. En tant que tel, il combine la variation des paramètres de Poisson et la variation de ces distributions de Poisson elles-mêmes. Ce qui est pertinent, c'est la distribution postérieure de lambda[N].
whuber

Réponses:

6

Ajoutez simplement la variable hà la liste des paramètres à surveiller. Si vous utilisez un package comme R2WinBUGS, ajoutez une variable hà la liste passée à l' parameters.to.saveargument de la bugsfonction. Ensuite, regardez votre dernière valeur dans h(celle avec NA) - vous obtiendrez une distribution postérieure là-bas.

C'est la manière habituelle de faire des prédictions dans l'inférence bayésienne ( voir aussi cette question ). C'est sympa et simple! Plus de séparation de l'évaluation et de la prévision des paramètres. Tout se fait en même temps. La distrubution postérieure des paramètres est donnée par les données réelles et propagée aux valeurs NA (sous forme de "prédictions").

Curieuse
la source
Tomas, merci pour votre aide. J'essaie de surveiller la variable de h dans Sample Monitor Tool mais cela ne fonctionne pas. Pourriez-vous m'aider à nouveau? Voici la procédure que j'ai effectuée dans WinBUGS (je ne sais pas comment utiliser R2WinBUGS): 1) sélectionnez Sample dans Sample Monitor Tool 2) tapez h dans la case blanche marquée node 3) cliquez sur le bouton marqué set 4) h is pas sur la liste des paramètres que je veux surveiller, tandis que d'autres paramètres (beta0, beta1, beta2, beta3, an p) sont affichés dans la liste. Savez-vous comment ajouter "h" à la liste des paramètres que je souhaite surveiller? Merci encore!
Bo Yu
@BoYu, je ne sais pas comment le faire directement dans WinBUGS puisque j'exécute WinBUGS à partir de R, en utilisant le package R2WinBUGS. C'est beaucoup plus pratique car vous pouvez simplement enregistrer le script R et l'exécuter en tant que lot, tout en produisant vos propres graphiques, etc. Regardez ici des exemples de scripts.
Curieux
Cela dit, cela sera certainement possible dans WinBUGS lui-même, mais je ne sais pas comment (et je suppose que la plupart des gens l'appellent depuis R).
Curieux
1
Tout d'abord, merci, whuber, atiretoo et Tomas! Comme déjà mentionné, oui, c'est un modèle linéaire généralisé, la variable de h est ajustée par la distribution de Poisson à taux variable (lambda) conditionné avec différents prédicteurs (CF01, CF02, CF03 et IND). La dernière valeur de h est ce que je dois savoir et n'est pas observée (marquée comme NA), tandis que toutes les autres valeurs de h sont observées. Je pense que whuber a raison, je dois définir lambda comme paramètre dans Sample Monitor Tool et vérifier les statistiques de la dernière valeur de lambda, et obtenir en outre quelle est ma prédiction de la dernière heure. Merci a tous.
Bo Yu
1
@ Tomas, merci beaucoup. Oui, tu as raison! WinBUGS fournit la prédiction de h [N], y compris les statistiques et la densité de probabilité. Je l'ai maintenant. Cordialement,
Bo Yu