Ü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.

  • Lassen Sie den Code der Funktion mysim() laufen (sie ist dann im workspace unter Functions ersichtlich)

  • 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. Berechnen Sie global und bezüglich sex, age und jeder sex:age-Kombination:
  • Durchschnittliche Anzahlen und Beobachtungszeiten
  • Summe der Anzahlen und der Beobachtungszeiten
  • Durchschnittsrate (MLE). Beachte hierzu: Der MLE für die Rate einer Gruppe ist das mit der Beobachtungszeit gewichtete arithmetische Mittel der Teilraten.
    • Mit \(\hat{\lambda}_i = \frac{Y_i}{t_i}\) haben wir

\[ \hat{\lambda} = \frac{\sum_{i=1}^n Y_i}{\sum_{i=1}^n t_i} = \frac{\sum_{i=1}^n t_i \hat{\lambda}_i}{\sum_{i=1}^n t_i} = \frac{\bar{Y}}{\bar{t}}\neq \frac{1}{n}\sum_{i=1}^n \hat{\lambda}_i \]

Brauche dazu

  • aggregate(cbind(count,time)~...,FUN=mean) mit ~1, ~sex, ~age und ~sex:age.

  • aggregate(cbind(count,time)~...,FUN=sum), analog.

  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, Beobachtungszeiten und Raten
  • global
aggregate(cbind(count,time)~1,FUN=mean,data=d.sim)#hatten wir schon aus summary()
##   count     time
## 1 10.18 2.722518
aggregate(cbind(count,time)~1,FUN=sum,data=d.sim)
##   count     time
## 1  5090 1361.259
sumsG<-aggregate(cbind(count,time)~1,FUN=sum,data=d.sim)
rateG<-sumsG$count/sumsG$time
cbind(sumsG,rateG)
##   count     time    rateG
## 1  5090 1361.259 3.739186
  • sex
aggregate(cbind(count,time)~sex,FUN=mean,data=d.sim)
##   sex     count     time
## 1   w  8.003831 2.771476
## 2   m 12.556485 2.669053
aggregate(cbind(count,time)~sex,FUN=sum,data=d.sim)
##   sex count     time
## 1   w  2089 723.3552
## 2   m  3001 637.9038
sumsS<-aggregate(cbind(count,time)~sex,FUN=sum,data=d.sim)
rateS<-sumsS$count/sumsS$time
cbind(sumsS,rateS)
##   sex count     time    rateS
## 1   w  2089 723.3552 2.887931
## 2   m  3001 637.9038 4.704471
  • age
aggregate(cbind(count,time)~age,FUN=mean,data=d.sim)
##     age     count     time
## 1 young 11.857143 2.717101
## 2   old  8.568627 2.727723
aggregate(cbind(count,time)~age,FUN=sum,data=d.sim)
##     age count     time
## 1 young  2905 665.6897
## 2   old  2185 695.5693
sumsA<-aggregate(cbind(count,time)~age,FUN=sum,data=d.sim)
rateA<-sumsA$count/sumsA$time
cbind(sumsA,rateA)
##     age count     time    rateA
## 1 young  2905 665.6897 4.363895
## 2   old  2185 695.5693 3.141312
  • age:sex
aggregate(cbind(count,time)~sex+age,FUN=mean,data=d.sim)
##   sex   age     count     time
## 1   w young  9.260163 2.784774
## 2   m young 14.475410 2.648872
## 3   w   old  6.884058 2.759623
## 4   m   old 10.555556 2.690097
aggregate(cbind(count,time)~sex+age,FUN=sum,data=d.sim)
##   sex   age count     time
## 1   w young  1139 342.5272
## 2   m young  1766 323.1624
## 3   w   old   950 380.8279
## 4   m   old  1235 314.7413
sumsAS<-aggregate(cbind(count,time)~sex+age,FUN=sum,data=d.sim)
rateAS<-sumsAS$count/sumsAS$time
cbind(sumsAS,rateAS)
##   sex   age count     time   rateAS
## 1   w young  1139 342.5272 3.325283
## 2   m young  1766 323.1624 5.464744
## 3   w   old   950 380.8279 2.494565
## 4   m   old  1235 314.7413 3.923857
  1. Beobachtete Rate Ratios
  • marginal
rateS[2]/rateS[1] #m vs w
## [1] 1.629011
  1. Rate Ratios bezüglich sex für beide Gruppen von age, bedingt
rateAS[2]/rateAS[1] #in young
## [1] 1.643392
rateAS[4]/rateAS[3] #in old
## [1] 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) (\(\beta_0=1.2\), \(\beta_1=0.5\), \(\beta_2=-0.3\) und \(\beta_3=0\)). Sind unsere Schätzungen verzerrt? Dazu ziehen wir \(k=1000\) mal eine neue Stichprobe, 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 (und Kontraste) bekommen wir mit

  • emmeans(m.poisM,revpairwise~age+sex,type="response",infer=TRUE,offset=0) und
  • emmeans(m.poisM,revpairwise~age+sex,infer=TRUE,offset=0).
  • Machen Sie eine Graphik mit emmeans::emmip(age~sex,...) 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. Das muss hier natürlich so sein, da wir aus einer bekannten Wahrheit simuliert haben. 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
## 1       497     538.28            
## 2       496     537.70  1  0.58022
  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
  • Es gibt keine Evidenz für übergrosse Streuung, 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üssen wir offset=0 setzen (was dann einer Zeiteinheit entspricht, da \(\log(1)=0\)). Auf der Skala des linearen Prädiktors haben wir

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

Damit haben wir für die Raten (hier mit kleinen Rundungsfehlern):

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

Wenn wir vor dem ~ noch pairwise oder revpairwise setzen, haben wir zusätzlich alle Kontraste. Mit type="response" kontrollieren wir obige Rechnungen.

emmeans(m.poisM,revpairwise~age+sex,type="response",infer=TRUE,offset=0)
## $emmeans
##  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.1180 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 
## 
## $contrasts
##  contrast          ratio     SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  old w / young w   0.731 0.0207 Inf      0.68     0.786    1 -11.049  <.0001
##  young m / young w 1.612 0.0460 Inf      1.50     1.735    1  16.755  <.0001
##  young m / old w   2.205 0.0872 Inf      1.99     2.441    1  19.998  <.0001
##  old m / young w   1.179 0.0481 Inf      1.06     1.309    1   4.031  0.0003
##  old m / old w     1.612 0.0460 Inf      1.50     1.735    1  16.755  <.0001
##  old m / young m   0.731 0.0207 Inf      0.68     0.786    1 -11.049  <.0001
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 4 estimates 
## Intervals are back-transformed from the log scale 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Wir sehen unter $contrasts, dass die Effekte von sex unabhängig sind vom level von age und umgekehrt. Im Interaktionsmodell wäre das nicht so.

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,CIs=TRUE,offset=0) # log Rates
emmip(m.poisM,age~sex,CIs=TRUE,type="response",offset=0) # Rates