Kapitel 78 Inzidenz Raten

Dieses Kapitel ist noch nicht fertig und ist unkorrigiert.

78.0.1 Wir benötigen folgende Pakete:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(Statamarkdown)
## Stata found at C:/Program Files/Stata17/StataMP-64.exe
## The 'stata' engine is ready to use.

Wir schauen uns die Inzidenzen mit dem Beispiel von Stürzen an.

Wir simuliern Daten, die an die Sturzverteilung einer Studie angelehnt sind Swiss Chef Trial

set.seed(666)
sample_size=400
number_of_falls<-sample(c(0,1,2,3,4,5,6,8,9,16,18,20), size=sample_size, rep=TRUE, prob=c(0.61, 0.18, 0.1, 0.049, 0.02, 0.017, 0.0049, 0.004733333, 0.004733333, 0.0025, 0.0049,0.002233333))
hist(number_of_falls)

Im ersten Beispiel nehmen wir an, alle Teilnehmern sind genau 365 Tage in der Studie gewesen.

time_at_risk_days=rep(365, times=sample_size)

id<-1:sample_size

data<-data.frame(id, number_of_falls, time_at_risk_days)
rm(time_at_risk_days, id, number_of_falls)

78.1 Berechnung der Inzidenz, Formel 1

\[\begin{equation} \frac{\sum{Stürze Aller Teilnehmenden}}{\sum{ZeitUnterBeobachtung}} \tag{78.1} \end{equation}\]

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 402
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 146000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.002753425
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) # Diese und die nächsten Formeln sind äquivalent, einfach anders ausformuliert
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.005 
cat("Zweite Formulierung, gleiches Resultat: ", (incidence_rate_per_year<- (Total_Number_of_Falls / Total_Exposure_Time) *365))
Zweite Formulierung, gleiches Resultat:  1.005
cat("\nDritte Formulierung, gleiches Resultat: ", (incidence_rate_per_year<- Total_Number_of_Falls / Total_Exposure_Time *365))

Dritte Formulierung, gleiches Resultat:  1.005

78.1.1 Berechnung mit individueller Zeit

Wir können jetzt auch die Berechnungen mit den individuellen Zeiten machen - die jedoch in diesem Beispiel alle gleich sind.

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.005

Die beiden Berechnungsvarianten ergeben genau das gleiche Resultat.

78.2 Im nächsten Beispiel wurden nicht alle Teilnehmenden über 365 Tage verfolgt.

Wir simulieren die Daten zunächst so, dass es keine Korrelation zwischen Anzahl Stürze und Beobachtungszeit gibt (was ja eigentlich nicht ganz realistisch ist). Danke an diese Webseite: https://stats.stackexchange.com/questions/15011/generate-a-random-variable-with-a-defined-correlation-to-an-existing-variables.

n     <- sample_size                    # length of vector
rho   <- 0.0                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho
[1] 1.952382e-17
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  338.4   339.6   340.0   340.0   340.4   341.6 
cor(data$number_of_falls, data$time_at_risk_days)
[1] 1.468732e-15
cor(data$number_of_falls, data$time_at_risk_days)
[1] 1.468732e-15
ggplot(data, aes(x=time_at_risk_days, y=number_of_falls))+
  geom_point()+
  geom_smooth(method="lm")
`geom_smooth()` using formula = 'y ~ x'

78.2.1 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 402
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.002955882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.078897

78.2.2 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.0789

Die beiden Berechnungen stimmen praktisch überein - nicht 100%, da ja die Korrelation nur in der Population 0 ist und nicht in der Stichprobe.

78.2.3 Korrelation zwischen Anzahl Stürze und Beobachtungszeit.

n     <- sample_size                    # length of vector
rho   <- 0.5                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)     
[1] 0.5
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  337.8   339.5   340.0   340.0   340.4   342.9 
cor(data$number_of_falls, data$time_at_risk_days)
[1] 0.5

78.2.4 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 402
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.002955882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.078897

78.2.5 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.076615

Nun sehen wir, dass die Werte nicht mehr exakt übereinstimmen.

78.2.6 negative Korrelation zwischen Anzahl Stürzen und Beobachtungszeit

n     <- sample_size                    # length of vector
rho   <- -0.5                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)     
[1] -0.5
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  336.4   339.6   340.1   340.0   340.5   341.9 
cor(data$number_of_falls, data$time_at_risk_days)
[1] -0.5
cor(data$number_of_falls, data$time_at_risk_days)
[1] -0.5
ggplot(data, aes(x=time_at_risk_days, y=number_of_falls))+
  geom_point()+
  geom_smooth(method="lm")
`geom_smooth()` using formula = 'y ~ x'

78.2.7 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 402
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.002955882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.078897

78.2.8 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.081207

Nun sehen wir, dass die Incidence Rate mit den individuellen Beobachtungszeiten höher sind, als die Berechnung mit der globalen Berechnung.

  • possitive Korrelation –> individuelle Berechnung tiefer als bei globaler Berechnung –> die höheren Beobachtungszeiten haben höhere absolute Sturzzahlen, da in der individuellen Formel die Sturzhäufigkeiten auf ein Jahr hochgerechnet werden, haben alle Personen das gleiche Gewicht. Hier möchten wir jedoch Personen mit längerer Beobachtungszeit mehr Gewicht geben.
  • negative Korrelation –> individuelle Berechnung höher als bei globaler Berechnung –> die höheren Beobachtungszeiten haben tiefere absolute Sturzzahlen, dies wird bei der globalen Berechnung berücksichtigt, bei der individuellen Berechnung habe alle Beobachtungen gleich viel Gewicht, d.h. die kürzeren Beobachtungszeiten bekommen gleich viel Gewicht wie die längeren Beobachtungszeiten, was relativ gesehen den kürzeren Zeiten (mit höheren Sturzhäufigkeiten) zu viel Gewicht gibt.

Welche Berechnungsform ist nun besser?

Das Problem liegt nicht wirklich in der individuellen Berechnung; die ist nämlich absolut in Ordnung.

Das Problem liegt darin, dass wir hier einen ungewichteten Mittelwert berechnet haben. Das ist natürlich nicht ok. Mit diesem ungewichteten Mittelwert bekommen alle Personen gleich viel Gewicht, was natürlich den Personen mit kürzeren Beobachtungsdauer gleich viel Gewicht gibt wie denen mit langer Beobachtungsdauer - was ja irgendwie nicht korrekt ist.

Wenn wir nun den gewichteten Mittelwert berechnen, so erhalten wir genau das gleiche Resultat wie mit der globalen Formel.

\[\begin{equation} \frac{\sum{(StürzeProZeit \cdot ZeitUnterBeobachtung)}}{\sum{ZeitUnterBeobachtung}} \tag{78.1} \end{equation}\]

Diese Formel ist äequivalent zu der Formel

\[\begin{equation} \frac{\sum{Stürze Aller Teilnehmenden}}{\sum{ZeitUnterBeobachtung}} \tag{78.2} \end{equation}\]

78.3 Berechnung mit einer negativen binomialen Regression.

Hierfür nehmen wir ein etwas realistischeres Beispiel und berechnen zuerst wieder mit beiden Ansätzen, danach mit einer negativ binomialen Regression und einer Poisson Regression.

n     <- sample_size                    # length of vector
rho   <- 0.3                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho
[1] 0.3
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  338.1   339.5   340.0   340.0   340.4   342.0 
cor(data$number_of_falls, data$time_at_risk_days)
[1] 0.3

78.3.1 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 402
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.002955882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.078897

78.3.2 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)
rio::export(data, "example_negative_binomial.dta")

mean(data$fall_per_year_per_person)
[1] 1.077653

78.3.3 Negative binomiale Regression

summary(negbinomial_model <- MASS::glm.nb(number_of_falls ~ 1+offset(log(time_at_risk_days/365)) , data = data))

Call:
MASS::glm.nb(formula = number_of_falls ~ 1 + offset(log(time_at_risk_days/365)), 
    data = data, init.theta = 0.4466212457, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0272  -1.0259  -1.0252  -0.0014   3.1703  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.07514    0.08993   0.836    0.403

(Dispersion parameter for Negative Binomial(0.4466) family taken to be 1)

    Null deviance: 340.27  on 399  degrees of freedom
Residual deviance: 340.27  on 399  degrees of freedom
AIC: 1084

Number of Fisher Scoring iterations: 1

              Theta:  0.4466 
          Std. Err.:  0.0608 

 2 x log-likelihood:  -1080.0340 
jtools::summ(negbinomial_model, exp=TRUE, confint=TRUE, digit=3)
Error in glm.control(...) : 
  unused arguments (offset = c(-0.0717326265478253, -0.0703539966177321, -0.0698133896732971, -0.0700385720495962, -0.0705653226017391, -0.0706654044286179, -0.069372895424496, -0.0701410615253942, -0.0726209145137364, -0.0696404729776312, -0.0711065510312744, -0.070175943230548, -0.0718867810282543, -0.0699419669830161, -0.0709586339620409, -0.0677520021384915, -0.0705149231821088, -0.0740780399555779, -0.0735002757815042, -0.0744832632898757, -0.0669833492063024, -0.0719117863169038, -0.0693819233477633, 
-0.0733025545320889, -0.0717655851621572, -0.072421609024525, -0.0717722870531041, -0.0703690248952421, -0.0689984452146951, -0.0677248367776067, -0.0725050173590278, -0.0730697713299173, -0.0720509763240432, -0.0703664142949842, -0.0730145512376663, -0.0703208307676104, -0.0706587317097451, -0.0718871120043004, -0.0749974183872817, -0.068621957606525, -0.0699927301682643, -0.0710180023355251, -0.072205508683922, -0.0652199790511748, -0.0712821573311864, -0.0727618108654913, -0.0713932073597814, 
-0.0701220515444283, -0.0737320186126656, -0.0671673644967413, -0.0694688123619422, -0.0766376486719216, -0.0681527899516639, -0.0701135649829617, -0.0711395697032553, -0.0723179612062762, -0.0737756603913992, -0.0727767930962053, -0.0716542867604467, -0.0729257141835503, -0.0725080897576667, -0.0711227937987278, -0.0732536428499788, -0.0699964966938734, -0.0679015996329161, -0.068339283133052, -0.0711887534273264, -0.0706750376837532, -0.0700350912887492, -0.0690137759972215, -0.0681994476909769, 
-0.0715979871437912, -0.0727677169018155, -0.0711429427723367, -0.069726606318111, -0.0685567037663104, -0.0715148948893871, -0.0695902240125851, -0.0707991176528878, -0.0697329792136626, -0.0705123062181056, -0.0719518450181477, -0.0722842568528479, -0.0732225888921247, -0.0680423289315448, -0.0707470738913164, -0.0706578287463643, -0.074579231924022, -0.0741003069199608, -0.0682721069117699, -0.0732315151909405, -0.0712034389698701, -0.0724765720970332, -0.0717072977840475, -0.0708953338251607, 
-0.0706207013729129, -0.0711027253424281, -0.0709531314207932, -0.0702376065866295, -0.0687317028613966, -0.0684269947396981, -0.0729264765220373, -0.0687146090771676, -0.0721479395400269, -0.0714495529047439, -0.0694920137488919, -0.0689594092917127, -0.0696103640480073, -0.0693334988290038, -0.0689536546650167, -0.0692923170736788, -0.070687908029419, -0.0705284815070228, -0.0723986489406819, -0.0719373860466306, -0.0703718556388417, -0.0702888682560468, -0.0728411836032829, -0.070908826000875, 
-0.0686158312569195, -0.0709283287185805, -0.0689608170130467, -0.0703210924533979, -0.0729994320555412, -0.0729172166577399, -0.0705100248320235, -0.0706834713826818, -0.0699777256874866, -0.0695399528016594, -0.0739406580869321, -0.0690515377018701, -0.0701262817529155, -0.0729907509778227, -0.0699462735137674, -0.0750658000227377, -0.0748185706400213, -0.0718089310578512, -0.0675640340418058, -0.0694957476580504, -0.0735017179017247, -0.0711936130290188, -0.0715563085585809, -0.0711969872372015, 
-0.068141565196159, -0.0717873400562467, -0.069551866478808, -0.0729630175871226, -0.0709583782064442, -0.0715712348979839, -0.0698422020430729, -0.0700939321504124, -0.0746217156741855, -0.0727877965103801, -0.0685659795903609, -0.0712456500576775, -0.0702266364436788, -0.0689543830253089, -0.0733920841605624, -0.0695441727853988, -0.0694282776280983, -0.0686299715372868, -0.0740045724424149, -0.0688184465640763, -0.0731417142611463, -0.0707263946810963, -0.0680494201751091, -0.0712502392416455, 
-0.0684112335146856, -0.0737629283562855, -0.071175605579386, -0.0691136231460267, -0.0707107809660952, -0.072237820504059, -0.0709728164517749, -0.0705206666974762, -0.0707695958332427, -0.0683057024990633, -0.0706971204671054, -0.0683189772093647, -0.0717285291044729, -0.0708846719640111, -0.0713813262258521, -0.071612024576405, -0.0698667601123456, -0.069972663171551, -0.0717321591842828, -0.072178379686446, -0.0718611172055841, -0.0728423535067146, -0.0736535210037957, -0.0673238990544608, -0.0723044387049933, 
-0.0725937737312853, -0.0723399085918712, -0.0736801647016733, -0.0702211240047069, -0.0729859891006385, -0.0719637190166357, -0.0734601404022343, -0.0709880698091787, -0.0732754839940353, -0.0711174840111001, -0.0715238741830342, -0.0716556311729074, -0.0702035405110388, -0.0704482999306964, -0.0702790947467147, -0.0686329151924903, -0.0686069162227465, -0.0713817895830631, -0.0679839692992444, -0.06802566739774, -0.07218794270634, -0.0725007320506683, -0.0692611576133595, -0.067508208809327, -0.0731407970515015, 
-0.0723277696624863, -0.0708401269977342, -0.0689876754706456, -0.0706936581191554, -0.0724350632737513, -0.0684465632684041, -0.0680283417806704, -0.0695167153971896, -0.0716100929426023, -0.0708888668386387, -0.0696013583981321, -0.0716896172678894, -0.0695236596904689, -0.0723766020571118, -0.0731168682437001, -0.0704649729945247, -0.0695836982488852, -0.0697489939839756, -0.0730145388802594, -0.0713451190125861, -0.0686689315023757, -0.0720452387879767, -0.0721679871385046, -0.07227879576735, 
-0.0683030675802553, -0.0717341360203211, -0.0693218619218664, -0.0698467126664045, -0.0709628966612516, -0.0706027192557368, -0.0703159962612394, -0.0708089362277375, -0.0695574323153742, -0.0690354861744911, -0.0700389917354469, -0.0686251827939006, -0.0708139061064606, -0.0675922736799527, -0.0726999840566022, -0.0707168310002694, -0.070913322660224, -0.0712276495176074, -0.075108910345449, -0.0720036924264521, -0.0721728584215209, -0.0714884771571677, -0.07250609674838, -0.0715331410934775, -0.0701004426785778, 
-0.0684020704838686, -0.0715634941209588, -0.071790190180531, -0.0729051113819281, -0.0694674998950155, -0.0725920838789434, -0.0705866818297941, -0.0684008029215127, -0.0665680094699465, -0.0733599102897859, -0.069686821165149, -0.0703195705940046, -0.0697091728350947, -0.0686667827941409, -0.0679880779792389, -0.0700356572575936, -0.0715239347228351, -0.0710852958750256, -0.0711854046562165, -0.0721480135322643, -0.0705312360175185, -0.0712651661737678, -0.0721411710970086, -0.0722446844514777, 
-0.0727063240624525, -0.0719773884793783, -0.0709633003452211, -0.0706074314436247, -0.0665693456657847, -0.067590367842326, -0.0714371816500099, -0.0734464103971011, -0.0681251416698172, -0.0717517548950393, -0.0709755781761124, -0.0688468649881978, -0.074017762323839, -0.0729237498130097, -0.0702042409764206, -0.072168775126775, -0.0724349245764692, -0.0694456599542266, -0.0700539984980535, -0.0719695458779297, -0.0706875049048342, -0.0708297148678662, -0.0741184430430564, -0.072837285984327, -0.0723805607212887, 
-0.0710637144240753, -0.0723783165383631, -0.0734841039989656, -0.0704933289693316, -0.0703565081428539, -0.0729035358648481, -0.0675100403118587, -0.0721131456867718, -0.0738087714527264, -0.0707899556438699, -0.0726345379397343, -0.0673176284140932, -0.0728155641057648, -0.0719793096414069, -0.0686332470185777, -0.0730643163485924, -0.071954934323424, -0.073670613294234, -0.0682018514750138, -0.0673909107978519, -0.0699582962011149, -0.0692805826728297, -0.073387936957518, -0.0732621569068001, 
-0.0726408594315515, -0.0743554148865078, -0.070381780483715, -0.0688676541078936, -0.0725119284218057, -0.0732588708497333, -0.0745775180820139, -0.0742995966046038, -0.0673062016367619, -0.0746150155101856, -0.0736021584506963, -0.0736528078962987, -0.0734186172131393, -0.0693213819845178, -0.070059999824277, -0.0732000865373183, -0.0729716391621944, -0.070526705959431, -0.0717785949086412, -0.0712441405596515, -0.0725039025324829, -0.0709771112812589, -0.0727275747045564, -0.0675241362544294, 
-0.0727120313077223, -0.0740844363690091, -0.0716687420938349, -0.0707360162986032, -0.0710289159302139, -0.0732837207979998, -0.0681063638599286, -0.0707524501562123, -0.071868327625517, -0.0698212651217147, -0.0694068760658432, -0.0711387487430101, -0.0713050390465972, -0.0705110065548295, -0.0722700186003683, -0.0702967927526921, -0.0697112427447397, -0.0726255916841179, -0.066161721071699, -0.0707546659011677, -0.0698692168977766, -0.0674835230719856, -0.0720686139841197, -0.0672977365825958, 
-0.069761901735935, -0.0725609993137609, -0.0700041165042845
Warning: Something went wrong when calculating the pseudo R-squared. Returning NA
instead.
Observations 400
Dependent variable number_of_falls
Type Generalized linear model
Family Negative Binomial(0.4466)
Link log
χ²(NA) NA
Pseudo-R² (Cragg-Uhler) NA
Pseudo-R² (McFadden) NA
AIC 1084.034
BIC 1092.017
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 1.078 0.904 1.286 0.836 0.403
Standard errors: MLE

78.3.4 Mit Stata

use example_negative_binomial.dta
gen time_at_risk_years=time_at_risk_days/365
nbreg number_of_falls , exposure(time_at_risk_years)
matrix b  = e(b)
dis exp(b[1,1])
Fitting Poisson model:

Iteration 0:   log likelihood = -708.37783  
Iteration 1:   log likelihood = -708.37783  

Fitting constant-only model:

Iteration 0:   log likelihood = -555.66995  
Iteration 1:   log likelihood =  -542.4481  
Iteration 2:   log likelihood = -540.01838  
Iteration 3:   log likelihood = -540.01681  
Iteration 4:   log likelihood = -540.01681  

Fitting full model:

Iteration 0:   log likelihood = -540.01681  
Iteration 1:   log likelihood = -540.01681  

Negative binomial regression                            Number of obs =    400
                                                        LR chi2(0)    =   0.00
Dispersion: mean                                        Prob > chi2   =      .
Log likelihood = -540.01681                             Pseudo R2     = 0.0000

------------------------------------------------------------------------------
number_of_~s | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   .0751405   .0899175     0.84   0.403    -.1010945    .2513756
ln(time_~rs) |          1  (exposure)
-------------+----------------------------------------------------------------
    /lnalpha |   .8060444   .1360249                      .5394405    1.072648
-------------+----------------------------------------------------------------
       alpha |   2.239034   .3045643                      1.715047     2.92311
------------------------------------------------------------------------------
LR test of alpha=0: chibar2(01) = 336.72               Prob >= chibar2 = 0.000


1.0780356

78.4 Berechnung mit einer Poisson Regression

poisson_model<-glm(formula = number_of_falls ~ 1+offset(log(time_at_risk_days/365)), family = poisson, data = data)
summary(poisson_model)

Call:
glm(formula = number_of_falls ~ 1 + offset(log(time_at_risk_days/365)), 
    family = poisson, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4209  -1.4178  -1.4161  -0.0033   8.3516  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.07594    0.04988   1.523    0.128

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 994.38  on 399  degrees of freedom
Residual deviance: 994.38  on 399  degrees of freedom
AIC: 1418.8

Number of Fisher Scoring iterations: 6
jtools::summ(poisson_model, exp=TRUE, confint=TRUE, digit=5)
Observations 400
Dependent variable number_of_falls
Type Generalized linear model
Family poisson
Link log
χ²(0) -0.00000
Pseudo-R² (Cragg-Uhler) 0.00000
Pseudo-R² (McFadden) 0.00000
AIC 1418.75567
BIC 1422.74713
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 1.07890 0.97842 1.18969 1.52258 0.12786
Standard errors: MLE

78.4.1 Poisson oder negative binomial Regression?

Wir können folgenden test durchführen siehe hier für mehr Erklärungen.

Wir können die beiden Log-Likelihoods der beiden Modelle vergleichen (wir dürfen dies, da das Poisson Modell im negative binomial model nested ist)

logLik(negbinomial_model)
'log Lik.' -540.0168 (df=2)
logLik(poisson_model)
'log Lik.' -708.3778 (df=1)

Die Log-Likelihood des negative binomialen Modells ist grösser und somit besser, deshalb sollten wir hier das negative biniomale Modell wählen. Wir können dies auch noch statistisch testen, hier der P-Wert des Log-Likelihood Tests:

pchisq(2 * (logLik(negbinomial_model) - logLik(poisson_model)), df = 1, lower.tail = FALSE)
'log Lik.' 3.301876e-75 (df=2)

78.5 Das Ganze mit Simulation nach einer negativ binomialen Verteilung

set.seed(666)
sample_size=400
number_of_falls<-rnbinom(sample_size, mu=1.4, size=1) # We could also simulate data based on a negative binomial distribution
hist(number_of_falls)

Im ersten Beispiel nehmen wir an, alle Teilnehmern sind genau 365 Tage in der Studie gewesen.

time_at_risk_days=rep(365, times=sample_size)

id<-1:sample_size

data<-data.frame(id, number_of_falls, time_at_risk_days)
rm(time_at_risk_days, id, number_of_falls)

78.6 Berechnung der Inzidenz, Formel 1

\[\begin{equation} \frac{\sum{Stürze Aller Teilnehmenden}}{\sum{ZeitUnterBeobachtung}} \tag{78.1} \end{equation}\]

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 572
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 146000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.003917808
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) # Diese und die nächsten Formeln sind equivalent, einfach anders ausformuliert
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.43 
cat("Zweite Formulierung, gleiches Resultat: ", (incidence_rate_per_year<- (Total_Number_of_Falls / Total_Exposure_Time) *365))
Zweite Formulierung, gleiches Resultat:  1.43
cat("\nDritte Formulierung, gleiches Resultat: ", (incidence_rate_per_year<- Total_Number_of_Falls / Total_Exposure_Time *365))

Dritte Formulierung, gleiches Resultat:  1.43

78.6.1 Berechnung mit individueller Zeit

Wir können jetzt auch die Berechnungen mit den individuellen Zeiten machen - die jedoch in diesem Beispiel alle gleich sind.

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.43

Die beiden Berechnungsvarianten ergeben genau das gleiche Resultat.

78.7 Im nächsten Beispiel wurden nicht alle Teilnehmenden über 365 Tage verfolgt.

Wir simulieren die Daten zunächst so, dass es keine Korrelation zwischen Anzahl Stürze und Beobachtungszeit gibt (was ja eigentlich nicht ganz realistisch ist). Danke an diese Webseite: https://stats.stackexchange.com/questions/15011/generate-a-random-variable-with-a-defined-correlation-to-an-existing-variables.

n     <- sample_size                    # length of vector
rho   <- 0.0                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho
[1] 4.019221e-16
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  338.1   339.6   340.0   340.0   340.4   341.8 
cor(data$number_of_falls, data$time_at_risk_days)
[1] 3.228183e-15
cor(data$number_of_falls, data$time_at_risk_days)
[1] 3.228183e-15
ggplot(data, aes(x=time_at_risk_days, y=number_of_falls))+
  geom_point()+
  geom_smooth(method="lm")
`geom_smooth()` using formula = 'y ~ x'

78.7.1 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 572
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.004205882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.535147

78.7.2 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.535152

Die beiden Berechnungen stimmen praktisch überein - nicht 100%, da ja die Korrelation nur in der Population 0 ist und nicht in der Stichprobe.

78.7.3 Korrelation zwischen Anzahl Stürze und Beobachtungszeit.

n     <- sample_size                    # length of vector
rho   <- 0.5                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)     
[1] 0.5
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  337.8   339.5   340.0   340.0   340.5   341.8 
cor(data$number_of_falls, data$time_at_risk_days)
[1] 0.5

78.7.4 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 572
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.004205882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.535147

78.7.5 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.533195

Nun sehen wir, dass die Werte nicht mehr exakt übereinstimmen.

78.7.6 negative Korrelation zwischen Anzahl Stürzen und Beobachtungszeit

n     <- sample_size                    # length of vector
rho   <- -0.5                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)     
[1] -0.5
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  337.6   339.5   340.1   340.0   340.5   341.9 
cor(data$number_of_falls, data$time_at_risk_days)
[1] -0.5
cor(data$number_of_falls, data$time_at_risk_days)
[1] -0.5
ggplot(data, aes(x=time_at_risk_days, y=number_of_falls))+
  geom_point()+
  geom_smooth(method="lm")
`geom_smooth()` using formula = 'y ~ x'

78.7.7 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 572
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.004205882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.535147

78.7.8 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)

mean(data$fall_per_year_per_person)
[1] 1.537117

Nun sehen wir, dass die Incidence Rate mit den individuellen Beobachtungszeiten höher sind, als die Berechnung mit der globalen Berechnung.

  • possitive Korrelation –> individuelle Berechnung tiefer als bei globaler Berechnung –> die höheren Beobachtungszeiten haben höhere absolute Sturzzahlen, da in der individuellen Formel die Sturzhäufigkeiten auf ein Jahr hochgerechnet werden, haben alle Personen das gleiche Gewicht. Hier möchten wir jedoch Personen mit längerer Beobachtungszeit mehr Gewicht geben.
  • negative Korrelation –> individuelle Berechnung höher als bei globaler Berechnung –> die höheren Beobachtungszeiten haben tiefere absolute Sturzzahlen, dies wird bei der globalen Berechnung berücksichtigt, bei der individuellen Berechnung habe alle Beobachtungen gleich viel Gewicht, d.h. die kürzeren Beobachtungszeiten bekommen gleich viel Gewicht wie die längeren Beobachtungszeiten, was relativ gesehen den kürzeren Zeiten (mit höheren Sturzhäufigkeiten) zu viel Gewicht gibt.

Welche Berechnungsform ist nun besser?

Das Problem liegt nicht wirklich in der individuellen Berechnung; die ist nämlich absolut in Ordnung.

Das Problem liegt darin, dass wir hier einen ungewichteten Mittelwert berechnet haben. Das ist natürlich nicht ok. Mit diesem ungewichteten Mittelwert bekommen alle Personen gleich viel Gewicht, was natürlich den Personen mit kürzeren Beobachtungsdauer gleich viel Gewicht gibt wie denen mit langer Beobachtungsdauer - was ja irgendwie nicht korrekt ist.

Wenn wir nun den gewichteten Mittelwert berechnen, so erhalten wir genau das gleiche Resultat wie mit der globalen Formel.

\[\begin{equation} \frac{\sum{(StürzeProZeit \cdot ZeitUnterBeobachtung)}}{\sum{ZeitUnterBeobachtung}} \tag{78.1} \end{equation}\]

Diese Formel ist äequivalent zu der Formel

\[\begin{equation} \frac{\sum{Stürze Aller Teilnehmenden}}{\sum{ZeitUnterBeobachtung}} \tag{78.2} \end{equation}\]

78.8 Berechnung mit einer negativen binomialen Regression.

Hierfür nehmen wir ein etwas realistischeres Beispiel und berechnen zuerst wieder mit beiden Ansätzen, danach mit einer negativ binomialen Regression und einer Poisson Regression.

n     <- sample_size                    # length of vector
rho   <- 0.3                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- data$number_of_falls       # fixed given data
x2    <- rnorm(n, 340, 12)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho
[1] 0.3
data$time_at_risk_days<-x*12+340
summary(data$time_at_risk_days)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  338.1   339.6   340.0   340.0   340.5   341.7 
cor(data$number_of_falls, data$time_at_risk_days)
[1] 0.3

78.8.1 Berechnung der Inzidenzrate mit der gleichen Formel wie oben:

(Total_Number_of_Falls<-sum(data$number_of_falls) )
[1] 572
(Total_Exposure_Time<-sum(data$time_at_risk_days) )
[1] 136000
incidence_rate_per_day<-Total_Number_of_Falls / Total_Exposure_Time
cat("Inzidenz-Rate pro Tag:", incidence_rate_per_day)
Inzidenz-Rate pro Tag: 0.004205882
incidence_rate_per_year<- Total_Number_of_Falls / (Total_Exposure_Time / 365) 
cat("\nInzidence rate per year:", incidence_rate_per_year, "\n")

Inzidence rate per year: 1.535147

78.8.2 Berechnung mit individueller Zeit

data<-data %>% 
  mutate(fall_per_year_per_person=number_of_falls/time_at_risk_days*365)
rio::export(data, "example_negative_binomial.dta")
mean(data$fall_per_year_per_person)
[1] 1.534086

78.9 Berechnung mit einer negativen binomialen Regression.

summary(negbinomial_model <- MASS::glm.nb(number_of_falls ~ 1+offset(log(time_at_risk_days/365)) , data = data))

Call:
MASS::glm.nb(formula = number_of_falls ~ 1 + offset(log(time_at_risk_days/365)), 
    data = data, init.theta = 1.218133634, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3773  -1.3748  -0.2658   0.2964   2.5452  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.42825    0.06165   6.946 3.76e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(1.2181) family taken to be 1)

    Null deviance: 420.72  on 399  degrees of freedom
Residual deviance: 420.72  on 399  degrees of freedom
AIC: 1318.7

Number of Fisher Scoring iterations: 1

              Theta:  1.218 
          Std. Err.:  0.183 

 2 x log-likelihood:  -1314.730 
jtools::summ(negbinomial_model, exp=TRUE, confint=TRUE, digit=3)
Error in glm.control(...) : 
  unused arguments (offset = c(-0.068968676944569, -0.0697991498529847, -0.0737634530365079, -0.0736600075581091, -0.0712781836646401, -0.0724104524645587, -0.0694740134529092, -0.068420938566686, -0.0671135164379177, -0.0711039134371102, -0.0709155339980465, -0.0716035942956918, -0.0699977579649627, -0.0733264142656758, -0.0718563890406658, -0.0707415059241687, -0.0699029929028674, -0.0695041798678721, -0.0707687798143977, -0.0714717394990612, -0.0709209049209726, -0.0709449242688688, -0.0690501797899614, 
-0.0702037198547886, -0.0692898465978144, -0.0667820099321513, -0.0689574158686012, -0.0732628682990992, -0.0677402908098551, -0.0699490067757682, -0.0717060155350405, -0.0713961814823604, -0.0681120068077467, -0.0687413724114093, -0.0704348689886693, -0.071385912529988, -0.0680272041386719, -0.0707521226021123, -0.0702998440082665, -0.0745829190273382, -0.0747919116845183, -0.0721920616506488, -0.074839689078549, -0.0710444105202445, -0.0736485183536446, -0.0687755739619554, -0.0704054261682056, 
-0.0737130729291081, -0.068111982048286, -0.0716688133217535, -0.0731067094313954, -0.0708576372213558, -0.0707277662145798, -0.0692235380518121, -0.071321395522611, -0.0704009288464188, -0.0703314005800436, -0.0719057351800459, -0.069079689350198, -0.071680456797708, -0.06862743476505, -0.0699815167499935, -0.0733990991325681, -0.0691157959010355, -0.0688770583872399, -0.0717666499354607, -0.0712654529557867, -0.070261656360919, -0.0676002150080896, -0.0714991483881176, -0.0691686543923234, -0.0727745825991636, 
-0.0704336409841636, -0.0699738060248732, -0.0668620204747172, -0.0708749004423805, -0.0710337440643337, -0.0708981681868791, -0.0703099977157452, -0.0717895103915943, -0.0709003175863271, -0.0714085833176496, -0.0709910726834786, -0.0713078137006905, -0.0713794243999824, -0.0723550130688308, -0.0723252419936409, -0.0717661379896117, -0.0713389018238879, -0.0692968816414703, -0.0703973457911674, -0.0716710165857714, -0.0723998039517327, -0.0694105049061337, -0.0688116842687969, -0.0713044343382779, 
-0.0722102342742904, -0.07115143553691, -0.0705477866514349, -0.0695188260971974, -0.0730290505715778, -0.0695051137147795, -0.0716131139168621, -0.0714690823986951, -0.0705983681854955, -0.0752615502263946, -0.0719446129411068, -0.0682403406805914, -0.0713255282514406, -0.0706172448503127, -0.0697419113218802, -0.0716457178544998, -0.0716506966355823, -0.068991491368312, -0.0726166463123554, -0.0689164524303162, -0.0723460959533967, -0.0719786520249702, -0.072042787064502, -0.072821689301542, -0.0705547203285491, 
-0.0677472389835865, -0.0718870083460832, -0.0711374813108231, -0.0733425227462382, -0.0719127394442125, -0.0734029411285251, -0.0722066519260133, -0.073179974745751, -0.0738612019804595, -0.0683651457370465, -0.069180630609965, -0.0716077603412429, -0.0731056026246956, -0.0694575317159066, -0.0703859456905282, -0.0727005420433626, -0.0682516014564817, -0.0733011932240735, -0.0733347655310397, -0.0690090326102153, -0.0734131910955479, -0.0697023209749197, -0.0747043343034983, -0.0711563526661717, 
-0.0705322821096657, -0.0743873139737635, -0.0691298998867998, -0.0710466683587174, -0.0707559718442274, -0.070811351885487, -0.0696043486829533, -0.0737705037542945, -0.0735651314036602, -0.0728516188120079, -0.0707717626338588, -0.0724902666338875, -0.0719247013432658, -0.0681311363639439, -0.0708907520759019, -0.0707717683924887, -0.071525900948052, -0.0724261415544107, -0.0705576818901771, -0.0708133786746341, -0.0718645567536771, -0.0703319325574823, -0.0736426597649454, -0.0698580768178586, 
-0.0711756990747274, -0.0692518074353296, -0.0720473755058481, -0.0700197890109865, -0.0738419511832062, -0.0667637429955855, -0.0710198348409803, -0.0701062830499874, -0.0682546394004673, -0.0686056325737789, -0.0722379683583665, -0.0704874714709859, -0.0691110356427772, -0.0704856308650722, -0.0737587483831758, -0.0681638117469476, -0.0718732517542751, -0.0705588265637541, -0.0706275870054968, -0.0742771153452343, -0.0694457006958183, -0.0743730067910033, -0.0733795537031658, -0.069054859762211, 
-0.0713653126374744, -0.0713021963905755, -0.0708899091054867, -0.0697487854470766, -0.0706492697353122, -0.071862780147052, -0.0709350822941242, -0.070896486635807, -0.0692761762700546, -0.070592973354744, -0.0725948879851424, -0.0700720787138864, -0.0681442612457376, -0.0722594695337577, -0.069866016482192, -0.0694964161448693, -0.0709230427179947, -0.0764223804356864, -0.0727680981225187, -0.0679867562783139, -0.070344768124649, -0.0693087399355016, -0.067388704813847, -0.0717517659882705, -0.0684888593168522, 
-0.0729272622039071, -0.0723271231930488, -0.0677108583340867, -0.066975688921538, -0.0725729490052155, -0.0686030101938042, -0.0730757983950904, -0.07221204359675, -0.0722337118211256, -0.0706397770480092, -0.0676197788358755, -0.0725749590552142, -0.0740668683612802, -0.0696243106049556, -0.0690320103704146, -0.0708179853447513, -0.0728014060333791, -0.0729584988820738, -0.0725776004269808, -0.0716878935971817, -0.0712930573919346, -0.072049479200909, -0.0693258442333984, -0.0742101465464432, -0.0739679650818743, 
-0.0717066186716532, -0.0688151881093708, -0.06957286842177, -0.0686286335483264, -0.0746855882976034, -0.0698577476003184, -0.0722461231426789, -0.0684146023679939, -0.0688456749413525, -0.0701887828770924, -0.0718515308076314, -0.0710920020585773, -0.0710438722323253, -0.0737577756879568, -0.0728986155065347, -0.0678816792173357, -0.0720207332156938, -0.0723503161372634, -0.0706587951621789, -0.070334330657614, -0.0686405614281925, -0.0677852984034521, -0.0703874275254048, -0.0712511189626344, 
-0.0685190068928634, -0.0725624105122154, -0.0744612282683064, -0.0716885307149118, -0.0720255487847911, -0.068269402672371, -0.0697681819699314, -0.0716311187496996, -0.0722848962540236, -0.0690657170488812, -0.0720152497054469, -0.0711562638027921, -0.0698802492466996, -0.0724015511712743, -0.0702719744634681, -0.0717465752290297, -0.0695044418227888, -0.069321106576074, -0.0675020077286621, -0.0704533805885467, -0.0719692117073386, -0.071653190882551, -0.0692651451703671, -0.0679362618260636, 
-0.0753465240104577, -0.06917601082622, -0.0733060640460414, -0.0709198734709542, -0.0674785925522283, -0.0702310385556253, -0.0716266236053315, -0.0707840098952257, -0.0707882129660232, -0.0696890119700824, -0.0693532163021932, -0.0749324070873818, -0.0718106193312743, -0.0719306764450022, -0.0702988209387536, -0.0717527364255986, -0.0729842583349061, -0.0711006961724182, -0.0715601995565389, -0.0691752526867375, -0.0726498598141964, -0.0732210126030912, -0.0718523117284745, -0.071124154883815, 
-0.0685641953209176, -0.0690577193704192, -0.0695818284457453, -0.0667581368621873, -0.0700109171550977, -0.0729690541956926, -0.0733682636379316, -0.0680861122733448, -0.0717150083898187, -0.0723524887342725, -0.0713471476082131, -0.0712053495519193, -0.0720399299286299, -0.0710279471121786, -0.069780914671708, -0.0708611027077629, -0.0692527388526587, -0.0713909081910807, -0.0722701104059075, -0.0719161071372388, -0.0735993039984405, -0.073044715818802, -0.067745136803127, -0.0701551264749507, 
-0.0721240171031963, -0.0712978908668382, -0.0698589806002051, -0.0761438585352019, -0.070307714268863, -0.0659720819406725, -0.0692805152127871, -0.0713589540334185, -0.0702035941681447, -0.0738129472693291, -0.0741982919044186, -0.0682993217206927, -0.0686291061161914, -0.0728790868797975, -0.0701249781712195, -0.0724439887491804, -0.0731657766642315, -0.0717189007412009, -0.0692957505051303, -0.0687335314655417, -0.0735916258099715, -0.0735275081577112, -0.0676694961625887, -0.0695301039866577, 
-0.0689505673267044, -0.0683115438923095, -0.0708427509567494, -0.0720914751378338, -0.0734834197839452, -0.0709379747617577, -0.0706292863101127, -0.0728350275794555, -0.0686899161345738, -0.0710541091391243, -0.0715822740963603, -0.0699943456392014, -0.0745670401233685, -0.0711711800473823, -0.0718843827872101, -0.0721863433315149, -0.069748954446649, -0.0713945692907949, -0.0725811511455235, -0.0706214435231983, -0.0725463748457465, -0.0698154045981532, -0.070615106007038, -0.0704968044146147, 
-0.0706131853067918, -0.0702524283862846, -0.0730870755006845, -
Warning: Something went wrong when calculating the pseudo R-squared. Returning NA
instead.
Observations 400
Dependent variable number_of_falls
Type Generalized linear model
Family Negative Binomial(1.2181)
Link log
χ²(NA) NA
Pseudo-R² (Cragg-Uhler) NA
Pseudo-R² (McFadden) NA
AIC 1318.730
BIC 1326.713
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 1.535 1.360 1.732 6.946 0.000
Standard errors: MLE 
use example_negative_binomial.dta
gen time_at_risk_years=time_at_risk_days/365
nbreg number_of_falls , exposure(time_at_risk_years)
matrix b  = e(b)
dis exp(b[1,1])
Fitting Poisson model:

Iteration 0:   log likelihood = -726.74043  
Iteration 1:   log likelihood = -726.74043  

Fitting constant-only model:

Iteration 0:   log likelihood = -658.27734  
Iteration 1:   log likelihood = -657.36715  
Iteration 2:   log likelihood = -657.36522  
Iteration 3:   log likelihood = -657.36522  

Fitting full model:

Iteration 0:   log likelihood = -657.36522  
Iteration 1:   log likelihood = -657.36522  

Negative binomial regression                            Number of obs =    400
                                                        LR chi2(0)    =   0.00
Dispersion: mean                                        Prob > chi2   =      .
Log likelihood = -657.36522                             Pseudo R2     = 0.0000

------------------------------------------------------------------------------
number_of_~s | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   .4282527   .0616488     6.95   0.000     .3074232    .5490821
ln(time_~rs) |          1  (exposure)
-------------+----------------------------------------------------------------
    /lnalpha |  -.1973199   .1501024                     -.4915152    .0968754
-------------+----------------------------------------------------------------
       alpha |    .820928   .1232233                      .6116989    1.101723
------------------------------------------------------------------------------
LR test of alpha=0: chibar2(01) = 138.75               Prob >= chibar2 = 0.000


1.5345738

78.10 Berechnung mit einer Poisson Regression

poisson_model<-glm(formula = number_of_falls ~ 1+offset(log(time_at_risk_days/365)), family = poisson, data = data)
summary(poisson_model)

Call:
glm(formula = number_of_falls ~ 1 + offset(log(time_at_risk_days/365)), 
    family = poisson, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6947  -1.6903  -0.3805   0.4505   4.6571  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.42863    0.04181   10.25   <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: 818.65  on 399  degrees of freedom
Residual deviance: 818.65  on 399  degrees of freedom
AIC: 1455.5

Number of Fisher Scoring iterations: 5
jtools::summ(poisson_model, exp=TRUE, confint=TRUE, digit=5)
Observations 400
Dependent variable number_of_falls
Type Generalized linear model
Family poisson
Link log
χ²(0) -0.00000
Pseudo-R² (Cragg-Uhler) 0.00000
Pseudo-R² (McFadden) 0.00000
AIC 1455.48086
BIC 1459.47232
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 1.53515 1.41436 1.66625 10.25150 0.00000
Standard errors: MLE

78.10.1 Poisson oder negative binomial Regression?

Wir können folgenden test durchführen siehe hier für mehr Erklärungen.

Wir können die beiden Log-Likelihoods der beiden Modelle vergleichen (wir dürfen dies, da das Poisson Modell im negative binomial model nested ist)

logLik(negbinomial_model)
'log Lik.' -657.3652 (df=2)
logLik(poisson_model)
'log Lik.' -726.7404 (df=1)

Die Log-Likelihood des negative binomialen Modells ist grösser und somit besser, deshalb sollten wir hier das negative biniomale Modell wählen. Wir können dies auch noch statistisch testen, hier der P-Wert des Log-Likelihood Tests:

pchisq(2 * (logLik(negbinomial_model) - logLik(poisson_model)), df = 1, lower.tail = FALSE)
'log Lik.' 4.99432e-32 (df=2)