Multilevel Report

Multilevel Analysis- Report

#library----
library(dplyr)

Attache Paket: 'dplyr'

Die folgenden Objekte sind maskiert von 'package:stats':

    filter, lag

Die folgenden Objekte sind maskiert von 'package:base':

    intersect, setdiff, setequal, union

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ ggplot2   3.5.1     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(forcats)
library(haven)
library(nlme) # preferred package for (basic) multilevel models 

Attache Paket: 'nlme'

Das folgende Objekt ist maskiert 'package:dplyr':

    collapse

library(lme4)

Lade nötiges Paket: Matrix

Attache Paket: 'Matrix'

Die folgenden Objekte sind maskiert von 'package:tidyr':

    expand, pack, unpack


Attache Paket: 'lme4'

Das folgende Objekt ist maskiert 'package:nlme':

    lmList

library(stats)
library(ggplot2)
library(plotly)

Warning: Paket 'plotly' wurde unter R Version 4.4.3 erstellt

Attache Paket: 'plotly'

Das folgende Objekt ist maskiert 'package:ggplot2':

    last_plot

Das folgende Objekt ist maskiert 'package:stats':

    filter

Das folgende Objekt ist maskiert 'package:graphics':

    layout

library(psych)

Attache Paket: 'psych'

Die folgenden Objekte sind maskiert von 'package:ggplot2':

    %+%, alpha

library(moments)#skew
library(stargazer)#visualization table

Please cite as: 

 Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
 R package version 5.2.3. https://CRAN.R-project.org/package=stargazer 

#wd
setwd("C:/Users/valen/OneDrive/Documents/traineeship_PD/PHASE")
phase = read_sav("Dataset originale.sav")

Data preparation

Creating the main variables:

#create diet----
names(phase)[names(phase) == "Q1"] <- "diet"
names(phase)[names(phase) == "Q105"] <- "edo"
table(phase$diet)

   1    2    3    4    5 
3458  389   73  116   30 

names(phase)[names(phase) == "Q5"] <- "meat_frequency"
names(phase)[names(phase) == "Q161"] <- "pol_orientation"
names(phase)[names(phase) == "Q170"] <- "age"
names(phase)[names(phase) == "Q172"] <- "education"

phase$gender <- factor(phase$Q169,
                       levels = c(1, 2, 3, 4),
                       labels = c("Women", "Men", "Non-binary", "Prefer not to say"))

phase$meat_identity = rowMeans(phase[,c("Q12_Q14r1", "Q12_Q14r2", "Q12_Q14r3")])
#reversing Q91_Q104r4, Q91_Q104r14
phase$Q91_Q104r4 = 8 - phase$Q91_Q104r4 #reverse item 4
phase$Q91_Q104r14 = 8 - phase$Q91_Q104r14 #reverse item 14
phase$connect_nature = rowMeans(phase[,c("Q91_Q104r1", "Q91_Q104r2", "Q91_Q104r3", "Q91_Q104r4", 
                                         "Q91_Q104r5", "Q91_Q104r6", "Q91_Q104r7", "Q91_Q104r8",
                                         "Q91_Q104r9", "Q91_Q104r10", "Q91_Q104r11", "Q91_Q104r12", 
                                         "Q91_Q104r13", "Q91_Q104r14")])

Cronbach’s alpha for connect nature and meat eater identity

scale_meatid = phase[, c("Q12_Q14r1", "Q12_Q14r2", "Q12_Q14r3")]
scale_connectnat = phase[, c("Q91_Q104r1", "Q91_Q104r2", "Q91_Q104r3", "Q91_Q104r4", 
                            "Q91_Q104r5", "Q91_Q104r6", "Q91_Q104r7", "Q91_Q104r8",
                            "Q91_Q104r9", "Q91_Q104r10", "Q91_Q104r11", "Q91_Q104r12", 
                            "Q91_Q104r13", "Q91_Q104r14")]
alpha(scale_meatid)  

Reliability analysis   
Call: alpha(x = scale_meatid)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean  sd median_r
      0.83      0.83    0.77      0.63   5 0.0045  3.8 1.6     0.63

    95% confidence boundaries 
         lower alpha upper
Feldt     0.82  0.83  0.84
Duhachek  0.82  0.83  0.84

 Reliability if an item is dropped:
          raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
Q12_Q14r1      0.74      0.75    0.60      0.60 3.0   0.0080    NA  0.60
Q12_Q14r2      0.79      0.79    0.66      0.66 3.8   0.0066    NA  0.66
Q12_Q14r3      0.77      0.77    0.63      0.63 3.4   0.0071    NA  0.63

 Item statistics 
             n raw.r std.r r.cor r.drop mean  sd
Q12_Q14r1 3825  0.88  0.88  0.79   0.72  4.6 1.8
Q12_Q14r2 3825  0.85  0.86  0.74   0.67  2.9 1.8
Q12_Q14r3 3825  0.88  0.87  0.76   0.69  3.8 2.0

Non missing response frequency for each item
             1    2    3    4    5    6    7 miss
Q12_Q14r1 0.07 0.08 0.12 0.18 0.18 0.16 0.21 0.06
Q12_Q14r2 0.31 0.17 0.15 0.16 0.12 0.05 0.04 0.06
Q12_Q14r3 0.19 0.13 0.13 0.15 0.15 0.11 0.14 0.06

alpha(scale_connectnat) 

Reliability analysis   
Call: alpha(x = scale_connectnat)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd median_r
      0.89       0.9    0.92      0.38 8.6 0.0025  4.8 0.86     0.46

    95% confidence boundaries 
         lower alpha upper
Feldt     0.88  0.89  0.89
Duhachek  0.88  0.89  0.89

 Reliability if an item is dropped:
            raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
Q91_Q104r1       0.87      0.88    0.90      0.37  7.5   0.0029 0.062  0.45
Q91_Q104r2       0.87      0.88    0.90      0.36  7.4   0.0029 0.061  0.45
Q91_Q104r3       0.87      0.88    0.91      0.37  7.6   0.0029 0.063  0.45
Q91_Q104r4       0.89      0.90    0.92      0.41  9.1   0.0024 0.060  0.54
Q91_Q104r5       0.87      0.88    0.91      0.37  7.6   0.0029 0.063  0.45
Q91_Q104r6       0.87      0.88    0.91      0.37  7.6   0.0029 0.063  0.45
Q91_Q104r7       0.87      0.88    0.90      0.36  7.4   0.0030 0.061  0.45
Q91_Q104r8       0.87      0.88    0.91      0.37  7.6   0.0029 0.064  0.45
Q91_Q104r9       0.87      0.88    0.91      0.37  7.6   0.0029 0.063  0.45
Q91_Q104r10      0.87      0.88    0.91      0.37  7.5   0.0029 0.062  0.45
Q91_Q104r11      0.87      0.88    0.90      0.36  7.4   0.0030 0.060  0.45
Q91_Q104r12      0.91      0.91    0.92      0.44 10.3   0.0022 0.040  0.54
Q91_Q104r13      0.88      0.89    0.91      0.39  8.3   0.0027 0.065  0.54
Q91_Q104r14      0.90      0.91    0.92      0.43  9.9   0.0022 0.046  0.54

 Item statistics 
               n raw.r std.r r.cor r.drop mean  sd
Q91_Q104r1  4066  0.78  0.78 0.773  0.728  4.9 1.3
Q91_Q104r2  4066  0.80  0.81 0.808  0.761  5.1 1.3
Q91_Q104r3  4066  0.75  0.76 0.739  0.701  5.4 1.2
Q91_Q104r4  4066  0.40  0.39 0.324  0.293  4.6 1.4
Q91_Q104r5  4066  0.76  0.76 0.743  0.709  5.0 1.3
Q91_Q104r6  4066  0.77  0.77 0.747  0.710  4.9 1.4
Q91_Q104r7  4066  0.82  0.82 0.812  0.773  5.0 1.3
Q91_Q104r8  4066  0.75  0.76 0.738  0.704  5.1 1.2
Q91_Q104r9  4066  0.75  0.75 0.736  0.696  4.8 1.3
Q91_Q104r10 4066  0.77  0.78 0.763  0.723  4.9 1.3
Q91_Q104r11 4066  0.82  0.82 0.822  0.779  4.9 1.3
Q91_Q104r12 4066  0.15  0.14 0.061  0.017  3.8 1.6
Q91_Q104r13 4066  0.59  0.59 0.543  0.507  4.8 1.4
Q91_Q104r14 4066  0.22  0.21 0.136  0.094  4.2 1.6

Non missing response frequency for each item
               1    2    3    4    5    6    7 miss
Q91_Q104r1  0.02 0.02 0.06 0.28 0.31 0.19 0.12    0
Q91_Q104r2  0.02 0.02 0.05 0.20 0.35 0.21 0.15    0
Q91_Q104r3  0.01 0.01 0.03 0.15 0.33 0.27 0.19    0
Q91_Q104r4  0.02 0.05 0.14 0.27 0.25 0.15 0.13    0
Q91_Q104r5  0.02 0.02 0.06 0.27 0.29 0.20 0.15    0
Q91_Q104r6  0.02 0.03 0.08 0.27 0.27 0.18 0.15    0
Q91_Q104r7  0.02 0.02 0.06 0.27 0.30 0.19 0.15    0
Q91_Q104r8  0.01 0.01 0.05 0.21 0.36 0.22 0.14    0
Q91_Q104r9  0.02 0.02 0.07 0.30 0.32 0.18 0.10    0
Q91_Q104r10 0.02 0.02 0.07 0.27 0.31 0.18 0.12    0
Q91_Q104r11 0.02 0.02 0.07 0.27 0.32 0.19 0.11    0
Q91_Q104r12 0.11 0.09 0.16 0.32 0.20 0.08 0.04    0
Q91_Q104r13 0.03 0.03 0.09 0.26 0.29 0.16 0.13    0
Q91_Q104r14 0.05 0.09 0.19 0.28 0.19 0.10 0.10    0

Create Province groups to account for the provinces with too less participants

# Count of observations per province
table(haven::as_factor(phase$h_Province))

                      Torino                     Vercelli 
                         164                           14 
                      Novara                        Cuneo 
                          17                           31 
                        Asti                  Alessandria 
                          21                           25 
                      Biella         Verbano-Cusio-Ossola 
                           9                           10 
Valle d'Aosta/Vallée d'Aoste                       Varese 
                          28                           57 
                        Como                      Sondrio 
                          25                            6 
                      Milano                      Bergamo 
                         292                           63 
                     Brescia                        Pavia 
                          72                           29 
                     Cremona                      Mantova 
                          13                           17 
                       Lecco                         Lodi 
                          14                           20 
       Monza e della Brianza                Bolzano/Bozen 
                          65                           20 
                      Trento                       Verona 
                          50                           63 
                     Vicenza                      Belluno 
                          40                           11 
                     Treviso                      Venezia 
                          49                           63 
                      Padova                       Rovigo 
                          84                           19 
                       Udine                      Gorizia 
                          26                           16 
                     Trieste                    Pordenone 
                          25                           17 
                     Imperia                       Savona 
                          15                           14 
                      Genova                    La Spezia 
                          66                           10 
                    Piacenza                        Parma 
                          21                           34 
          Reggio nell'Emilia                       Modena 
                          18                           57 
                     Bologna                      Ferrara 
                          74                           20 
                     Ravenna                 Forlì-Cesena 
                          31                           22 
                      Rimini                Massa-Carrara 
                          24                           13 
                       Lucca                      Pistoia 
                          31                           15 
                     Firenze                      Livorno 
                          72                           25 
                        Pisa                       Arezzo 
                          25                           15 
                       Siena                     Grosseto 
                          22                           19 
                       Prato                      Perugia 
                          13                           45 
                       Terni              Pesaro e Urbino 
                          14                           23 
                      Ancona                     Macerata 
                          36                           24 
               Ascoli Piceno                        Fermo 
                           8                           12 
                     Viterbo                        Rieti 
                          18                            8 
                        Roma                       Latina 
                         316                           27 
                   Frosinone                     L'Aquila 
                          18                           13 
                      Teramo                      Pescara 
                          17                           32 
                      Chieti                   Campobasso 
                          26                           45 
                     Isernia                      Caserta 
                          25                           41 
                   Benevento                       Napoli 
                          13                          207 
                    Avellino                      Salerno 
                          19                           85 
                      Foggia                         Bari 
                          27                          115 
                     Taranto                     Brindisi 
                          34                           15 
                       Lecce        Barletta-Andria-Trani 
                          46                           29 
                     Potenza                       Matera 
                          27                           11 
                     Cosenza                    Catanzaro 
                          43                           29 
             Reggio Calabria                      Crotone 
                          42                            8 
               Vibo Valentia                      Trapani 
                           5                           29 
                     Palermo                      Messina 
                          77                           49 
                   Agrigento                Caltanissetta 
                          15                           17 
                        Enna                      Catania 
                           7                           86 
                      Ragusa                     Siracusa 
                          20                           22 
                     Sassari                        Nuoro 
                          20                            6 
                    Cagliari                     Oristano 
                          49                           14 
                Sud Sardegna 
                          21 

# Get original province labels
province_labels <- attr(phase$h_Province, "labels")
label_lookup <- setNames(names(province_labels), province_labels)

# Convert h_Province to character using labels
phase <- phase %>%
  mutate(h_Province_name = factor(h_Province, labels = names(province_labels)))

# Create a reverse lookup for group mapping
group_map <- list(
  "Teramo_pescara" = c(71, 72),
  "Matera_Potenza" = c(88, 87),
  "Crotone_Catanzaro_Vibo" = c(92, 90, 93),
  "Benevento_Avellino" = c(77, 79),
  "Modena_Reggio" = c(42, 41),
  "Udine_Gorizia" = c(31, 32),
  "Rieti_Viterbo" = c(66, 65),
  "Latina_Frosinone" = c(68, 69),
  "Genova_La_Spezia" = c(37, 38),
  "Savona_Imperia" = c(36, 35),
  "Sondrio_Bergamo" = c(12, 14),
  "Mantova_Cremona" = c(18, 17),
  "Lecco_Como" = c(19, 11),
  "Lodi_Pavia" = c(20, 16),
  "Ascoli_Piceno_Fermo" = c(63, 64),
  "Biella_Vercelli" = c(7, 2),
  "Verbano-Cusio-Ossola_Novara" = c(8, 3),
  "Brindisi_Lecce" = c(84, 85),
  "Nuoro_Sassari" = c(104, 103),
  "Oristano_Sud_Sardegna" = c(106, 107),
  "Enna_Caltanissetta" = c(99, 98),
  "Agrigento_Trapani" = c(97, 94),
  "Ragusa_Siracusa" = c(101, 102),
  "Arezzo_Siena" = c(54, 55),
  "Massa-Carrara_Lucca" = c(48, 49),
  "Pistoia_Prato" = c(50, 57),
  "Belluno_Treviso" = c(26, 27)
)

# Reverse group_map: value = group name, for each province code
group_lookup <- setNames(rep(names(group_map), lengths(group_map)), unlist(group_map))

# Create the new unified province name column
phase <- phase %>%
  mutate(h_Province_grouped = ifelse(
    h_Province %in% as.integer(names(group_lookup)),
    group_lookup[as.character(h_Province)],
    as.character(h_Province_name)
  ))


phase$h_Province_grouped = as.factor(phase$h_Province_grouped)
length(levels(phase$h_Province_grouped))

[1] 79

table(phase$h_Province_grouped)

           Agrigento_Trapani                  Alessandria 
                          44                           25 
                      Ancona                 Arezzo_Siena 
                          36                           37 
         Ascoli_Piceno_Fermo                         Asti 
                          20                           21 
                        Bari        Barletta-Andria-Trani 
                         115                           29 
             Belluno_Treviso           Benevento_Avellino 
                          60                           32 
             Biella_Vercelli                      Bologna 
                          23                           74 
               Bolzano/Bozen                      Brescia 
                          20                           72 
              Brindisi_Lecce                     Cagliari 
                          61                           49 
                  Campobasso                      Caserta 
                          45                           41 
                     Catania                       Chieti 
                          86                           26 
                     Cosenza       Crotone_Catanzaro_Vibo 
                          43                           42 
                       Cuneo           Enna_Caltanissetta 
                          31                           24 
                     Ferrara                      Firenze 
                          20                           72 
                      Foggia                 Forlì-Cesena 
                          27                           22 
            Genova_La_Spezia                     Grosseto 
                          76                           19 
                     Isernia                     L'Aquila 
                          25                           13 
            Latina_Frosinone                   Lecco_Como 
                          45                           39 
                     Livorno                   Lodi_Pavia 
                          25                           49 
                    Macerata              Mantova_Cremona 
                          24                           30 
         Massa-Carrara_Lucca               Matera_Potenza 
                          44                           38 
                     Messina                       Milano 
                          49                          292 
               Modena_Reggio        Monza e della Brianza 
                          75                           65 
                      Napoli                Nuoro_Sassari 
                         207                           26 
       Oristano_Sud_Sardegna                       Padova 
                          35                           84 
                     Palermo                        Parma 
                          77                           34 
                     Perugia              Pesaro e Urbino 
                          45                           23 
                    Piacenza                         Pisa 
                          21                           25 
               Pistoia_Prato                    Pordenone 
                          28                           17 
             Ragusa_Siracusa                      Ravenna 
                          42                           31 
             Reggio Calabria                Rieti_Viterbo 
                          42                           26 
                      Rimini                         Roma 
                          24                          316 
                      Rovigo                      Salerno 
                          19                           85 
              Savona_Imperia              Sondrio_Bergamo 
                          29                           69 
                     Taranto               Teramo_pescara 
                          34                           49 
                       Terni                       Torino 
                          14                          164 
                      Trento                      Trieste 
                          50                           25 
               Udine_Gorizia Valle d'Aosta/Vallée d'Aoste 
                          42                           28 
                      Varese                      Venezia 
                          57                           63 
 Verbano-Cusio-Ossola_Novara                       Verona 
                          27                           63 
                     Vicenza 
                          40 

which(table(phase$h_Province_grouped) < 30)

                 Alessandria          Ascoli_Piceno_Fermo 
                           2                            5 
                        Asti        Barletta-Andria-Trani 
                           6                            8 
             Biella_Vercelli                Bolzano/Bozen 
                          11                           13 
                      Chieti           Enna_Caltanissetta 
                          20                           24 
                     Ferrara                       Foggia 
                          25                           27 
                Forlì-Cesena                     Grosseto 
                          28                           30 
                     Isernia                     L'Aquila 
                          31                           32 
                     Livorno                     Macerata 
                          35                           37 
               Nuoro_Sassari              Pesaro e Urbino 
                          46                           52 
                    Piacenza                         Pisa 
                          53                           54 
               Pistoia_Prato                    Pordenone 
                          55                           56 
               Rieti_Viterbo                       Rimini 
                          60                           61 
                      Rovigo               Savona_Imperia 
                          63                           65 
                       Terni                      Trieste 
                          69                           72 
Valle d'Aosta/Vallée d'Aoste  Verbano-Cusio-Ossola_Novara 
                          74                           77 

which(table(phase$h_Province_grouped) < 20)

 Grosseto  L'Aquila Pordenone    Rovigo     Terni 
       30        32        56        63        69 

which(table(phase$h_Province_grouped) < 15)

L'Aquila    Terni 
      32       69 

Missing Values

They must be missing by design, because participants could skip ths question if they are vegetarian/vegan/pescetarian.

colSums(is.na(phase)) 

            record              Intro               Q169           Q169r3oe 
                 0                  0                  0                  0 
               age          Age_group               Q171         h_Province 
                 0                  0                  0                  0 
         h_Region2            h_Area2          education                 Q0 
                 0                  0                  0                  0 
              diet                 Q2                 Q3                 Q4 
                 0                  0                  0                219 
    meat_frequency                 Q6                 Q7                 Q8 
               241                694                241                530 
                Q9                Q10                Q11          Q12_Q14r1 
               241                388                241                241 
         Q12_Q14r2          Q12_Q14r3          Q15_Q17r1          Q15_Q17r2 
               241                241                241                241 
         Q15_Q17r3                Q18          Q19_Q29r1          Q19_Q29r2 
               241                241                241                241 
         Q19_Q29r3          Q19_Q29r4          Q19_Q29r5          Q19_Q29r6 
               241                241                241                241 
         Q19_Q29r7          Q19_Q29r8          Q19_Q29r9         Q19_Q29r10 
               241                241                241                241 
        Q19_Q29r11          Q30_Q32r1          Q30_Q32r2          Q30_Q32r3 
               241                241                241                241 
         Q33_Q35r1          Q33_Q35r2          Q33_Q35r3          Q36_Q47r1 
               241                241                241                  0 
         Q36_Q47r2          Q36_Q47r3          Q36_Q47r4          Q36_Q47r5 
                 0                  0                  0                  0 
         Q36_Q47r6          Q36_Q47r7          Q36_Q47r8          Q36_Q47r9 
                 0                  0                  0                  0 
        Q36_Q47r10         Q36_Q47r11         Q36_Q47r12          Q48_Q53r1 
                 0                  0                  0                  0 
         Q48_Q53r2          Q48_Q53r3          Q48_Q53r4          Q48_Q53r5 
                 0                  0                  0                  0 
         Q48_Q53r6          Q54_Q59r1          Q54_Q59r2          Q54_Q59r3 
                 0                  0                  0                  0 
         Q54_Q59r4          Q54_Q59r5          Q54_Q59r6          Q60_Q62r1 
                 0                  0                  0                219 
         Q60_Q62r2          Q60_Q62r3          Q63_Q65r1          Q63_Q65r2 
               219                219                219                219 
         Q63_Q65r3        Q66_Q69_1r1        Q66_Q69_1r2        Q66_Q69_1r3 
               219                219                219                219 
       Q66_Q69_1r4        Q66_Q69_1r5          Q70_Q90r1          Q70_Q90r2 
               219                219                  0                  0 
         Q70_Q90r3          Q70_Q90r4          Q70_Q90r5          Q70_Q90r6 
                 0                  0                  0                  0 
         Q70_Q90r7          Q70_Q90r8          Q70_Q90r9         Q70_Q90r10 
                 0                  0                  0                  0 
        Q70_Q90r11         Q70_Q90r12         Q70_Q90r13         Q70_Q90r14 
                 0                  0                  0                  0 
        Q70_Q90r15         Q70_Q90r16         Q70_Q90r17         Q70_Q90r18 
                 0                  0                  0                  0 
        Q70_Q90r19         Q70_Q90r20         Q70_Q90r21         Q91_Q104r1 
                 0                  0                  0                  0 
        Q91_Q104r2         Q91_Q104r3         Q91_Q104r4         Q91_Q104r5 
                 0                  0                  0                  0 
        Q91_Q104r6         Q91_Q104r7         Q91_Q104r8         Q91_Q104r9 
                 0                  0                  0                  0 
       Q91_Q104r10        Q91_Q104r11        Q91_Q104r12        Q91_Q104r13 
                 0                  0                  0                  0 
       Q91_Q104r14                edo        Q106_Q111r1        Q106_Q111r2 
                 0                  0                  0                  0 
       Q106_Q111r3        Q106_Q111r4        Q106_Q111r5        Q106_Q111r6 
                 0                  0                  0                  0 
       Q112_Q119r1        Q112_Q119r2        Q112_Q119r3        Q112_Q119r4 
                 0                  0                  0                  0 
       Q112_Q119r5        Q112_Q119r6        Q112_Q119r7        Q112_Q119r8 
                 0                  0                  0                  0 
           Q120pre        Q120_Q124r1        Q120_Q124r2        Q120_Q124r3 
                 0                  0                  0                  0 
       Q120_Q124r4        Q120_Q124r5               Q125               Q126 
                 0                  0                  0                138 
       Q127_Q131r1        Q127_Q131r2        Q127_Q131r3        Q127_Q131r4 
               138                138                138                138 
       Q127_Q131r5        Q132_Q137r1        Q132_Q137r2        Q132_Q137r3 
               138                138                138                138 
       Q132_Q137r4        Q132_Q137r5        Q132_Q137r6               Q138 
               138                138                138                138 
              Q139               Q140        Q141_Q145r1        Q141_Q145r2 
               138               3227               3227               3227 
       Q141_Q145r3        Q141_Q145r4        Q141_Q145r5        Q146_Q151r1 
              3227               3227               3227               3227 
       Q146_Q151r2        Q146_Q151r3        Q146_Q151r4        Q146_Q151r5 
              3227               3227               3227               3227 
       Q146_Q151r6               Q152               Q153               Q154 
              3227               3227                  0                  0 
       Q155_Q160r1        Q155_Q160r2        Q155_Q160r3        Q155_Q160r4 
                 0                  0                  0                  0 
       Q155_Q160r5        Q155_Q160r6    pol_orientation        Q162_Q168r1 
                 0                  0                  0                  0 
       Q162_Q168r2        Q162_Q168r3        Q162_Q168r4        Q162_Q168r5 
                 0                  0                  0                  0 
       Q162_Q168r6        Q162_Q168r7               Q173             gender 
                 0                  0                  0                  0 
     meat_identity     connect_nature    h_Province_name h_Province_grouped 
               241                  0                  0                  0 

table(phase$diet[is.na(phase$meat_identity)]) 

  1   2   3   4   5 
 11  11  73 116  30 

#219 non omnivores & 22 omniovres&flexitarians, strange!
#sanity check
table(is.na(phase$meat_frequency), is.na(phase$meat_identity))

       
        FALSE TRUE
  FALSE  3825    0
  TRUE      0  241

##very strange..
#let's investigate
#create missingness variable 0 = not missing 1= missing (success)
phase$missing <- ifelse(is.na(phase$meat_identity), 1, 0)
#log regression for missing meat identity
pred1 = glm(missing ~ age , data = phase, family = binomial)
pred2 = glm(missing ~ h_Province_grouped, data = phase, family = binomial) 
pred3 = glm(missing ~ gender, data = phase, family = binomial)
pred4 = glm(missing ~ edo, data = phase, family = binomial)
pred5 = glm(missing ~ connect_nature, data = phase, family = binomial)
pred6 = glm(missing ~ pol_orientation, data = phase, family = binomial)
pred7 = glm(missing ~ education, data = phase, family = binomial)

stargazer(pred1, pred2, pred3, pred4, pred5, pred6, pred7,
          type = "text",
          title = "Logistic Regression of Missingness",
          dep.var.labels = "Missing",
          digits = 3,
          star.cutoffs = c(0.1, 0.05, 0.01))

Logistic Regression of Missingness
======================================================================================================================
                                                                         Dependent variable:                          
                                               -----------------------------------------------------------------------
                                                                               Missing                                
                                                  (1)        (2)        (3)       (4)       (5)       (6)       (7)   
----------------------------------------------------------------------------------------------------------------------
age                                             -0.004                                                                
                                                (0.004)                                                               
                                                                                                                      
h_Province_groupedAlessandria                              16.574                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedAncona                                   15.011                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedArezzo_Siena                             14.983                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedAscoli_Piceno_Fermo                      0.00000                                                    
                                                         (1,759.025)                                                  
                                                                                                                      
h_Province_groupedAsti                                     15.570                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBari                                     16.215                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBarletta-Andria-Trani                    16.733                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBelluno_Treviso                          14.489                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBenevento_Avellino                       15.858                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBiella_Vercelli                          15.475                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBologna                                  16.307                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBolzano/Bozen                            16.831                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBrescia                                  15.011                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedBrindisi_Lecce                           15.909                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCagliari                                 15.409                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCampobasso                               15.927                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCaserta                                  15.596                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCatania                                  15.546                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedChieti                                   16.081                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCosenza                                  15.546                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCrotone_Catanzaro_Vibo                   15.570                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedCuneo                                    0.00000                                                    
                                                         (1,529.490)                                                  
                                                                                                                      
h_Province_groupedEnna_Caltanissetta                       16.168                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedFerrara                                  15.622                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedFirenze                                  15.733                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedFoggia                                   16.040                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedForlì-Cesena                             0.00000                                                    
                                                         (1,703.168)                                                  
                                                                                                                      
h_Province_groupedGenova_La_Spezia                         15.913                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedGrosseto                                 0.00000                                                    
                                                         (1,790.567)                                                  
                                                                                                                      
h_Province_groupedIsernia                                  0.00000                                                    
                                                         (1,633.622)                                                  
                                                                                                                      
h_Province_groupedL'Aquila                                 16.081                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedLatina_Frosinone                         15.498                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedLecco_Como                               15.648                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedLivorno                                  15.388                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedLodi_Pavia                               17.326                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMacerata                                 16.620                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMantova_Cremona                          15.199                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMassa-Carrara_Lucca                      15.522                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMatera_Potenza                           16.109                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMessina                                  16.597                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMilano                                   16.059                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedModena_Reggio                            15.690                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedMonza e della Brianza                    15.538                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedNapoli                                   15.475                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedNuoro_Sassari                            15.347                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedOristano_Sud_Sardegna                    15.040                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedPadova                                   15.270                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedPalermo                                  15.662                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedParma                                    15.070                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedPerugia                                  14.782                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedPesaro e Urbino                          15.475                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedPiacenza                                 0.00000                                                    
                                                         (1,729.992)                                                  
                                                                                                                      
h_Province_groupedPisa                                     16.124                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedPistoia_Prato                            0.00000                                                    
                                                         (1,576.828)                                                  
                                                                                                                      
h_Province_groupedPordenone                                15.793                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedRagusa_Siracusa                          16.565                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedRavenna                                  15.892                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedReggio Calabria                          16.001                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedRieti_Viterbo                            16.081                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedRimini                                   16.168                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedRoma                                     15.871                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedRovigo                                   0.00000                                                    
                                                         (1,790.567)                                                  
                                                                                                                      
h_Province_groupedSalerno                                  15.558                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedSavona_Imperia                           15.963                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedSondrio_Bergamo                          14.347                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedTaranto                                  15.070                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedTeramo_pescara                           15.409                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedTerni                                    0.00000                                                    
                                                         (2,001.460)                                                  
                                                                                                                      
h_Province_groupedTorino                                   16.027                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedTrento                                   17.180                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedTrieste                                  0.00000                                                    
                                                         (1,633.622)                                                  
                                                                                                                      
h_Province_groupedUdine_Gorizia                            16.315                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedValle d'Aosta/Vallée d'Aoste             15.270                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedVarese                                   15.252                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedVenezia                                  16.315                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedVerbano-Cusio-Ossola_Novara              0.00000                                                    
                                                         (1,594.573)                                                  
                                                                                                                      
h_Province_groupedVerona                                   15.875                                                     
                                                          (983.325)                                                   
                                                                                                                      
h_Province_groupedVicenza                                  15.622                                                     
                                                          (983.325)                                                   
                                                                                                                      
genderMen                                                            -0.718***                                        
                                                                      (0.142)                                         
                                                                                                                      
genderNon-binary                                                       1.380                                          
                                                                      (1.158)                                         
                                                                                                                      
genderPrefer not to say                                                0.533                                          
                                                                      (1.072)                                         
                                                                                                                      
edo                                                                            -0.228***                              
                                                                                (0.038)                               
                                                                                                                      
connect_nature                                                                           0.745***                     
                                                                                          (0.081)                     
                                                                                                                      
pol_orientation                                                                                    -0.215***          
                                                                                                    (0.048)           
                                                                                                                      
education                                                                                                     -0.014  
                                                                                                              (0.091) 
                                                                                                                      
Constant                                       -2.585***   -18.566   -2.479*** -2.033*** -6.519*** -1.960*** -2.739***
                                                (0.211)   (983.325)   (0.082)   (0.127)   (0.433)   (0.184)   (0.185) 
                                                                                                                      
----------------------------------------------------------------------------------------------------------------------
Observations                                     4,066      4,066      4,066     4,066     4,066     4,066     4,066  
Log Likelihood                                 -914.292   -863.735   -900.133  -895.057  -870.804  -904.490  -914.675 
Akaike Inf. Crit.                              1,832.583  1,885.470  1,808.266 1,794.114 1,745.607 1,812.980 1,833.350
======================================================================================================================
Note:                                                                                      *p<0.1; **p<0.05; ***p<0.01

####odds and prob
odds_to_prob <- function(model, coef_name) {
  coef_value <- coef(model)[coef_name]
  
  if (is.na(coef_value)) {
    stop("Coefficient not found in model.")
  }
  
  odds <- exp(coef_value)
  prob <- odds / (1 + odds)
  
  return(list(odds_ratio = odds, probability = prob))
}

results <- list(
  Age = odds_to_prob(pred1, "age"),
  Province = odds_to_prob(pred2, "h_Province_groupedAsti"),
  GenderMen = odds_to_prob(pred3, "genderMen"),
  EDO = odds_to_prob(pred4, "edo"),
  Connect_Nature = odds_to_prob(pred5, "connect_nature"),
  Pol_Orientation = odds_to_prob(pred6, "pol_orientation"),
  Education <- odds_to_prob(pred7, "education")
)

print(results)

$Age
$Age$odds_ratio
      age 
0.9964563 

$Age$probability
      age 
0.4991125 


$Province
$Province$odds_ratio
h_Province_groupedAsti 
               5782440 

$Province$probability
h_Province_groupedAsti 
             0.9999998 


$GenderMen
$GenderMen$odds_ratio
genderMen 
0.4876773 

$GenderMen$probability
genderMen 
0.3278112 


$EDO
$EDO$odds_ratio
      edo 
0.7963601 

$EDO$probability
      edo 
0.4433188 


$Connect_Nature
$Connect_Nature$odds_ratio
connect_nature 
      2.107455 

$Connect_Nature$probability
connect_nature 
     0.6781933 


$Pol_Orientation
$Pol_Orientation$odds_ratio
pol_orientation 
      0.8064374 

$Pol_Orientation$probability
pol_orientation 
      0.4464242 


[[7]]
[[7]]$odds_ratio
education 
0.9865534 

[[7]]$probability
education 
0.4966156 

#missing at random but making up a very small part of the data


sum(phase$missing)

[1] 241

phase = phase[phase$missing != 1, ]

Descriptive

#sample stats
#age (item)
table(phase$Q139)#born in the same region as they live in

   1    2 
3043  658 

sum(is.na(phase$Q139))#na = not born in Italy

[1] 124

table(phase$h_Province_grouped)

           Agrigento_Trapani                  Alessandria 
                          44                           22 
                      Ancona                 Arezzo_Siena 
                          35                           36 
         Ascoli_Piceno_Fermo                         Asti 
                          20                           20 
                        Bari        Barletta-Andria-Trani 
                         105                           25 
             Belluno_Treviso           Benevento_Avellino 
                          59                           30 
             Biella_Vercelli                      Bologna 
                          22                           67 
               Bolzano/Bozen                      Brescia 
                          17                           70 
              Brindisi_Lecce                     Cagliari 
                          57                           47 
                  Campobasso                      Caserta 
                          42                           39 
                     Catania                       Chieti 
                          82                           24 
                     Cosenza       Crotone_Catanzaro_Vibo 
                          41                           40 
                       Cuneo           Enna_Caltanissetta 
                          31                           22 
                     Ferrara                      Firenze 
                          19                           68 
                      Foggia                 Forlì-Cesena 
                          25                           22 
            Genova_La_Spezia                     Grosseto 
                          71                           19 
                     Isernia                     L'Aquila 
                          25                           12 
            Latina_Frosinone                   Lecco_Como 
                          43                           37 
                     Livorno                   Lodi_Pavia 
                          24                           38 
                    Macerata              Mantova_Cremona 
                          21                           29 
         Massa-Carrara_Lucca               Matera_Potenza 
                          42                           35 
                     Messina                       Milano 
                          43                          270 
               Modena_Reggio        Monza e della Brianza 
                          71                           62 
                      Napoli                Nuoro_Sassari 
                         198                           25 
       Oristano_Sud_Sardegna                       Padova 
                          34                           81 
                     Palermo                        Parma 
                          73                           33 
                     Perugia              Pesaro e Urbino 
                          44                           22 
                    Piacenza                         Pisa 
                          21                           23 
               Pistoia_Prato                    Pordenone 
                          28                           16 
             Ragusa_Siracusa                      Ravenna 
                          37                           29 
             Reggio Calabria                Rieti_Viterbo 
                          39                           24 
                      Rimini                         Roma 
                          22                          296 
                      Rovigo                      Salerno 
                          19                           81 
              Savona_Imperia              Sondrio_Bergamo 
                          27                           68 
                     Taranto               Teramo_pescara 
                          33                           47 
                       Terni                       Torino 
                          14                          152 
                      Trento                      Trieste 
                          40                           25 
               Udine_Gorizia Valle d'Aosta/Vallée d'Aoste 
                          38                           27 
                      Varese                      Venezia 
                          55                           57 
 Verbano-Cusio-Ossola_Novara                       Verona 
                          27                           59 
                     Vicenza 
                          38 

table(phase$h_Province_grouped, phase$gender)

                              
                               Women Men Non-binary Prefer not to say
  Agrigento_Trapani               26  18          0                 0
  Alessandria                     13   9          0                 0
  Ancona                          16  19          0                 0
  Arezzo_Siena                    18  17          0                 1
  Ascoli_Piceno_Fermo              9  11          0                 0
  Asti                             9  11          0                 0
  Bari                            51  54          0                 0
  Barletta-Andria-Trani           11  14          0                 0
  Belluno_Treviso                 34  25          0                 0
  Benevento_Avellino              18  12          0                 0
  Biella_Vercelli                 13   9          0                 0
  Bologna                         30  37          0                 0
  Bolzano/Bozen                    8   9          0                 0
  Brescia                         40  30          0                 0
  Brindisi_Lecce                  29  28          0                 0
  Cagliari                        25  22          0                 0
  Campobasso                      26  16          0                 0
  Caserta                         21  17          1                 0
  Catania                         36  46          0                 0
  Chieti                          13  11          0                 0
  Cosenza                         18  23          0                 0
  Crotone_Catanzaro_Vibo          19  21          0                 0
  Cuneo                           16  15          0                 0
  Enna_Caltanissetta              14   8          0                 0
  Ferrara                          8  11          0                 0
  Firenze                         37  31          0                 0
  Foggia                          14  11          0                 0
  Forlì-Cesena                     7  15          0                 0
  Genova_La_Spezia                36  35          0                 0
  Grosseto                        10   9          0                 0
  Isernia                          8  17          0                 0
  L'Aquila                         7   5          0                 0
  Latina_Frosinone                15  28          0                 0
  Lecco_Como                      18  19          0                 0
  Livorno                         10  14          0                 0
  Lodi_Pavia                      19  18          0                 1
  Macerata                        14   7          0                 0
  Mantova_Cremona                 15  14          0                 0
  Massa-Carrara_Lucca             23  19          0                 0
  Matera_Potenza                  16  19          0                 0
  Messina                         19  24          0                 0
  Milano                         126 144          0                 0
  Modena_Reggio                   33  38          0                 0
  Monza e della Brianza           30  32          0                 0
  Napoli                          99  98          1                 0
  Nuoro_Sassari                   13  12          0                 0
  Oristano_Sud_Sardegna           17  17          0                 0
  Padova                          42  39          0                 0
  Palermo                         35  37          0                 1
  Parma                           14  19          0                 0
  Perugia                         16  27          1                 0
  Pesaro e Urbino                 12  10          0                 0
  Piacenza                        10  11          0                 0
  Pisa                            12  11          0                 0
  Pistoia_Prato                   14  14          0                 0
  Pordenone                        5  11          0                 0
  Ragusa_Siracusa                 17  20          0                 0
  Ravenna                         15  14          0                 0
  Reggio Calabria                 19  20          0                 0
  Rieti_Viterbo                   15   9          0                 0
  Rimini                          13   9          0                 0
  Roma                           157 138          0                 1
  Rovigo                           7  12          0                 0
  Salerno                         39  42          0                 0
  Savona_Imperia                  13  14          0                 0
  Sondrio_Bergamo                 38  29          0                 1
  Taranto                         17  16          0                 0
  Teramo_pescara                  29  18          0                 0
  Terni                           10   4          0                 0
  Torino                          88  64          0                 0
  Trento                          19  20          0                 1
  Trieste                         15  10          0                 0
  Udine_Gorizia                   17  21          0                 0
  Valle d'Aosta/Vallée d'Aoste    13  14          0                 0
  Varese                          30  24          0                 1
  Venezia                         28  29          0                 0
  Verbano-Cusio-Ossola_Novara     18   9          0                 0
  Verona                          32  27          0                 0
  Vicenza                         16  22          0                 0

#tapply(phase$age, phase$h_Province_grouped, mean, na.rm = TRUE)
#tapply(phase$age, phase$h_Province_grouped, median, na.rm = TRUE)

Bivariate Descriptives

#level 2 province
boxplot(phase$meat_identity~ phase$h_Province, xlab = "Province", ylab = "Meat eater Identity")

boxplot(phase$edo~ phase$h_Province, xlab = "Province", ylab = "edo")

boxplot(phase$connect_nature~ phase$h_Province, xlab = "Province", ylab = "connectedness to nature")

#level1
##cormatrix
compute_lower_tri_corr = function(data) {
  cor_matrix = cor(data, use = "complete.obs")  # Compute correlation matrix
  cor_matrix[upper.tri(cor_matrix)] = NA  # Set upper triangle to NA
  return(cor_matrix)
}

cor_m = compute_lower_tri_corr(phase[,c("meat_identity","edo", "connect_nature")])
print("Full Sample Correlation (Lower Triangle):")

[1] "Full Sample Correlation (Lower Triangle):"

print(cor_m, na.print = "")#no paradox

               meat_identity      edo connect_nature
meat_identity      1.0000000                        
edo                0.1678507  1.00000               
connect_nature    -0.1575645 -0.22201              1

Multilevel Models Data centering (mean centering) and the Empty model

phase$edo_gmc = scale(phase$edo, center = TRUE, scale = FALSE) 

phase$connect_gmc = scale(phase$connect_nature, center = TRUE, scale = FALSE)

#Empty model
model0 = lme(
  meat_identity ~ 1,
  data = phase,
  method = "ML",
  na.action = "na.omit",
  random = ~1 | h_Province_grouped
)


summary(model0)

Linear mixed-effects model fit by maximum likelihood
  Data: phase 
       AIC      BIC    logLik
  14612.33 14631.08 -7303.167

Random effects:
 Formula: ~1 | h_Province_grouped
        (Intercept) Residual
StdDev:  0.09477886 1.630519

Fixed effects:  meat_identity ~ 1 
               Value Std.Error   DF  t-value p-value
(Intercept) 3.792183 0.0296878 3746 127.7354       0

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.80906569 -0.85617946 -0.03844385  0.74802629  2.00589518 

Number of Observations: 3825
Number of Groups: 79 

sigma0 <- VarCorr(model0);sigma0

h_Province_grouped = pdLogChol(1) 
            Variance    StdDev    
(Intercept) 0.008983032 0.09477886
Residual    2.658591736 1.63051885

var_u0 <- as.numeric(sigma0[1, "Variance"]);var_u0   # Between-group variance

[1] 0.008983032

var_e  <- as.numeric(sigma0[2, "Variance"]); var_e # Within-group variance

[1] 2.658592

# ICC
ICC <- var_u0 / (var_u0 + var_e)
ICC  # Intra-class correlation

[1] 0.00336749

Model 1- adding EDO as a level 1 predictor

model1 <- lme(
  meat_identity ~ edo_gmc,
  data = phase,
  method = "ML",
  na.action = "na.omit",
  random = ~1 |h_Province_grouped
)

summary(model1)

Linear mixed-effects model fit by maximum likelihood
  Data: phase 
       AIC      BIC    logLik
  14504.93 14529.93 -7248.467

Random effects:
 Formula: ~1 | h_Province_grouped
        (Intercept) Residual
StdDev:   0.0959436 1.607251

Fixed effects:  meat_identity ~ edo_gmc 
               Value  Std.Error   DF   t-value p-value
(Intercept) 3.792523 0.02941518 3745 128.93082       0
edo_gmc     0.143034 0.01358075 3745  10.53215       0
 Correlation: 
        (Intr)
edo_gmc 0.002 

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.0814017 -0.7939752  0.0121789  0.7661948  2.2605481 

Number of Observations: 3825
Number of Groups: 79 

sigma1 <- VarCorr(model1)
var1_u0 <- as.numeric(sigma1[1, "Variance"])
var1_e <- as.numeric(sigma1[2, "Variance"])

# R-squared
R2_level1 <- (var_e - var1_e) / var_e
R2_level2 <- (var_u0 - var1_u0) / var_u0
R2_total <- ((var_u0 + var_e) - (var1_u0 + var1_e)) / (var_u0 + var_e)

Mdoel 2

model2 <- lme(
  meat_identity ~ edo_gmc,
  data = phase,
  method = "ML",
  na.action = na.omit,
  random = ~1 + edo_gmc | h_Province_grouped,
  control = lmeControl(opt = "optim", maxIter = 100, msMaxIter = 100)
)

summary(model2)

Linear mixed-effects model fit by maximum likelihood
  Data: phase 
       AIC      BIC   logLik
  14508.46 14545.96 -7248.23

Random effects:
 Formula: ~1 + edo_gmc | h_Province_grouped
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev     Corr  
(Intercept) 0.10008718 (Intr)
edo_gmc     0.03895789 -0.294
Residual    1.60538365       

Fixed effects:  meat_identity ~ edo_gmc 
               Value  Std.Error   DF   t-value p-value
(Intercept) 3.791472 0.02963702 3745 127.93026       0
edo_gmc     0.144546 0.01470290 3745   9.83115       0
 Correlation: 
        (Intr)
edo_gmc -0.046

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.09324405 -0.79412174  0.01461931  0.76246745  2.27694524 

Number of Observations: 3825
Number of Groups: 79 

VarCorr(model2)

h_Province_grouped = pdLogChol(1 + edo_gmc) 
            Variance    StdDev     Corr  
(Intercept) 0.010017444 0.10008718 (Intr)
edo_gmc     0.001517717 0.03895789 -0.294
Residual    2.577256664 1.60538365       

Model 3

model3 <- lme(
  meat_identity ~ edo_gmc* connect_gmc,
  data = phase,
  method = "ML",
  na.action = "na.omit",
  random = ~1 + edo_gmc| h_Province_grouped,
  control = lmeControl(opt = "optim", maxIter = 100, msMaxIter = 100)
)

summary(model3)

Linear mixed-effects model fit by maximum likelihood
  Data: phase 
       AIC      BIC    logLik
  14451.71 14501.71 -7217.856

Random effects:
 Formula: ~1 + edo_gmc | h_Province_grouped
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev     Corr  
(Intercept) 0.10081309 (Intr)
edo_gmc     0.04510783 -0.474
Residual    1.59215567       

Fixed effects:  meat_identity ~ edo_gmc * connect_gmc 
                        Value  Std.Error   DF   t-value p-value
(Intercept)          3.793283 0.02994437 3743 126.67766  0.0000
edo_gmc              0.121399 0.01522746 3743   7.97235  0.0000
connect_gmc         -0.244429 0.03137277 3743  -7.79113  0.0000
edo_gmc:connect_gmc  0.004279 0.01490499 3743   0.28705  0.7741
 Correlation: 
                    (Intr) ed_gmc cnnct_
edo_gmc             -0.089              
connect_gmc          0.000  0.200       
edo_gmc:connect_gmc  0.179 -0.028  0.008

Standardized Within-Group Residuals:
          Min            Q1           Med            Q3           Max 
-2.1620976713 -0.7926821294 -0.0008223576  0.7585414647  2.5674098250 

Number of Observations: 3825
Number of Groups: 79 

VarCorr(model3)

h_Province_grouped = pdLogChol(1 + edo_gmc) 
            Variance    StdDev     Corr  
(Intercept) 0.010163280 0.10081309 (Intr)
edo_gmc     0.002034716 0.04510783 -0.474
Residual    2.534959663 1.59215567       

sigma3 <- VarCorr(model3)
var3_u0 <- as.numeric(sigma3[1, "Variance"])
var3_e <- as.numeric(sigma3[2, "Variance"])

Plotting the models (interactive plots just because it is fun and informative)

#plot Model0----
library(plotly)

#plot model1
# Predict for Model 1
phase$pred_model1 <- predict(model1, level = 1)

# Plot Model 1
onemodel =ggplot(phase, aes(x = edo_gmc, y = pred_model1, group = h_Province_grouped, color = as.factor(h_Province_grouped))) +
  geom_line(alpha = 0.6) +
  labs(
    x = "Ecological Dominance Orientation",
    y = "Predicted Meat Identity",
    color = "Province",
    title = "Model 1: Fixed Effect of Level 1 Predictor Edo"
  ) +
  theme_minimal()

p <- ggplotly(onemodel)
p  # questo forza la stampa se `ggplotly(onemodel)` non si mostra subito

# plot Model2-----
# Predict for Model 2
phase$pred_model2 <- predict(model2, level = 1)

# Plot Model 2
twomodel =ggplot(phase, aes(x = edo_gmc, y = pred_model2, group = h_Province_grouped, color = as.factor(h_Province_grouped))) +
  geom_line(alpha = 0.6) +
  labs(
    x = "Ecological Dominance Orientation",
    y = "Predicted Meat Identity",
    color = "Province",
    title = "Model 2: Random Slope Model with EDO at Level 1"
  ) +
  theme_minimal()

o = ggplotly(twomodel)
o

#plot Model3----
# Create grid of edo_gmc values for each province
new_data <- expand.grid(
  edo_gmc = seq(min(phase$edo_gmc, na.rm = TRUE), 
                max(phase$edo_gmc, na.rm = TRUE), 
                length.out = 100),
  h_Province_grouped = unique(phase$h_Province_grouped)
)

# Merge province-specific connect_gmc values
new_data <- merge(new_data, 
                  phase[!duplicated(phase$h_Province_grouped), 
                        c("h_Province_grouped", "connect_gmc")],
                  by = "h_Province_grouped")

# Generate predictions (including random effects)
new_data$pred <- predict(model3, newdata = new_data)

threemodel = ggplot(new_data, aes(x = edo_gmc, y = pred, 
                     group = h_Province_grouped,
                     color = connect_gmc)) +
  geom_line(alpha = 0.7) +
  scale_color_gradient2(name = "Connect Nature",
                        low = "darkblue", mid = "plum", high = "red",
                        midpoint = mean(new_data$connect_gmc)) +
  labs(x = "Ecological Dominance Orientation", 
       y = "Predicted Meat Identity",
       title = "Cross-Level Interaction: Province-Specific Slopes") +
  theme_minimal() +
  theme(legend.position = "bottom")

q = ggplotly(threemodel)
q

Assumption checks

# --- Level 1: Residual Diagnostics ---

# Predicted values
phase$predict <- predict(model3)

# Raw residuals
phase$resid <- phase$meat_identity - phase$predict

# Standardized residuals (Level 1)
phase$zresid <- scale(phase$resid)

# Histogram: Standardized residuals distribution (should be ~normal)
hist(phase$zresid, main = "Histogram of Level 1 Standardized Residuals",
     xlab = "Standardized Residuals")

# Residuals vs predictor (check for linearity & homoscedasticity)
plot(phase$edo_gmc, phase$zresid,
     main = "Level 1 Residuals vs edo_gmc",
     xlab = "edo_gmc", ylab = "Standardized Residuals")
abline(h = 0, col = "red", lty = 2)

# Residuals vs fitted values
plot(phase$predict, phase$zresid,
     main = "Level 1 Residuals vs Fitted Values",
     xlab = "Fitted Values", ylab = "Standardized Residuals")
abline(h = 0, col = "blue", lty = 2)

# --- Level 2: Random Effects Diagnostics ---

# Extract random effects (intercept + slope for edo_gmc)
randeff <- ranef(model3)  # Returns a data frame

# Optional: rename columns
colnames(randeff) <- c("u0", "u1")  # Intercept and slope

# Add grouping variable as a column
randeff$h_Province_grouped <- rownames(randeff)

# Merge with main data
phase <- merge(phase, randeff, by = "h_Province_grouped")

# Histograms
hist(phase$u0, main = "Random Intercepts (u0)",
     xlab = "Intercept (u0)", prob = TRUE)

hist(phase$u1, main = "Random Slopes (u1)",
     xlab = "Slope for edo_gmc (u1)", prob = TRUE)

# Random intercepts vs level-2 predictor (e.g., connect_gmc)
plot(phase$connect_gmc, phase$u0,
     main = "Random Intercepts vs connect_gmc",
     xlab = "connect_gmc", ylab = "u0")
abline(lm(u0 ~ connect_gmc, data = phase), col = "darkgreen")

# Random slopes vs level-2 predictor
plot(phase$connect_gmc, phase$u1,
     main = "Random Slopes vs connect_gmc",
     xlab = "connect_gmc", ylab = "u1")
abline(lm(u1 ~ connect_gmc, data = phase), col = "darkred")

Investigating possible reasons why the ICC is so low

############discussion
table(phase$Q125)

   1    2 
3701  124 

phase$Q141_Q145r1

<labelled<double>[3825]>: Q141_Q145r1: Mi sento parte della regione in cui vivo - Pensando alla regione in cui vivi, indica il tuo grado di accordo con le seguenti affermazioni, su una scala da 1= Completamente in disaccordo a 7= Completamente d’accordo.
   [1] NA NA NA NA NA NA NA  6 NA NA NA NA NA NA  6 NA NA NA NA NA NA NA NA NA
  [25] NA NA NA NA NA NA NA NA NA NA NA NA  4 NA NA NA NA NA NA NA NA  3 NA  5
  [49] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  2 NA NA  1 NA NA  5  6
  [73] NA NA NA NA NA NA  1 NA  4 NA  5 NA NA NA NA  7 NA NA NA NA NA NA NA  4
  [97] NA  6 NA NA NA  4  5 NA  4 NA NA NA  5 NA  2  4 NA  5  5 NA  4  7 NA NA
 [121]  3  7 NA NA NA  1 NA NA  5 NA NA NA  4 NA NA  5  5 NA NA  5 NA NA NA NA
 [145]  3 NA NA NA NA NA NA NA NA NA  5 NA NA NA NA NA NA NA  5 NA  6 NA NA NA
 [169] NA NA NA NA NA  4  5 NA  4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [193] NA NA NA  5 NA  5 NA NA NA NA NA NA NA NA NA NA NA  6 NA NA NA NA NA NA
 [217] NA NA NA  5 NA NA NA NA  4 NA NA  5 NA NA NA NA NA NA NA NA NA NA  7 NA
 [241] NA NA NA NA NA  7 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [265] NA NA NA NA NA NA NA NA NA NA NA NA NA  1 NA NA NA NA NA NA NA NA NA NA
 [289]  6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  6
 [313] NA NA NA NA NA  5 NA NA  7 NA NA NA  5 NA  7 NA NA NA NA NA NA  7 NA NA
 [337] NA NA NA NA NA  5 NA NA NA NA  7 NA NA  6 NA  7 NA NA NA NA NA NA NA NA
 [361] NA  5  7 NA NA NA NA NA NA NA NA NA  6 NA NA NA NA  1 NA NA NA NA NA NA
 [385] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  5
 [409] NA NA NA NA  7 NA NA NA NA NA  5 NA  7 NA  5 NA  2 NA NA NA  4 NA NA  4
 [433] NA NA  7 NA NA NA NA NA NA NA  3 NA NA NA NA NA NA NA NA NA NA NA NA NA
 [457] NA NA NA  5 NA NA  5 NA  3  6 NA  7 NA  5 NA NA NA NA  4  5 NA  4  1  4
 [481] NA  5  5 NA NA NA NA  7 NA NA NA NA NA NA NA NA NA NA NA NA  6  3 NA NA
 [505] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [529] NA NA NA NA NA NA  5 NA NA NA  5  4 NA NA NA NA NA NA NA NA NA NA  2 NA
 [553] NA NA NA NA NA NA NA  7 NA NA  5 NA NA NA NA NA NA NA NA NA NA NA NA NA
 [577] NA NA  6 NA  3 NA NA NA NA NA NA NA NA  5 NA  7 NA  5 NA NA NA NA NA  1
 [601]  6  5 NA NA NA NA  6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  5  6
 [625] NA NA NA  5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [649] NA  5 NA NA  5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [673] NA  7 NA NA  3  5 NA  6 NA  5  6 NA  7  3 NA  5  6 NA NA NA NA NA NA NA
 [697] NA  4 NA NA NA NA NA NA NA NA NA NA NA NA  3  1 NA  4 NA NA  5 NA NA NA
 [721] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [745] NA  5 NA NA NA NA NA  5 NA NA NA NA NA NA NA  4 NA NA NA NA NA NA  7 NA
 [769] NA  6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [793] NA NA NA NA NA NA NA NA NA NA NA NA  1 NA  6 NA NA  2 NA NA NA NA NA NA
 [817] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [841] NA NA NA NA NA NA NA NA  4  1 NA NA NA NA  6 NA NA NA NA  7 NA NA NA NA
 [865] NA NA NA NA  6 NA NA NA NA NA  7  6 NA NA NA NA NA NA NA NA NA NA NA NA
 [889] NA NA NA NA NA NA  4 NA NA NA  4 NA  5 NA NA NA NA NA NA NA NA NA NA NA
 [913] NA NA NA NA NA NA NA NA NA NA NA  5  1 NA NA  6 NA NA  5 NA  3 NA  7 NA
 [937] NA NA NA NA NA NA  2  7 NA NA NA NA  2 NA  5 NA NA NA  5 NA  4 NA NA NA
 [961] NA NA NA NA  3  7 NA NA NA  5 NA NA NA NA NA NA NA NA NA  7 NA NA  5  4
 [985] NA NA NA NA NA NA NA NA NA  7 NA NA NA NA NA  5 NA NA NA NA NA NA NA NA
[1009] NA NA  6  3  5 NA NA NA NA NA NA  6 NA NA  6  4 NA NA  5 NA  4 NA NA NA
[1033] NA NA NA NA NA  6 NA  5 NA  4 NA  7  5  2 NA NA NA  4 NA  5  5 NA  2  5
[1057]  5 NA  6  6  4 NA NA NA NA NA NA  2 NA  5 NA NA NA  4 NA NA NA NA NA NA
[1081] NA NA NA NA NA NA  3 NA NA NA NA NA NA NA NA  2 NA NA NA NA NA NA NA NA
[1105] NA NA NA NA NA  3 NA NA NA NA NA NA NA NA NA NA NA NA  6  7 NA NA  4  2
[1129] NA  7  6  6 NA NA NA NA NA NA NA NA NA NA  6 NA NA NA NA NA NA NA NA NA
[1153] NA NA NA NA NA NA NA NA NA  3  7 NA NA NA NA NA  5 NA NA NA NA NA  7 NA
[1177] NA NA NA NA NA NA NA NA NA NA NA  4 NA  5 NA NA  5 NA  6 NA NA NA NA  2
[1201] NA NA NA NA NA  5  7  4 NA NA  5  3  4 NA NA NA NA  6 NA  6 NA NA  3  7
[1225] NA NA NA NA NA  4 NA NA NA NA NA NA NA NA NA  3  5  5 NA NA NA  3 NA NA
[1249] NA NA  4 NA NA  4 NA  4  1 NA NA NA NA NA  5 NA NA NA NA NA NA NA NA NA
[1273] NA NA  7 NA NA NA NA NA NA  5 NA  5 NA  4 NA NA NA NA NA NA NA NA NA  5
[1297] NA NA NA  5 NA NA NA NA NA NA  4 NA  6 NA NA  4 NA NA NA NA  7 NA  7  3
[1321] NA NA  5 NA NA NA NA NA  4 NA NA NA NA NA  5  4  4  5 NA NA NA NA NA NA
[1345] NA NA  1 NA NA  1  5  6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[1369] NA  4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA  7 NA NA  5  3  4 NA  6
[1393]  1  6 NA NA  1 NA  4  5 NA NA NA NA  5 NA NA  7 NA NA NA  5 NA NA NA  4
[1417] NA NA NA NA NA NA  3  5  3 NA NA NA NA NA NA NA NA NA  7  4  5 NA NA NA
[1441] NA NA NA NA NA NA  4 NA NA NA  5 NA NA NA  5 NA NA NA NA NA NA NA NA  6
[1465]  5 NA NA NA  5 NA NA NA NA NA NA  6 NA NA  5 NA NA NA NA NA  5  5 NA NA
[1489] NA  5 NA NA  7 NA NA NA  5  7  7 NA NA  6  7 NA NA NA NA  7 NA NA NA  5
[1513] NA NA NA NA NA NA NA NA NA NA  5 NA NA NA NA  4  7 NA NA NA NA NA NA NA
[1537] NA NA  3 NA NA NA NA NA  4 NA NA NA NA NA NA  2 NA NA NA  4 NA NA NA NA
[1561]  5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  4  4  4 NA
[1585] NA  7 NA NA  5  7 NA  6  6 NA NA  4  5 NA NA NA  6 NA NA NA  5 NA NA  7
[1609] NA  1  5 NA NA NA  6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA  2 NA  6
[1633] NA NA  4 NA NA NA NA NA  4 NA  7  4 NA  3 NA NA NA NA NA NA NA  5 NA  5
[1657] NA NA  5  5 NA NA NA NA NA  6 NA  6 NA NA NA NA NA  4  6  4 NA NA  5 NA
[1681]  5 NA NA NA NA  6  5 NA  4 NA NA  5  4 NA NA NA  4 NA NA NA NA NA  4 NA
[1705] NA  5 NA NA NA  5 NA NA  5 NA  5  4 NA NA  5  4 NA NA NA NA  6 NA NA NA
[1729] NA NA  6 NA  3 NA NA  5 NA  5  7 NA NA NA NA NA  3 NA NA NA NA NA NA NA
[1753]  1 NA  5 NA NA  4 NA NA  4 NA NA NA  4 NA  4 NA NA  5  7  4  4 NA NA NA
[1777] NA NA NA NA NA NA  7 NA  6 NA NA NA NA  5 NA NA NA NA NA NA  5  4 NA  4
[1801] NA  2 NA NA NA NA NA NA  5  6  2  5 NA NA NA NA NA  5 NA NA NA  3 NA NA
[1825] NA  5  4 NA NA  5 NA NA NA  5  7 NA NA NA NA NA  4 NA NA NA NA  3 NA NA
[1849]  4 NA NA  4  4 NA NA NA  4  2 NA  5 NA  5 NA  5 NA NA  3 NA NA NA  6 NA
[1873] NA NA NA NA  5 NA  4 NA NA NA NA  5 NA  4  1 NA NA  6 NA NA NA NA NA  5
[1897] NA NA NA  5 NA  5 NA NA NA NA NA NA NA NA NA NA  6 NA NA NA NA NA  5 NA
[1921] NA NA NA  5 NA NA  4 NA  6 NA NA  5 NA  7 NA NA NA NA  4  5 NA NA  4 NA
[1945]  4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  4 NA NA NA
[1969] NA NA NA NA NA NA NA NA NA NA  5  5 NA NA NA  5  6 NA NA NA NA  4 NA  5
[1993] NA NA  5 NA NA NA NA NA NA NA NA NA NA NA  5 NA NA NA NA NA NA NA NA NA
[2017] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[2041] NA NA NA NA NA NA NA NA NA NA NA NA NA NA  5  7 NA NA NA NA NA NA NA NA
[2065] NA NA NA NA NA NA NA NA NA NA NA  5 NA NA NA NA NA NA NA NA NA NA NA NA
[2089] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  4
[2113] NA NA NA NA NA NA NA NA NA NA NA NA NA NA  6 NA NA NA NA NA NA NA NA NA
[2137] NA NA NA NA NA NA NA NA NA NA NA NA  3 NA NA NA NA NA NA NA NA NA NA NA
[2161] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  5 NA NA  6 NA NA NA NA NA
[2185]  6 NA  7  7 NA NA  1 NA NA NA NA NA NA  4 NA NA NA NA NA NA NA NA NA NA
[2209] NA NA NA NA NA NA NA NA NA NA  5 NA  5 NA NA NA NA  7 NA NA NA NA NA  4
[2233]  1 NA NA NA NA  3 NA  5  2 NA NA NA NA NA NA  5 NA NA NA NA NA NA NA  4
[2257] NA NA NA NA  7  5 NA NA NA NA NA NA  1 NA  4 NA  7 NA NA NA NA  4 NA NA
[2281] NA NA NA  5 NA  5 NA NA NA  5 NA  5 NA NA  4 NA NA NA NA  4 NA NA NA  4
[2305] NA  4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[2329] NA NA NA NA NA NA  5 NA NA NA NA  4 NA NA  5 NA NA NA NA NA  7 NA NA NA
[2353]  5 NA NA  2 NA NA NA NA NA NA NA NA  4 NA NA NA NA NA NA NA NA NA NA NA
[2377] NA NA NA NA NA NA NA NA NA  6 NA NA  6 NA NA NA NA  4 NA NA  4 NA  7  5
[2401]  5 NA  5  6 NA  4 NA  5 NA NA NA  5  5 NA  4 NA  6 NA NA NA NA NA NA NA
[2425] NA NA NA  6 NA NA NA NA  4  5 NA NA NA NA NA NA  7 NA NA NA NA NA  5  3
[2449] NA NA NA NA  6  3 NA NA NA NA  4 NA NA  5 NA NA NA NA NA NA NA NA NA NA
[2473]  4 NA NA NA NA  7 NA NA  5 NA  3  6 NA NA NA NA NA NA  6 NA  5  3  4 NA
[2497] NA  5 NA NA NA  5 NA NA NA NA NA NA NA  5 NA NA  5 NA NA NA  6 NA NA  4
[2521]  6  3 NA NA NA NA NA NA NA NA NA NA  6 NA NA NA NA NA  7 NA NA NA NA NA
[2545] NA NA NA  5 NA NA NA NA NA NA NA NA NA NA NA  5 NA  1 NA NA NA  7 NA NA
[2569] NA NA NA NA NA  6 NA NA NA NA NA NA NA NA NA NA NA NA  7 NA  7 NA NA NA
[2593] NA  5 NA NA NA NA NA NA NA NA NA  6 NA NA  5 NA NA NA NA NA NA NA NA NA
[2617] NA NA NA  4 NA NA NA NA NA  6 NA NA NA NA NA  7 NA NA NA NA NA NA NA NA
[2641] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[2665] NA NA NA NA  4  5 NA NA NA NA NA NA  4 NA NA NA NA  3 NA  5  3 NA NA NA
[2689] NA NA NA NA NA NA NA NA NA NA  4  5  5 NA NA  3  6 NA  6 NA  4 NA NA NA
[2713] NA NA NA NA  6 NA NA  6 NA  5 NA  3  6 NA NA NA  5  3  2 NA  5 NA NA  6
[2737] NA NA  3 NA  6 NA NA  6 NA NA NA NA NA NA  7 NA NA  6 NA NA NA NA NA NA
[2761] NA NA NA NA NA NA  1 NA NA NA  6  4 NA NA NA NA  5 NA NA  4 NA NA  4 NA
[2785] NA NA NA NA NA NA NA  4  7  4 NA NA  7  5 NA NA NA NA NA NA NA  4 NA NA
[2809] NA NA NA NA NA NA  4 NA  5 NA NA NA NA  4  5  4 NA NA NA NA NA  6 NA NA
[2833] NA NA NA NA NA NA NA NA NA NA NA  1 NA  4  6 NA  4 NA NA NA  4 NA NA NA
[2857] NA NA NA  5 NA NA NA NA  6 NA NA NA NA NA NA NA NA NA NA NA NA  5  5 NA
[2881]  5 NA NA  4 NA NA NA  6 NA NA NA  5 NA  6  6 NA  4  2 NA NA NA NA NA NA
[2905] NA NA NA NA NA  5 NA NA  3  7 NA NA NA  1 NA NA NA NA  4 NA NA NA NA NA
[2929] NA NA  5  2 NA NA NA  1 NA NA NA NA  5 NA NA  5 NA NA NA NA NA NA NA NA
[2953] NA NA  6  6 NA NA NA NA NA NA  5 NA  3  4 NA  5 NA  4  4 NA NA NA NA NA
[2977] NA  7 NA  7  4 NA  3 NA NA NA NA NA NA NA  5  7 NA NA NA NA NA NA  5 NA
[3001] NA  7 NA NA NA NA NA NA NA NA NA NA NA NA NA  5 NA NA NA NA NA NA NA  4
[3025] NA NA NA NA NA NA NA NA NA  6 NA NA NA NA NA NA NA NA NA NA NA NA NA  6
[3049]  5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[3073] NA NA NA NA NA NA NA NA NA NA NA NA  5  6 NA NA  5 NA NA NA NA NA NA NA
[3097] NA  5  5 NA  3 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[3121] NA NA NA  4 NA NA NA  5 NA NA NA NA NA  3  6  6 NA NA  4 NA NA NA NA NA
[3145] NA  4 NA  5 NA NA NA NA NA  7 NA NA NA NA  1 NA NA  5 NA NA NA NA  6 NA
[3169] NA NA NA NA NA NA NA NA NA NA NA NA NA  4 NA NA  4 NA NA  5 NA  3 NA NA
[3193]  1 NA  4 NA NA NA NA NA NA NA NA NA NA  4 NA  5 NA NA NA NA NA NA NA NA
[3217] NA NA NA NA  1 NA NA NA NA NA NA NA NA NA  5 NA NA NA NA NA  6 NA NA NA
[3241] NA NA NA  6 NA NA NA NA NA NA  7  5 NA  5 NA NA NA NA NA NA NA NA NA  6
[3265] NA  5 NA NA  4  4 NA NA NA  5 NA  3  7 NA NA NA  5 NA NA  5 NA  6  5 NA
[3289] NA NA NA  5 NA NA NA NA NA NA NA  3  5 NA  4 NA NA NA NA  6 NA NA NA NA
[3313]  2 NA NA NA NA NA  4 NA NA NA NA NA NA  5 NA NA NA NA  6 NA NA NA  7 NA
[3337] NA NA NA NA NA NA NA NA  4 NA NA NA NA NA NA  5  4 NA NA NA NA NA NA NA
[3361] NA NA NA NA NA NA NA NA NA  5 NA NA NA NA  5  5  5 NA NA NA NA NA NA  6
[3385] NA  5 NA NA  3  4 NA NA NA NA NA NA NA  3  5  1 NA NA NA NA NA  4  5  7
[3409]  5 NA  4  5  4 NA NA NA NA NA NA NA NA NA  4 NA NA NA NA NA  5 NA NA NA
[3433] NA NA NA NA NA NA NA NA  5 NA NA NA NA NA NA NA NA  6 NA NA  6 NA NA NA
[3457]  2 NA NA NA  4 NA NA NA NA NA NA NA NA NA  5 NA  7 NA  2  6 NA  6  7 NA
[3481] NA  4  5 NA  6 NA  3  5 NA NA  6 NA NA NA NA NA NA NA  5 NA  5 NA  6 NA
[3505] NA NA  5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA  3 NA NA NA NA NA NA
[3529] NA NA NA NA NA  7 NA  5 NA NA NA NA NA NA NA  5 NA NA NA NA  5 NA NA NA
[3553]  4 NA NA  6 NA NA NA  1 NA  1  5 NA NA  6 NA  2 NA NA NA NA NA NA  4  6
[3577]  5 NA NA NA NA  4 NA NA  4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[3601] NA NA NA  4 NA NA  6  1 NA  4  5  4 NA NA  4 NA NA NA NA NA NA NA  5 NA
[3625]  7 NA NA NA NA NA NA NA NA NA NA  4 NA  5 NA NA NA  7 NA  5 NA NA NA  6
[3649] NA NA NA  7 NA NA NA NA NA NA NA NA NA NA NA NA NA  6 NA NA NA  6 NA NA
[3673] NA  3 NA NA NA NA NA NA NA NA NA NA NA NA  4  6 NA NA NA NA NA NA NA NA
[3697] NA NA  4 NA NA NA NA  3 NA NA NA  5 NA NA NA NA NA NA NA NA NA NA NA NA
[3721] NA NA NA NA  6 NA NA  4  6 NA  3 NA NA NA NA NA  5 NA NA NA NA  4 NA NA
[3745] NA NA NA  4  4 NA  3 NA NA  3 NA NA NA NA  5  4 NA NA  5  1 NA  1 NA NA
[3769] NA NA NA NA NA  5  5  1 NA NA NA NA NA NA NA NA  5 NA NA NA  5 NA NA NA
[3793] NA NA NA NA NA  6 NA  5 NA NA NA NA NA NA NA  5 NA NA  5 NA NA NA NA NA
[3817] NA NA NA NA NA NA  5  4 NA

Labels:
 value                          label
     1   1Completamente in disaccordo
     2           2Molto in disaccordo
     3      3Abbastanza in disaccordo
     4 4Né d’accordo né in disaccordo
     5          5Abbastanza d’accordo
     6               6Molto d’accordo
     7       7Completamente d’accordo

phase$legame_nato = rowMeans(phase[, c("Q127_Q131r1", "Q127_Q131r2", "Q127_Q131r3", "Q127_Q131r4", "Q127_Q131r5")])
phase$legame_abita = rowMeans(phase[, c("Q141_Q145r1", "Q141_Q145r2", "Q141_Q145r3", "Q141_Q145r4", "Q141_Q145r5")])
#how good is the connection
# Count of scores 4 and below
sum(phase$legame_nato <= 4, na.rm = TRUE)

[1] 859

sum(phase$legame_abita <= 4, na.rm = TRUE)

[1] 233

sum(phase$legame_nato > 4, na.rm = TRUE)

[1] 2842

sum(phase$legame_abita > 4, na.rm = TRUE)

[1] 549

sum(phase$legame_nato >= 6, na.rm = TRUE)

[1] 833

sum(phase$legame_abita >= 6, na.rm = TRUE)

[1] 114