Statistik 6: Übung

Datensatz novanimal.csv

Hintergrund Daten

Im Forschungsprojekt NOVANIMAL wird u.a. der Frage nachgegangen, was es braucht, damit Menschen freiwillig weniger tierische Produkte konsumieren? Ein interessanter Ansatzpunkt ist die Ausser-Haus-Verpflegung. Gemäss der ersten in den Jahren 2014/2015 durchgeführten nationalen Ernährungserhebung menuCH essen 70 % der Bevölkerung zwischen 18 und 75 Jahren am Mittag auswärts (Bochud u. a. 2017). Daher rückt die Gastronomie als zentraler Akteur einer innovativen und nachhaltigen Ernährungswirtschaft ins Blickfeld. Welche Innovationen in der Gastronomie könnten dazu beitragen, den Pro-Kopf-Verbrauch an tierischen Nahrungsmitteln zu senken?

Dazu wurde u.a. ein Experiment in zwei Hochschulmensen durchgeführt. Forschungsleitend war die Frage, wie die Gäste dazu bewogen werden können, häufiger vegetarische oder vegane Gerichte zu wählen. Konkret wurde untersucht, wie die Gäste auf ein verändertes Menü-Angebot mit einem höheren Anteil an vegetarischen und veganen Gerichten reagieren. Das Experiment fand während 12 Wochen statt und bestand aus zwei Mensazyklen à 6 Wochen. Über den gesamten Untersuchungszeitraum werden insgesamt 90 verschiedene Gerichte angeboten. In den 6 Referenz- bzw. Basiswochen wurden zwei fleisch- oder fischhaltige Menüs und ein vegetarisches Menü angeboten. In den 6 Interventionswochen wurde das Verhältnis umgekehrt und es wurden ein veganes, ein vegetarisches und ein fleisch- oder fischhaltiges Gericht angeboten. Basis- und Interventionsangebote wechselten wöchentlich ab. Während der gesamten 12 Wochen konnten die Gäste jeweils auf ein Buffet ausweichen und ihre Mahlzeit aus warmen und kalten Komponenten selber zusammenstellen. Die Gerichte wurden über drei vorgegebene Menülinien (F, K, W) randomisiert angeboten.

Die Abbildung zeigt das Versuchsdesign der ersten 6 Experimentwochen (Kalenderwoche 40 bis 45).

Mehr Informationen über das Forschungsprojekt NOVANIMAL findet ihr auf der Webpage.

Aufgaben

Führt mit dem novanimal Datensatz (inviduelle Daten) eine logistische Regression durch, wobei ihr die einzelnen Käufer (single campus_card holder “ccrs”) als randomisierte Variable mitberücksichtigt. Kann der Fleischkonsum durch das Geschlecht, die Hochschulzugehörigkeit und das Alter erklärt werden? Berücksichtigt auch mögliche Interaktionen zwischen dem Geschlecht und dem Alter sowie dem Geshchlecht und der Hochschulzugehörigkeit

• Bestimmt das minimal adäquate Modell
• Stellt die Ergebnisse dar

Musterlösung Übung 6: GLMM

Demoscript herunterladen (.R)

Demoscript herunterladen (.qmd)

• Lösungstext als Download

# lade daten
nova <- read_delim("datasets/stat/novanimal.csv", delim = ";") |>
  mutate(across(where(is.character), as.factor))
str(nova)

tibble [17,900 × 6] (S3: tbl_df/tbl/data.frame)
 $ ccrs       : num [1:17900] 1e+09 1e+09 1e+09 1e+09 1e+09 ...
 $ meat       : num [1:17900] 1 0 1 0 0 0 0 0 0 0 ...
 $ gender     : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
 $ member     : Factor w/ 2 levels "Mitarbeitende",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ age_group  : Factor w/ 4 levels "AK1","AK2","AK3",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ age_group_2: Factor w/ 4 levels "15 - 25-jaehrig",..: 1 1 1 1 1 1 1 1 1 1 ...

# Die Variable“ccrs” (single campus_card holder) sollte ein facor sein
nova$ccrs <- as.factor(nova$ccrs)
summary(nova)

         ccrs            meat        gender              member      age_group 
 1000564796:   57   Min.   :0.0000   F: 5863   Mitarbeitende: 6868   AK1:7651  
 1000621354:   57   1st Qu.:0.0000   M:12037   Studierende  :11032   AK2:5720  
 1000594415:   54   Median :1.0000                                   AK3:3104  
 1000597794:   54   Mean   :0.6041                                   AK4:1425  
 1000608246:   54   3rd Qu.:1.0000                                             
 1000621447:   54   Max.   :1.0000                                             
 (Other)   :17570                                                              
          age_group_2  
 15 - 25-jaehrig:7651  
 26 - 34-jaehrig:5720  
 35 - 49-jaehrig:3104  
 50 - 64-jaehrig:1425  
                       
                       
                       

# sieht euch die Verteilung zwischen Fleisch und  kein Fleisch an,
# beide Kategorien kommen nicht gleich häufig vor, ist aber nicht tragisch
table(nova$meat) 

    0     1 
 7087 10813 

prop.table(table(nova$meat)) # Werte in Prozent

        0         1 
0.3959218 0.6040782 

# Definiert das logistische Modell mit ccrs als random intercept und
# wendet es auf den Datensatz an
glmm_1 <- glmmTMB(meat ~ gender + age_group + member + gender:age_group + gender:member + 
                    (1 | ccrs), 
                 family = binomial, 
                 data = nova)

Anova(glmm_1)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: meat
                    Chisq Df Pr(>Chisq)    
gender           125.9903  1     <2e-16 ***
age_group          4.3512  3     0.2260    
member             1.7696  1     0.1834    
gender:age_group   0.8800  3     0.8302    
gender:member      1.1107  1     0.2919    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(glmm_1)

 Family: binomial  ( logit )
Formula:          
meat ~ gender + age_group + member + gender:age_group + gender:member +  
    (1 | ccrs)
Data: nova

     AIC      BIC   logLik deviance df.resid 
 21128.1  21213.8 -10553.0  21106.1    17889 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 ccrs   (Intercept) 1.478    1.216   
Number of obs: 17900, groups:  ccrs, 1427

Conditional model:
                           Estimate Std. Error z value Pr(>|z|)   
(Intercept)               -0.060517   0.221841  -0.273  0.78501   
genderM                    0.825594   0.281970   2.928  0.00341 **
age_groupAK2              -0.076795   0.171426  -0.448  0.65417   
age_groupAK3              -0.023239   0.241615  -0.096  0.92338   
age_groupAK4               0.209746   0.313142   0.670  0.50298   
memberStudierende         -0.336658   0.203570  -1.654  0.09817 . 
genderM:age_groupAK2      -0.152603   0.213397  -0.715  0.47454   
genderM:age_groupAK3       0.002457   0.319039   0.008  0.99385   
genderM:age_groupAK4      -0.211615   0.410679  -0.515  0.60636   
genderM:memberStudierende  0.274622   0.260581   1.054  0.29194   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## erste Interpretation: Geschlecht (Referenzkategorie: Mann) scheint den Fleischkonsum positiv zu beeinflussen, schauen wir ob wir das Model vereinfachen können:

# Model optimierung
drop1(glmm_1, test = "Chi")

Single term deletions

Model:
meat ~ gender + age_group + member + gender:age_group + gender:member + 
    (1 | ccrs)
                 Df   AIC     LRT Pr(>Chi)
<none>              21128                 
gender:age_group  3 21123 0.87972   0.8303
gender:member     1 21127 1.11067   0.2919

glmm_2 <- update(glmm_1,~. -gender:age_group)

drop1(glmm_2, test = "Chi")

Single term deletions

Model:
meat ~ gender + age_group + member + (1 | ccrs) + gender:member
              Df   AIC    LRT Pr(>Chi)  
<none>           21123                  
age_group      3 21121 4.3315  0.22782  
gender:member  1 21124 3.2749  0.07035 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glmm_3 <- update(glmm_2,~. -gender:member)

drop1(glmm_3, test = "Chi")

Single term deletions

Model:
meat ~ gender + age_group + member + (1 | ccrs)
          Df   AIC     LRT Pr(>Chi)    
<none>       21124                     
gender     1 21245 123.112   <2e-16 ***
age_group  3 21123   4.256   0.2351    
member     1 21124   1.968   0.1607    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glmm_4 <- update(glmm_3,~. -age_group)

drop1(glmm_4, test = "Chi")

Single term deletions

Model:
meat ~ gender + member + (1 | ccrs)
       Df   AIC     LRT Pr(>Chi)    
<none>    21123                     
gender  1 21242 121.295   <2e-16 ***
member  1 21126   4.963   0.0259 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Model validierung
# Test overdisepersion
check_overdispersion(glmm_4)

# Overdispersion test

 dispersion ratio = 0.976
          p-value = 0.024

# der Wert ist nahe 1 daher i.o

set.seed(123)
simulationOutput <- simulateResiduals(fittedModel = glmm_4, plot = TRUE, n = 500)

# Plot residuals vs covariates des models
plotResiduals(simulationOutput, form = nova$gender)

plotResiduals(simulationOutput, form = nova$member)

# Die formalen Tests zeigen, dass es Probleme gibt. Die visuelle Inspektion zeigt jedoch, dass es überhaupt keine Probleme gibt und die Residiuen beinahe perfekt normalverteilt sind (dies zeigt, dass die Testergebnisse stark von der Anzahl der Beobachtungen (in diesem Fall viele) abhängen).

# Modelresultat
Anova(glmm_4)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: meat
          Chisq Df Pr(>Chisq)    
gender 123.6924  1    < 2e-16 ***
member   4.9512  1    0.02607 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(glmm_4)

 Family: binomial  ( logit )
Formula:          meat ~ gender + member + (1 | ccrs)
Data: nova

     AIC      BIC   logLik deviance df.resid 
 21122.5  21153.7 -10557.3  21114.5    17896 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 ccrs   (Intercept) 1.499    1.224   
Number of obs: 17900, groups:  ccrs, 1427

Conditional model:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -0.15918    0.08504  -1.872   0.0612 .  
genderM            0.93860    0.08439  11.122   <2e-16 ***
memberStudierende -0.19512    0.08769  -2.225   0.0261 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Pseudo R^2
r2(glmm_1)

# R2 for Mixed Models

  Conditional R2: 0.338
     Marginal R2: 0.040

# das marginale R^2 (r2m) gibt uns die erklärte Varianz der fixen Effekte: hier 4% (das ist sehr wenig)
# das conditionale R^2 (r2c) gibt uns die erklärte Varianz für das ganze Modell
# (mit fixen und variablen Effekten): hier 34% 

# Odds Ratios

# Extrahiere Estimates
estimates <- fixef(glmm_4)

# Berechne Odds Ratios
odsr <- rbind( 
  "Mitarbeiterinnen" = exp( estimates$cond["(Intercept)"] ), 
  "Mitarbeiter" = exp( estimates$cond["(Intercept)"] + estimates$cond["genderM"]),
  "Studentinnen" =  exp( estimates$cond["(Intercept)"] + estimates$cond["memberStudierende"]),
  "Studenten" = exp( estimates$cond["(Intercept)"] + estimates$cond["genderM"] + estimates$cond["memberStudierende"])
                 )
colnames(odsr) <- c("Odds_ratios")
odsr <- as.data.frame(odsr)

# Wahrscheinlichkeiten in %
odsr$Wahrscheinlichkeit <- ( odsr$Odds_ratios / (odsr$Odds_ratios + 1 ) )
odsr

                 Odds_ratios Wahrscheinlichkeit
Mitarbeiterinnen    0.852846          0.4602897
Mitarbeiter         2.180207          0.6855550
Studentinnen        0.701665          0.4123403
Studenten           1.793729          0.6420555

# Visualisierung
p_load("sjPlot")
plot_model(glmm_4, 
           type = "pred", pred.type = "re", 
           terms = c("member", "gender")) +
  theme_classic()

Bochud, Murielle, Angéline Chatelan, Juan-Manuel Blanco, und Sigrid Maria Beer-Borst. 2017. „Anthropometric characteristics and indicators of eating and physical activity behaviors in the Swiss adult population: results from menuCH 2014-2015“.