Générer des échantillons de données à partir de la régression de Poisson

13

Je me demandais comment générer des données à partir d'une équation de régression de Poisson dans R? Je suis un peu confus sur la façon d'aborder le problème.

Donc, si je suppose que nous avons deux prédicteurs et qui sont distribués . Et l'ordonnée à l'origine est 0 et les deux coefficients sont égaux à 1. Alors mon estimation est simplement:X1X2N(0,1)

Journal(Oui)=0+1X1+1X2

Mais une fois que j'ai calculé le log (Y) - comment générer le nombre de poissons en fonction de cela? Quel est le paramètre de taux pour la distribution de Poisson?

Si quelqu'un pouvait écrire un bref script R qui génère des échantillons de régression de Poisson, ce serait génial!

Michael
la source

Réponses:

24

Le modèle de régression du poisson suppose une distribution de Poisson pour et utilise la fonction de liaison . Donc, pour une seule variable explicative , on suppose que (de sorte que ) et que . La génération de données selon ce modèle suit facilement. Voici un exemple que vous pouvez adapter selon votre propre scénario.OuiJournalXOuiP(μ)E(Oui)=V(Oui)=μJournal(μ)=β0+β1X

>   #sample size
> n <- 10
>   #regression coefficients
> beta0 <- 1
> beta1 <- 0.2
>   #generate covariate values
> x <- runif(n=n, min=0, max=1.5)
>   #compute mu's
> mu <- exp(beta0 + beta1 * x)
>   #generate Y-values
> y <- rpois(n=n, lambda=mu)
>   #data set
> data <- data.frame(y=y, x=x)
> data
   y         x
1  4 1.2575652
2  3 0.9213477
3  3 0.8093336
4  4 0.6234518
5  4 0.8801471
6  8 1.2961688
7  2 0.1676094
8  2 1.1278965
9  1 1.1642033
10 4 0.2830910
ocram
la source
3

Si vous vouliez générer un ensemble de données parfaitement adapté au modèle, vous pourriez faire quelque chose comme ceci dans R:

# y <- exp(B0 + B1 * x1 + B2 * x2)

set.seed(1234)

B0 <-  1.2                # intercept
B1 <-  1.5                # slope for x1
B2 <- -0.5                # slope for x2

y <- rpois(100, 6.5)

x2 <- seq(-0.5, 0.5,,length(y))
x1 <- (log(y) - B0 - B2 * x2) / B1

my.model <- glm(y ~ x1 + x2, family = poisson(link = log))
summary(my.model)

Qui retourne:

Call:
glm(formula = y ~ x1 + x2, family = poisson(link = log))

Deviance Residuals: 
       Min          1Q      Median          3Q         Max  
-2.581e-08  -1.490e-08   0.000e+00   0.000e+00   4.215e-08  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.20000    0.08386  14.309  < 2e-16 ***
x1           1.50000    0.16839   8.908  < 2e-16 ***
x2          -0.50000    0.14957  -3.343 0.000829 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 8.8619e+01  on 99  degrees of freedom
Residual deviance: 1.1102e-14  on 97  degrees of freedom
AIC: 362.47

Number of Fisher Scoring iterations: 3
Mark Miller
la source