Übung zur Poisson-Regression

Simulation und Analyse von Poisson-Daten

Das Ziel dieser Übung ist das Verstehen von Poisson-Daten und einer Poisson-Regression. Wir werden hier Daten gemäss einem als bekannt vorausgesetzten Modell (zweifaktorielles Design) simulierten und dann analysieren.

Zusätzlich üben wir den Umgang mit Interaktionseffekten. In einem zweifaktoriellen Design haben Faktoren \(A\) und \(B\) einen Interaktionseffekt, wenn die erwartete Veränderung in der Zielgrösse bei einer Veränderung der Stufe auf dem Faktor \(A\) von der Stufe des Faktors \(B\) abhängt, und umgekehrt. Wenn jede Veränderung der Stufe des Faktors \(A\) über alle Stufen vom Faktor \(B\) konstant ist, dann haben die beiden Faktoren keinen Interaktionseffekt, sie agieren dann \(additiv\) und das Modell hat nur Haupteffekte.

Simuliere mit untigem Code Poisson-Daten (\(n=500\)) aus dem folgendem Modell:

\[ \boxed{\log \mu_i=\beta_0+\beta_1 sex_i+\beta_2 age_i+\beta_3 (sex:age)_i+\log(time_i), \qquad Y_i\sim \mathcal{P}(\mu_i)} \]

  • mit den Indikator-Variablen \[ sex_{i}\begin{cases} 1 &\mbox{if }sex_{i}=m\\ 0 & \mbox{if } sex_{i}=w \end{cases}, \quad age_{i}\begin{cases} 1 &\mbox{if }age_{i}=old\\ 0 & \mbox{if } age_{i}=young \end{cases}, \quad sex:age_i\begin{cases} 1 &\mbox{if }sex_{i}=m,age_i=old\\ 0 & else \end{cases} \]

  • sex: biologisches Geschlecht, kategorielle Eingangsgrösse

  • age: Altersgrupppe, kategorielle Eingangsgrösse

  • time: Zeitdauer der Beobachtung

  • Count: Anzahl Ereignisse

    Simulation bedeutet, dass wir allen Parametern einen wahren Wert zuordnen. Sie müssen folgenden Code natürlich nicht selber schreiben können, es geht hier darum, zu verstehen, wie Daten gemäss einem Modell “entstehen”. Gleichzeitig können Sie das Verständnis von vielen Aspekten wiederholen, die wir in diesem Modul eingeführt haben (Parameter, lineare Prädiktor, Link-Funktion, Verteilung usw.)

    Versuchen Sie also, die Schritte der Simulation nachzuvollziehen.

    Speichern Sie die Funktion mysim() ab (sie ist dann im workspace unter Functions ersichtlich) und definieren Sie dann das data frame d.sim mit d.sim<-mysim() (ohne Argument). Betrachten Sie die Struktur der Daten.

mysim<-function(seed=123) {
  set.seed(seed) # for reproducibility
  n<-500 # Number of observations
  # Parameters                         
  beta0 <- 1.2                      # Intercept
  beta1 <- 0.5                      # main sex effect (m vs. w)
  beta2 <- -0.3                     # main age effect (old vs. young)
  beta3 <- 0                        # interaction effect
  beta<-c(beta0,beta1,beta2,beta3)
  # Simulate grouping variables
  sex <- factor(sample(c("w", "m"), size=n,replace=TRUE),levels=c("w","m"))  # First grouping variable
  age <- factor(sample(c("young", "old"), size=n,replace = TRUE),levels=c("young","old"))  # Second grouping variable
  # Create the design matrix for modelling
  G1_num <- ifelse(sex == "m", 1, 0)  # Convert to numeric for modeling
  G2_num <- ifelse(age == "old", 1, 0)
  G1G2_num <- G1_num * G2_num
  # Simulate Exposure times
  time <- runif(n, min = 0.5, max = 5)  # Exposure times from 0.5 to 5
  # Compute the linear predictor
  linear_predictor <- beta0 + beta1 * G1_num + beta2 * G2_num + beta3 * G1G2_num + log(time)
  # Compute the expectation
  mu<-exp(linear_predictor)
  # Simulate Poisson response variable
  count <- rpois(n, lambda = mu)
  # Combine into a data frame
  data.frame(count = count,sex = sex,age = age,time = time)
}
d.sim<-mysim()
str(d.sim)
## 'data.frame':    500 obs. of  4 variables:
##  $ count: int  8 15 3 14 15 13 21 13 2 8 ...
##  $ sex  : Factor w/ 2 levels "w","m": 1 1 1 2 1 2 2 2 1 1 ...
##  $ age  : Factor w/ 2 levels "young","old": 2 2 2 2 1 1 1 1 2 1 ...
##  $ time : num  1.73 3.17 1.22 4.34 4.31 ...
psych::headTail(d.sim,10,10)
##     count  sex   age time
## 1       8    w   old 1.73
## 2      15    w   old 3.17
## 3       3    w   old 1.22
## 4      14    m   old 4.34
## 5      15    w young 4.31
## 6      13    m young 2.65
## 7      21    m young 3.98
## 8      13    m young 1.83
## 9       2    w   old  0.8
## 10      8    w young 2.48
## ...   ... <NA>  <NA>  ...
## 491     0    w young 0.56
## 492     0    m   old 0.55
## 493    13    w young 4.98
## 494    16    m   old 3.95
## 495     8    m young 1.92
## 496    19    m young 4.15
## 497    22    m   old 4.96
## 498     6    w   old  1.4
## 499     8    w   old 3.41
## 500     4    w   old 2.47

Deskriptive Statistik

  1. Beschreiben Sie die Variablen mit summary()
  2. Betrachten Sie eine Kreuztabelle von sex und age. (mit table()). Sind die Daten balanciert, d.h. haben alle Zellen gleich viele Beobachtungen?
  3. Berechne die durchschnittliche Anzahlen und Raten bezüglich sex, age und sex:age-Kombination. Für die Raten berechnen wir harmonische Mittel:
  • sum(d.sim$count)/sum(d.sim$time) statt
  • mean(d.sim$count/d.sim$time) (für overall) und
  • tapply(d.sim,list(d.sim$sex,d.sim$age),FUN=function(x){sum(x$count)/sum(x$time)}) statt
  • tapply(d.sim$count/d.sim$time,list(d.sim$sex,d.sim$age),mean) (für sex:age-Kombination)

  1. Wie gross ist die beobachtete Rate Ratio bezüglich sex (\(m/w\)), unabhängig von age (“marginal”)?

  2. Wie gross ist die beobachtete Rate Ratio bezüglich sex (\(m/w\)), in beiden Gruppen von age? (“conditional”). Ist die Rate Ratio klar verschieden in beiden age-Gruppen?

  1. Beschreiben

summary()

summary(d.sim)
##      count       sex        age           time       
##  Min.   : 0.00   w:261   young:245   Min.   :0.5052  
##  1st Qu.: 5.00   m:239   old  :255   1st Qu.:1.5992  
##  Median : 9.00                       Median :2.7398  
##  Mean   :10.18                       Mean   :2.7225  
##  3rd Qu.:14.00                       3rd Qu.:3.8133  
##  Max.   :38.00                       Max.   :4.9979
  1. Kreuztabelle
table(d.sim$sex,d.sim$age)
##    
##     young old
##   w   123 138
##   m   122 117

Die Daten sind nicht ganz balanciert.

  1. Anzahlen und Raten
  • Anzahlen
mean(d.sim$count)
## [1] 10.18
tapply(d.sim$count,d.sim$sex,mean)
##         w         m 
##  8.003831 12.556485
tapply(d.sim$count,d.sim$age,mean)
##     young       old 
## 11.857143  8.568627
tapply(d.sim$count,list(d.sim$sex,d.sim$age),mean)
##       young       old
## w  9.260163  6.884058
## m 14.475410 10.555556
  • Raten (Harmonic means!)
ratesOverall<-sum(d.sim$count)/sum(d.sim$time)
ratesAge<-tapply(d.sim,list(d.sim$age),FUN=function(x){sum(x$count)/sum(x$time)})
ratesSex<-tapply(d.sim,list(d.sim$sex),FUN=function(x){sum(x$count)/sum(x$time)})
ratesSexAge<-tapply(d.sim,list(d.sim$sex,d.sim$age),FUN=function(x){sum(x$count)/sum(x$time)})
ratesOverall
## [1] 3.739186
ratesAge
##    young      old 
## 4.363895 3.141312
ratesSex
##        w        m 
## 2.887931 4.704471
ratesSexAge
##      young      old
## w 3.325283 2.494565
## m 5.464744 3.923857
  1. Beobachtete Rate Ratios
  • marginal
ratesSex[2]/ratesSex[1]
##        m 
## 1.629011
  1. Rate Ratios bezüglich sex für beide Gruppen von age, bedingt
ratesSexAge
##      young      old
## w 3.325283 2.494565
## m 5.464744 3.923857
apply(ratesSexAge,2, function(x)x[2]/x[1])
##    young      old 
## 1.643392 1.572962

Die Rate Ratio bezüglich sex scheint nur leicht vom level auf age abzuhängen. Wir vermuten hier schon, dass es keinen Interaktionseffekt gibt.

Modelle

  1. Passen Sie ein glm an mit Interaktionseffekt (mit Namen mPoisI) (das heisst, mit Termen sex+age+sex:age oder kurz sex*age) und mit dem Offset-Term. Passen Sie gleichzeitig auch ein Modell an ohne Interaktion (mit Namen mPoisM) und ein Modell ohne sex und age (mit Namen mPois0) (nur mit Intercept und Offset.)

  2. Fakultativ: Vergleichen Sie die geschätzten Koeffizienten von mPoisI mit den wahren Parametern (die wir natürlich nur hier kennen, weil wir simuliert haben). Sind unsere Schätzungen verzerrt? Dazu ziehen wir \(k=1000\) mal eine neue Stichprobe und passen das Modell an (mit untigem Code). Studieren Sie die Matrix betasim (Beschreiben Sie die vier Kolonnen). Verzerrt ist eine Schätzung, wenn sie im Schnitt nicht gleich dem wahren Wert ist.

k<-1000
betasim<-matrix(NA,nrow=k,ncol=4) #Matrix zum Auffüllen 
for (i in 1:k) {
  betasim[i,]<-glm(count~sex*age+offset(log(time)),family=poisson(link="log"), data=mysim(seed=i))$coef
}
##
head(betasim)
##          [,1]      [,2]       [,3]         [,4]
## [1,] 1.189445 0.5506043 -0.3096503 -0.084067866
## [2,] 1.238503 0.4839592 -0.3221731 -0.045097397
## [3,] 1.242606 0.4548178 -0.3334876  0.003967023
## [4,] 1.207219 0.4919022 -0.3138161  0.070311851
## [5,] 1.140620 0.5888954 -0.1601007 -0.209345945
## [6,] 1.215929 0.4800045 -0.2435946 -0.006623491

  1. Wenden Sie drop1 an, um den sex:age Interaktionseffekt zu testen. Kann man ihn weglassen?

  2. Betrachten Sie das Modell ohne Interaktion (Haupteffektmodell). Berechnen Sie das Verhältnis der Residualdevianz zu den Residualfreiheitsgraden. Ist dieses klar grösser als 1? (was übergrosser Streuung entsprechen würde).

  3. Berechnen Sie eine Punkt- und Intervallschätzung (im Haupteffektmodell) für den sex und age-Effekt auf Rate Ratio-Skala.

  4. Wie gross ist die geschätzte relative Ratenreduktion von alt versus jung (im Haupteffektmodell)?

  5. Um welchen Faktor ist die geschätzte Rate gemäss Modell bei sex=w grösser als bei sex=m (im Haupteffektmodell)?

  6. Wie gross ist die geschätzte Rate gemäss Modell für eine weiblich junge Person (im Haupteffektmodell)?

  7. Die geschätzten Raten und log Raten gemäss Modell bekommen wir mit

  • emmeans(m.poisM,~age+sex,type="response",infer=TRUE,offset=0) und
  • emmeans(m.poisM,~age+sex,infer=TRUE,offset=0).
  • Machen Sie eine Graphik mit emmeans::emmip() zum Visualisieren der vom Modell vorhergesagten Raten für die vier sex:age Gruppen, einmal auf dem linearen Prädiktor und einmal bezüglich der Raten (type="response").
  • Setzen Sie immer offset=0 (das entspricht ja einer Zeiteinheit), damit man Raten statt Anzahlen hat.
  1. Wir passen ein Poisson Modell an mit glm() mit den zwei Haupteffekte sex und age und dem sex:age-Interaktionseffekt (sex*age ist äquivalent mit sex+age+sex:age). Gleichzeitig passen wir noch ein Modell an ohne Interaktionseffekt und ein Modell ohne Eingangsgrössen.
m.poisI<-glm(count~sex*age+offset(log(time)),family=poisson(link="log"), data=d.sim)
m.poisM<-update(m.poisI,.~.-sex:age)
m.pois0<-update(m.poisI,.~.-sex*age)
coef(m.poisI)
## (Intercept)        sexm      ageold sexm:ageold 
##  1.20155480  0.49676249 -0.28744050 -0.04380172
coef(m.poisM)
## (Intercept)        sexm      ageold 
##   1.2131174   0.4776709  -0.3130446
coef(m.pois0)
## (Intercept) 
##    1.318868
summary(m.poisI)
## 
## Call:
## glm(formula = count ~ sex * age + offset(log(time)), family = poisson(link = "log"), 
##     data = d.sim)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.20155    0.02963  40.551  < 2e-16 ***
## sexm         0.49676    0.03800  13.072  < 2e-16 ***
## ageold      -0.28744    0.04394  -6.542 6.07e-11 ***
## sexm:ageold -0.04380    0.05750  -0.762    0.446    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 960.71  on 499  degrees of freedom
## Residual deviance: 537.70  on 496  degrees of freedom
## AIC: 2505.1
## 
## Number of Fisher Scoring iterations: 4
  1. Fakultativ: Die Schätzungen sind hier unverzerrt (d.h. nicht systematisch verschieden von den wahren Werten), weil wir aus der Wahrheit simuliert haben. Zufallsstreuung bleibt natürlich immer.
  • Kennwerte der Verteilung von \(k=1000\) simulierten \(\hat{\beta}_0\),\(\hat{\beta}_1\),\(\hat{\beta}_2\),\(\hat{\beta}_3\):
psych::describe(betasim)
##    vars    n mean   sd median trimmed  mad   min   max range  skew kurtosis se
## X1    1 1000  1.2 0.03    1.2     1.2 0.03  1.11  1.28  0.17 -0.05    -0.12  0
## X2    2 1000  0.5 0.04    0.5     0.5 0.04  0.37  0.62  0.24 -0.05     0.18  0
## X3    3 1000 -0.3 0.05   -0.3    -0.3 0.05 -0.46 -0.13  0.33  0.02     0.14  0
## X4    4 1000  0.0 0.06    0.0     0.0 0.06 -0.21  0.18  0.39 -0.03     0.02  0

Die mean sind praktisch identisch mit den wahren \(\beta\)’s. Die Schätzungen sind also unverzerrt. sd sind die Standardfehler. Vergleiche mit den Standardfehlern unter Std.Error.

  1. Das reduzierte Modell ohne Interaktionseffekt kann nicht verworfen werden:
drop1(m.poisI,test="Chisq")
## Single term deletions
## 
## Model:
## count ~ sex * age + offset(log(time))
##         Df Deviance    AIC     LRT Pr(>Chi)
## <none>       537.70 2505.1                 
## sex:age  1   538.28 2503.6 0.58022   0.4462
anova(m.poisM,m.poisI)# $äquivalent
## Analysis of Deviance Table
## 
## Model 1: count ~ sex + age + offset(log(time))
## Model 2: count ~ sex * age + offset(log(time))
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       497     538.28                     
## 2       496     537.70  1  0.58022   0.4462
  1. Wir betrachten das Haupteffekt-Modell
summary(m.poisM)
## 
## Call:
## glm(formula = count ~ sex + age + offset(log(time)), family = poisson(link = "log"), 
##     data = d.sim)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.21312    0.02530   47.95   <2e-16 ***
## sexm         0.47767    0.02851   16.75   <2e-16 ***
## ageold      -0.31304    0.02833  -11.05   <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: 960.71  on 499  degrees of freedom
## Residual deviance: 538.28  on 497  degrees of freedom
## AIC: 2503.6
## 
## Number of Fisher Scoring iterations: 4
  • Übergrosse Streuung haben wir nicht, da das Verhältnis von Residualdevianz zu den Residualfreiheitsgraden nahe bei 1 ist.
deviance(m.poisM)/df.residual(m.poisM)
## [1] 1.08305
  1. Punkt und Intervallschätzung der Rate-Ratios für sex und age
  • Punktschätzung und Likelihood Konfidenzintervalle mit confint():
exp(m.poisM$coef)##Punktschätzung
## (Intercept)        sexm      ageold 
##   3.3639552   1.6123148   0.7312173
exp(confint(m.poisM))##Intervallschätzung
##                 2.5 %    97.5 %
## (Intercept) 3.2001898 3.5338550
## sexm        1.5248312 1.7051354
## ageold      0.6916662 0.7729176
  • Alternative: parameters() mit exponentiate=TRUE gibt ebenfalls Likelihood-CI’s für Rate Ratios.

\(RR<1\) entsprechen negativen Effekten auf log Skala, \(RR>1\) entsprechen positiven Effekten auf log Skala.

parameters::parameters(m.poisM,exponentiate=TRUE)
## Parameter   |  IRR |   SE |       95% CI |      z |      p
## ----------------------------------------------------------
## (Intercept) | 3.36 | 0.09 | [3.20, 3.53] |  47.95 | < .001
## sex [m]     | 1.61 | 0.05 | [1.52, 1.71] |  16.76 | < .001
## age [old]   | 0.73 | 0.02 | [0.69, 0.77] | -11.05 | < .001
parameters::parameters(m.poisM)
## Parameter   | Log-Mean |   SE |         95% CI |      z |      p
## ----------------------------------------------------------------
## (Intercept) |     1.21 | 0.03 | [ 1.16,  1.26] |  47.95 | < .001
## sex [m]     |     0.48 | 0.03 | [ 0.42,  0.53] |  16.76 | < .001
## age [old]   |    -0.31 | 0.03 | [-0.37, -0.26] | -11.05 | < .001

  1. Die geschätzte Rate Ratio für age ist 0.7312173, also ist die geschätzte relative Ratenreduktion (alt zu jung) -0.2687827 oder -26.8782737 Prozent.

  2. Die geschätzte Rate Ratio für sex=m versus sex=w ist 1.6123148, also ist der Faktor, um den die Rate für sex=w grösser ist als für sex=m gleich \(\frac{1}{1.6123148}\)=0.6202263

  3. Die geschätzte Event-Rate für eine weibliche junge Person ist 3.3639552 Events pro Zeiteinheit.

  4. Geschätzte Raten und log Raten gemäss Modell: Damit wir Raten statt log Raten haben, müssen wir type="response" setzen. Damit wir Raten statt Anzahlen haben, müssten wir offset=0 setzen (was dann einer Zeiteinheit entspricht, da \(\log(1)=0\)).

library(emmeans)
emmeans(m.poisM,~age+sex,infer=TRUE,offset=0)
##  age   sex emmean     SE  df asymp.LCL asymp.UCL z.ratio p.value
##  young w     1.21 0.0253 Inf     1.164     1.263  47.948  <.0001
##  old   w     0.90 0.0269 Inf     0.847     0.953  33.476  <.0001
##  young m     1.69 0.0217 Inf     1.648     1.733  77.816  <.0001
##  old   m     1.38 0.0246 Inf     1.329     1.426  55.918  <.0001
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
emmeans(m.poisM,~age+sex,type="response",infer=TRUE,offset=0)
##  age   sex rate     SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  young w   3.36 0.0851 Inf      3.20      3.53    1  47.948  <.0001
##  old   w   2.46 0.0661 Inf      2.33      2.59    1  33.476  <.0001
##  young m   5.42 0.1178 Inf      5.20      5.66    1  77.816  <.0001
##  old   m   3.97 0.0977 Inf      3.78      4.16    1  55.918  <.0001
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## Tests are performed on the log scale

Diese Grössen reproduzieren wir auch “von Hand” aus den geschätzten Parametern mit:

3.36            #young woman, oder: exp(1.21)
## [1] 3.36
3.36*0.73       #old woman, oder: exp(1.21-0.31)=exp(1.21)exp(-0.31) 
## [1] 2.4528
3.36*1.61       #young man, oder: exp(1.21+0.48)=exp(1.21)exp(0.48)
## [1] 5.4096
3.36*0.73*1.61  #old man, oder: exp(1.21+0.48-0.31)=exp(1.21)exp(0.48)exp(-0.31)
## [1] 3.949008

Mit emmeans::emmip() können wir Vorhersagen gemäss Modell plotten, auf dem linearen Prädiktor und auf “Response” (Rate).

emmip(m.poisM,age~sex,infer=TRUE,offset=0) # log Rates
emmip(m.poisM,age~sex,infer=TRUE,type="response",offset=0) # Rates