Kapitel 79 Epidemiologische Analysen

Dieses Kapitel ist noch nicht fertig und ist unkorrigiert.

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

library(epiR)

## Loading required package: survival

## Package epiR 2.0.53 is loaded

## Type help(epi.about) for summary information

## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses

##

library(sjPlot)

Lesen Sie bitte diesen Text durch: https://cran.r-project.org/web/packages/epiR/vignettes/epiR_measures_of_association.html.

data<-rio::import("titanic.RData")
head(data)

  reisende_id geschlecht ueberlebt langehaare died died.factor men men.factor
1           1          0         0          1    1        Died   1        Men
2           2          0         0          1    1        Died   1        Men
3           3          0         0          1    1        Died   1        Men
4           4          0         0          1    1        Died   1        Men
5           5          0         0          1    1        Died   1        Men
6           6          0         0          1    1        Died   1        Men
  short_hair short_hair.factor
1          0         Long Hair
2          0         Long Hair
3          0         Long Hair
4          0         Long Hair
5          0         Long Hair
6          0         Long Hair

tabulardata<-table(data[,c( "men.factor","died.factor")])
tabulardata

          died.factor
men.factor Survived Died
     Women      339  128
     Men        163  680

Jetzt müssen wir noch bestimmen, was als Outcome + bezeichnet wird. Schauen wir uns die Levels an:

levels(data$died.factor)

[1] "Survived" "Died"

levels(data$men.factor)

[1] "Women" "Men"

Bei der Variable died.factor ist Survided als erstes Level definiert, bei der Variable men.factor ist Women als erstes Level definiert. Das ist genau so, wie wir es nicht möchten.

Outcome+ und Exposed+ wird vom ersten Level definiert, deswegen wollen wir die Levels tauschen.

data$died.factor<-factor(data$died.factor, levels=c("Died", "Survived"))

tabulardata<-table(data[,c( "men.factor","died.factor")])
tabulardata

          died.factor
men.factor Died Survived
     Women  128      339
     Men    680      163

levels(data$died.factor)

[1] "Died"     "Survived"

data$men.factor<-factor(data$men.factor, levels=c("Men", "Women"))

tabulardata<-table(data[,c( "men.factor","died.factor")])
tabulardata

          died.factor
men.factor Died Survived
     Men    680      163
     Women  128      339

levels(data$men.factor)

[1] "Men"   "Women"

Nun können wie die Odds zu sterben für die Männer und die Odds zu sterben für die Frauen berechnen. Damit berechnen wir auch die Odds Ratio. Da die Daten aus einer Kohorte stammen und wir die Anzahl der Todesfälle haben, geben wir als Methode cohort.count an.

epi.2by2(dat = tabulardata, method = "cohort.count", conf.level = 0.95, units = 100, 
   interpret = FALSE, outcome = "as.columns")

             Outcome +    Outcome -      Total        Inc risk *        Odds
Exposed +          680          163        843              80.7       4.172
Exposed -          128          339        467              27.4       0.378
Total              808          502       1310              61.7       1.610

Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio                                 2.94 (2.53, 3.42)
Odds ratio                                     11.05 (8.47, 14.41)
Attrib risk in the exposed *                   53.26 (48.41, 58.10)
Attrib fraction in the exposed (%)            66.02 (60.47, 70.79)
Attrib risk in the population *                34.27 (29.44, 39.10)
Attrib fraction in the population (%)         55.56 (49.43, 60.95)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 360.601 Pr>chi2 = <0.001
Fisher exact test that OR = 1: Pr>chi2 = <0.001
 Wald confidence limits
 CI: confidence interval
 * Outcomes per 100 population units

Wir können die Odds Ratio auch mit einer logistischen Regression berechnen.

glm.or<-glm(died.factor~men.factor, data=data, family="binomial")
tab_model(glm.or)

	died factor
Predictors	Odds Ratios	CI	p
(Intercept)	0.24	0.20 – 0.28	<0.001
men factor [Women]	11.05	8.50 – 14.46	<0.001
Observations	1310
R² Tjur	0.275

Wir sind aber hier erstaunt, dass die Odds Ratio für die Frauen im Vergleich zu den Männern gilt, eigentlich sollte es umgekehrt sein. Hier haben wir das Problem wegen den Faktoren. Am besten ändern wir die abhängige Variable wieder in eine numerische Variable um.

data$died.num<-as.numeric(data$died.factor)
data$died.num<-recode(data$died.num, '2'=0)
table(data$died.num, data$died.factor)

   
    Died Survived
  0    0      502
  1  808        0

Died ist jetzt als 1 kodiert.

glm.or<-glm(died.num~men.factor, data=data, family="binomial")
tab_model(glm.or)

	died num
Predictors	Odds Ratios	CI	p
(Intercept)	4.17	3.53 – 4.96	<0.001
men factor [Women]	0.09	0.07 – 0.12	<0.001
Observations	1310
R² Tjur	0.275

Nun müssen wir wieder das Geschlecht anders kodieren.

data$men.factor<-factor(data$men.factor, levels=c("Women", "Men"))

Voilà, jetzt sind wir zufrieden.

glm.or<-glm(died.num~men.factor, data=data, family="binomial")
tab_model(glm.or)

	died num
Predictors	Odds Ratios	CI	p
(Intercept)	0.38	0.31 – 0.46	<0.001
men factor [Men]	11.05	8.50 – 14.46	<0.001
Observations	1310
R² Tjur	0.275

Wir machen den letzten Schritt wieder rückgängig und zeigen eine neue Variante.

data$men.factor<-factor(data$men.factor, levels=c("Men", "Women"))

glm.or<-glm(died.num~relevel(men.factor, ref="Women"), data=data, family="binomial")
tab_model(glm.or)

	died num
Predictors	Odds Ratios	CI	p
(Intercept)	0.38	0.31 – 0.46	<0.001
relevel(men factor, ref = “Women”) [Men]	11.05	8.50 – 14.46	<0.001
Observations	1310
R² Tjur	0.275

Die Risk Ratio (die im Befehl epi.2by2 als Incidence Risk Ratio angegeben wird), können wir auch mit einem generalisierten linearen Modell berechnen. Dieser Befehl akzeptiert jedoch die Faktorvariable died.factor nicht.

glm.rr<-glm(died.factor~men.factor, data=data, family="poisson")

Warning in Ops.factor(y, 0): '<' not meaningful for factors

Error in if (any(y < 0)) stop("negative values not allowed for the 'Poisson' family"): missing value where TRUE/FALSE needed

tab_model(glm.rr)

Error in tab_model(glm.rr): object 'glm.rr' not found

Wir müssen die Faktor Variable in eine numerische Variable umtauschen, und die 2 muss zur 0 werden.

data$died.num<-as.numeric(data$died.factor)
data$died.num<-recode(data$died.num, '2'=0)

glm.rr<-glm(died.num~men.factor, data=data, family="poisson")
tab_model(glm.rr)

	died num
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	0.81	0.75 – 0.87	<0.001
men factor [Women]	0.34	0.28 – 0.41	<0.001
Observations	1310
R² Nagelkerke	0.252