---
title: "Uni-Impresa 'You Drink You Lose' - Analisi dei dati Pre vs. Post & Pre vs. Post vs. Follow-up"
author:
- name: Enrico Rubaltelli
  affiliation: "DPSS - Università di Padova"
- name: Marta Caserotti
  affiliation: "DPSS - Università di Padova"
- name: Beatrice Conte
  affiliation: "DPSS - Università di Padova"
- name: Lorella Lotto
  affiliation: "DPSS - Università di Padova"
date: "`r Sys.Date()`"
output:
  html_document:
    include:
      before_body: Header.html
    css: Style.css
    toc: yes
    toc_float: no
    fig_height: 8
    fig_width: 8
    number_sections: yes
    code_folding: hide
    df_print: paged
    theme: journal
  word_document:
    toc: yes
header-includes:
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
- \usepackage{rotating}
- \usepackage{xcolor}
- \usepackage{longtable}
- \usepackage{gensymb}
- \usepackage{amsfonts}
- \usepackage{amsmath}
- \pretitle{\begin{left} \includegraphics[width=2in,height=2in]{logo.jpg}\LARGE\\}
fontsize: 10pt
always_allow_html: yes
editor_options:
  markdown:
    wrap: sentence
---

```{r knitr_init, echo=FALSE, cache=FALSE}
library(knitr)
library(rmdformats)

## Global options
options(scipen=999)
options(digits=8)
options(round=6)
options(max.print=10000)
opts_chunk$set(echo=FALSE,
               cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE)
opts_knit$set(width=75)
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# ***INTRODUZIONE E SUMMARY RISULTATI PRINCIPALI***

# ***STATISTICHE DESCRITTIVE***

```{r td1, echo=TRUE}
library(tidyverse) 
library(sandwich)
library(gtsummary)
library(flextable)
library(tables)
#library(openxlsx)

#Change directory to folder where folder with files is located
#library(readbulk)
#dsSIM <- read_bulk(directory = "csv_guide",
#    extension = ".csv")

#write.table(dsSIM, file = "dsSIM.csv", sep = ",",
#            na = " ", dec = ".", row.names = FALSE,
#            col.names = FALSE)

## Upload of dataset from PRE-TEST
dsPRE <- rio::import(file = "UDUL_Questionario+PRE_January+22,+2024_18.38.xlsx",
                     which = 1L, skip = 1L)

## I prepare the dataset from the PRE-TEST for the analyses
## Delete incomplete responses
dsPRE_ok <- subset(dsPRE, Finished1 == 1)
## Delete participants with multiple complete responses
dsPRE_ok <- subset(dsPRE_ok, Code3 != "en710os")
dsPRE_ok <- subset(dsPRE_ok, Code3 != "gr383ot")
dsPRE_ok <- subset(dsPRE_ok, Code3 != "go020be")
dsPRE_ok <- subset(dsPRE_ok, Code3 != "jo399io")
dsPRE_ok <- subset(dsPRE_ok, Code3 != "va203on")
dsPRE_ok <- subset(dsPRE_ok, Code3 != "en877zo")
dsPRE_ok <- subset(dsPRE_ok, Code3 != "pa753an")
dsPRE_ok1 <- dsPRE_ok%>%select(!c(Code1,Code21))

## Upload of dataset from POST-TEST
dsPOST <- rio::import(file = "UDUL_Questionario+POST_January+22,+2024_19.38.xlsx",
                      which = 1L, skip = 1L)

## I prepare the dataset from the PRE-TEST for the analyses
## Delete incomplete responses
dsPOST_ok <- subset(dsPOST, Finished2 == 1)
dsPOST_ok <- subset(dsPOST_ok, Code3 != "va203on")
dsPOST_ok1 <- dsPOST_ok%>%select(!c(Code2,Code22))

## Upload of dataset with driving conditions (under the influence vs. control)
dsSIM <- read.table(file = "dsSIM.csv", sep = ",", dec = ".", header = TRUE)


## Select only Code and Test
#dsSIM1 <- dsSIM %>% select(Code3,Test)

## Create full datset
dsPRE_POST <- left_join(dsPRE_ok1,dsPOST_ok1, by = "Code3")
dsPRE_POST_ok <- left_join(dsPRE_POST,dsSIM, by = "Code3")
## Delete training sessions
dsPRE_POST_ok <- subset(dsPRE_POST_ok, Test != "Training")
```

```{r td3, echo=TRUE, results='markup', fig.height=5, fig.width=8}
library(tidyr)
library(ggplot2)
library(ggthemes)
## Tranform Gender into a factor with labels
dsPRE_POST_ok$Genere <- factor(dsPRE_POST_ok$Genere1, levels = c(1,2,3,4),
                       labels = c("Uomo", "Donna", 
                                  "Altro","Preferisco non risponder"))

## Check that Age is 18 to 24 years
#summary(dsPRE_POST_ok$Età1)
### Delete participants above 24 years
dsPRE_POST_ok <- subset(dsPRE_POST_ok, Età1 <= 24)
dsPRE_POST_ok$Età <- dsPRE_POST_ok$Età1

## Transform Education into a factor with labels
dsPRE_POST_ok$Istruzione <- factor(dsPRE_POST_ok$Istr1, levels = c(1,2,3,4),
                       labels = c("Scuola media", "Scuola superiore", 
                                  "Laurea triennale", "Laurea magistrale"))

## Transform Work into a factor with labels
dsPRE_POST_ok$Lavoro <- factor(dsPRE_POST_ok$Profes1, levels = c(1,2,3,4),
                       labels = c("Studente", "Occupato",
                                  "Disoccupato", "In cerca di occupazione"))

## Transform Residence into a factor with labels
dsPRE_POST_ok$Residenza <- factor(dsPRE_POST_ok$Resid1, levels = c(1,2),
                       labels = c("Comune PD", "Provincia PD"))

## Clean answers to the question about the place where they live
dsPRE_POST_ok$Comune <- dsPRE_POST_ok$Resid2

## Transform Driving licence into a factor with labels
### CHECK WHETHER THE LABELS ARE CODED CORRECTLY 
### It looks strange that there are so many 4!
dsPRE_POST_ok$Patente <- factor(dsPRE_POST_ok$Patente1, levels = c(1,2,3,4,5),
                       labels = c("No patente", "A", "B1", "B", "Altro"))

## Transform the first question on Fines into a factor with labels
dsPRE_POST_ok$Sanzioni_alcool <- factor(dsPRE_POST_ok$sanz11, levels = c(1,2,3,4),
                       labels = c("0,00 g/l", "0,2 g/l", 
                                  "0,5 g/l", "0,8 g/l"))

## Transform the second question on Fines into a factor with labels
dsPRE_POST_ok$Sanzioni_neo_patentati <- factor(dsPRE_POST_ok$sanz21, levels = c(1,2,3),
                       labels = c("Multa €164-664 + 10 punti", 
                                  "Multa €708-2883.34 + 20 punti", 
                                  "Nessuna sanzione"))

## Recode to have the correct answer vs. the incorrect ones
dsPRE_POST_ok$Sanzioni_neo_patentati2 <- car::recode(dsPRE_POST_ok$sanz21,
                                                     "1=1;2:3=0")


## Compute average score for Risk Perception in the before the driving test
dsPRE_POST_ok$Percezione_rischio <- rowMeans(dsPRE_POST_ok[,c("RP_scale_11",
                                                              "RP_scale_21",
                                                              "RP_scale_31",
                                                              "RP_scale_41",
                                                              "RP_scale_51")],
                                             na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFArisk <- cronbach(cbind(dsPRE_POST_ok[,c("RP_scale_11",
                                               "RP_scale_21",
                                               "RP_scale_31",
                                               "RP_scale_41",
                                               "RP_scale_51")]))
ALFArisk

## Compute average score for Vankov scale in the before the driving test
dsPRE_POST_ok$Vankov <- rowMeans(dsPRE_POST_ok[,c("DUI_int1",
                                                  "Inst_att1",
                                                  "Per_cont1",
                                                  "Self_eff1",
                                                  "Des_norm1",
                                                  "Sub_norm1",
                                                  "Aff_att1")],
                                             na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAvankov <- cronbach(cbind(dsPRE_POST_ok[,c("DUI_int1",
                                                  "Inst_att1",
                                                  "Per_cont1",
                                                  "Self_eff1",
                                                  "Des_norm1",
                                                  "Sub_norm1",
                                                  "Aff_att1")]))
ALFAvankov

## Compute average score for Dejoy scale (involvement) in the before the driving test
dsPRE_POST_ok$Dejoy_coinvolgimento <- rowMeans(dsPRE_POST_ok[,c(
  "involv_DJ_11","involv_DJ_21",
  "involv_DJ_31","involv_DJ_41",
  "involv_DJ_51","involv_DJ_61",
  "involv_DJ_71","involv_DJ_81",
  "involv_DJ_91","involv_DJ_101")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyC <- cronbach(cbind(dsPRE_POST_ok[,c(
  "involv_DJ_11","involv_DJ_21",
  "involv_DJ_31","involv_DJ_41",
  "involv_DJ_51","involv_DJ_61",
  "involv_DJ_71","involv_DJ_81",
  "involv_DJ_91","involv_DJ_101")]))
ALFAdejoyC

## Compute average score for Dejoy scale (gravity) in the before the driving test
dsPRE_POST_ok$Dejoy_gravità <- rowMeans(dsPRE_POST_ok[,c(
  "grav_DJ_11","grav_DJ_21","grav_DJ_31","grav_DJ_41",
  "grav_DJ_51","grav_DJ_61","grav_DJ_71","grav_DJ_81",
  "grav_DJ_91","grav_DJ_101")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyG <- cronbach(cbind(dsPRE_POST_ok[,c(
  "grav_DJ_11","grav_DJ_21","grav_DJ_31","grav_DJ_41",
  "grav_DJ_51","grav_DJ_61","grav_DJ_71","grav_DJ_81",
  "grav_DJ_91","grav_DJ_101")]))
ALFAdejoyG

## Compute average score for Dejoy scale (control) in the before the driving test
dsPRE_POST_ok$Dejoy_controllo <- rowMeans(dsPRE_POST_ok[,c(
  "contol_DJ_11","contol_DJ_21","contol_DJ_31","contol_DJ_41",
  "contol_DJ_51","contol_DJ_61","contol_DJ_71","contol_DJ_81",
  "contol_DJ_91","contol_DJ_101")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyCon <- cronbach(cbind(dsPRE_POST_ok[,c(
  "contol_DJ_11","contol_DJ_21","contol_DJ_31","contol_DJ_41",
  "contol_DJ_51","contol_DJ_61","contol_DJ_71","contol_DJ_81",
  "contol_DJ_91","contol_DJ_101")]))
ALFAdejoyCon


## Compute average score for Risk Perception in after the driving test
dsPRE_POST_ok$Percezione_rischio2 <- rowMeans(dsPRE_POST_ok[,c("RP_scale_12",
                                                              "RP_scale_22",
                                                              "RP_scale_32",
                                                              "RP_scale_42",
                                                              "RP_scale_52")],
                                             na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFArisk2 <- cronbach(cbind(dsPRE_POST_ok[,c("RP_scale_12",
                                               "RP_scale_22",
                                               "RP_scale_32",
                                               "RP_scale_42",
                                               "RP_scale_52")]))
ALFArisk2
```

```{r}
## Compute average score for Vankov scale in after the driving test
dsPRE_POST_ok$Vankov2 <- rowMeans(dsPRE_POST_ok[,c("DUI_int2",
                                                  "Inst_att2",
                                                  "Per_cont2",
                                                  "Self_eff2",
                                                  "Des_norm2",
                                                  "Sub_norm2",
                                                  "Aff_att2")],
                                             na.rm=TRUE)
```

```{r}
## Cronbach's alfa
library(psy)
ALFAvankov2 <- cronbach(cbind(dsPRE_POST_ok[,c("DUI_int2",
                                                  "Inst_att2",
                                                  "Per_cont2",
                                                  "Self_eff2",
                                                  "Des_norm2",
                                                  "Sub_norm2",
                                                  "Aff_att2")]))
ALFAvankov2
```

```{r}
## Compute average score for Dejoy scale (involvement) in after the driving test
dsPRE_POST_ok$Dejoy_coinvolgimento2 <- rowMeans(dsPRE_POST_ok[,c(
  "involv_DJ_12","involv_DJ_22",
  "involv_DJ_32","involv_DJ_42",
  "involv_DJ_52","involv_DJ_62",
  "involv_DJ_72","involv_DJ_82",
  "involv_DJ_92","involv_DJ_102")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyC2 <- cronbach(cbind(dsPRE_POST_ok[,c(
  "involv_DJ_12","involv_DJ_22",
  "involv_DJ_32","involv_DJ_42",
  "involv_DJ_52","involv_DJ_62",
  "involv_DJ_72","involv_DJ_82",
  "involv_DJ_92","involv_DJ_102")]))
ALFAdejoyC2
```

```{r}
## Compute average score for Dejoy scale (gravity) in the before the driving test
dsPRE_POST_ok$Dejoy_gravità2 <- rowMeans(dsPRE_POST_ok[,c(
  "grav_DJ_12","grav_DJ_22","grav_DJ_32","grav_DJ_42",
  "grav_DJ_52","grav_DJ_62","grav_DJ_72","grav_DJ_82",
  "grav_DJ_92","grav_DJ_102")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyG2 <- cronbach(cbind(dsPRE_POST_ok[,c(
  "grav_DJ_12","grav_DJ_22","grav_DJ_32","grav_DJ_42",
  "grav_DJ_52","grav_DJ_62","grav_DJ_72","grav_DJ_82",
  "grav_DJ_92","grav_DJ_102")]))
ALFAdejoyG2
```

```{r}
## Compute average score for Dejoy scale (control) in the before the driving test
dsPRE_POST_ok$Dejoy_controllo2 <- rowMeans(dsPRE_POST_ok[,c(
  "contol_DJ_12","contol_DJ_22","contol_DJ_32","contol_DJ_42",
  "contol_DJ_52","contol_DJ_62","contol_DJ_72","contol_DJ_82",
  "contol_DJ_92","contol_DJ_102")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyCon2 <- cronbach(cbind(dsPRE_POST_ok[,c(
  "contol_DJ_12","contol_DJ_22","contol_DJ_32","contol_DJ_42",
  "contol_DJ_52","contol_DJ_62","contol_DJ_72","contol_DJ_82",
  "contol_DJ_92","contol_DJ_102")]))
ALFAdejoyCon2
```

## *Campione*

```{r, fig.width=8, fig.height=5, echo=TRUE}
dsPRE_POST_ok2 <- subset(dsPRE_POST_ok, Vankov2 != "NaN")
dsPRE_POST_ok2$Test <- factor(dsPRE_POST_ok2$Test)
library(table1)
T1 <- table1(~ Genere + Età + Istruzione + Lavoro +
         Residenza + Patente + Sanzioni_alcool +
         Sanzioni_neo_patentati + game1 + Fr_alcol1 +
         M_drink1 + Fam_guida1 | Test, data = dsPRE_POST_ok2)
T1
```

## *Misure di percezione del rischio*

```{r, echo=TRUE}
levels(dsPRE_POST_ok2$Test)[match("Guida da ubriaco",
                                  levels(dsPRE_POST_ok2$Test))] <- "Driving under the influence"
levels(dsPRE_POST_ok2$Test)[match("Guida da sobrio",
                                  levels(dsPRE_POST_ok2$Test))] <- "Driving sober"
T2 <- table1(~ Percezione_rischio + Percezione_rischio2 +
               Vankov + Vankov2 + 
               Dejoy_coinvolgimento + Dejoy_coinvolgimento2 +
               Dejoy_gravità + Dejoy_gravità2 +
               Dejoy_controllo + Dejoy_controllo2 | 
               Test, data = dsPRE_POST_ok2, overall = F)
T2
```

## *Illusione del controllo e bias ottimistico*

```{r, echo=TRUE}
dsPRE_POST_ok2$IllusionC1 <- rowMeans(dsPRE_POST_ok2[,c("ill_co11","ill_co12")],
                                             na.rm=TRUE)
dsPRE_POST_ok2$IllusionC2 <- rowMeans(dsPRE_POST_ok2[,c("ill_co21","ill_co22")],
                                             na.rm=TRUE)
dsPRE_POST_ok2$OptimismB1 <- rowMeans(dsPRE_POST_ok2[,c("Ob11","Ob21")],
                                             na.rm=TRUE)
dsPRE_POST_ok2$OptimismB2 <- rowMeans(dsPRE_POST_ok2[,c("Ob21","Ob22")],
                                             na.rm=TRUE)

dsPRE_POST_ok2$game_cat <- factor(dsPRE_POST_ok2$game1, levels = c(1,2,3,4,5),
                                  labels = c("Never","Very rarely","Only 1/2 times a week","3/4 times a week","Very often/every day"))
dsPRE_POST_ok2$FrAlcohol_cat <- factor(dsPRE_POST_ok2$Fr_alcol1, levels = c(1,2,3,4,5,6),
                                  labels = c("Never","Once","2/4 times","2/3 times a week","5/6 times a week","6 or more times a week"))
dsPRE_POST_ok2$NDrink_cat <- factor(dsPRE_POST_ok2$M_drink1, levels = c(1,2,3,4,5,6),
                                  labels = c("None","1/2 drinks","3/4 drinks","5/6 drinks","7/9 drinks","10 or more drinks"))
dsPRE_POST_ok2$FamGuida_cat <- factor(dsPRE_POST_ok2$Fam_guida1, levels = c(1,2,3,4),
                                  labels = c("Never","Sometimes","Often","Very often"))

dsPRE_POST_ok2$Passed <- factor(dsPRE_POST_ok2$Passed, levels = c(0,1),
                                labels = c("not passed","passed"))

library(lubridate)
dsPRE_POST_ok2$Duration2 <- ms(dsPRE_POST_ok2$Duration)
dsPRE_POST_ok2$Duration3 <- as.numeric(dsPRE_POST_ok2$Duration2)
dsPRE_POST_ok2$Durata <- dsPRE_POST_ok2$Duration3/60

T3 <- table1(~ Percezione_rischio + Percezione_rischio2 +
               Dejoy_coinvolgimento + Dejoy_coinvolgimento2 +
               Dejoy_gravità + Dejoy_gravità2 +
               IllusionC1 + IllusionC2 + 
               OptimismB1 + OptimismB2 + Ob31 + Ob32 +
               game_cat + FrAlcohol_cat + NDrink_cat + FamGuida_cat | 
               Test, data = dsPRE_POST_ok2, overall = F)
T3
```

# ***CORRELAZIONI***

## *Solo prima prova*

```{r, echo=TRUE}
library(psych)
pairs.panels(dsPRE_POST_ok2[, c("Percezione_rischio", 
                                "Dejoy_coinvolgimento","Dejoy_gravità",
                                "IllusionC1","OptimismB1","Ob31")],
             smooth = TRUE,      # If TRUE, draws loess smooths
             scale = FALSE,      # If TRUE, scales the correlation text font
             density = TRUE,     # If TRUE, adds density plots and histograms
             ellipses = FALSE,    # If TRUE, draws ellipses
             method = "pearson", # Correlation method (also "spearman" or "kendall")
             pch = 21,           # pch symbol
             lm = FALSE,         # If TRUE, plots linear fit rather than the LOESS (smoothed) fit
             cor = TRUE,         # If TRUE, reports correlations
             jiggle = TRUE,     # If TRUE, data points are jittered
             factor = 2,         # Jittering factor
             hist.col = 4,       # Histograms color
             stars = TRUE,       # If TRUE, adds significance level with stars
             ci = TRUE) 

```

## *Solo seconda prova - guida da sobri*

```{r, echo=TRUE}
Sober <- subset(dsPRE_POST_ok2,
                  Test == "Driving sober")

library(psych)
pairs.panels(Sober[, c("Percezione_rischio2", 
                                "Dejoy_coinvolgimento2","Dejoy_gravità2",
                                "IllusionC2","OptimismB2","Ob32",
                                "Durata","Wpt")],
             smooth = TRUE,      # If TRUE, draws loess smooths
             scale = FALSE,      # If TRUE, scales the correlation text font
             density = TRUE,     # If TRUE, adds density plots and histograms
             ellipses = FALSE,    # If TRUE, draws ellipses
             method = "pearson", # Correlation method (also "spearman" or "kendall")
             pch = 21,           # pch symbol
             lm = FALSE,         # If TRUE, plots linear fit rather than the LOESS (smoothed) fit
             cor = TRUE,         # If TRUE, reports correlations
             jiggle = TRUE,     # If TRUE, data points are jittered
             factor = 2,         # Jittering factor
             hist.col = 4,       # Histograms color
             stars = TRUE,       # If TRUE, adds significance level with stars
             ci = TRUE) 

```

## *Solo seconda prova - guida da ubriachi*

```{r, echo=TRUE}
DUI <- subset(dsPRE_POST_ok2,
                  Test == "Driving under the influence")

library(psych)
pairs.panels(DUI[, c("Percezione_rischio2", 
                                "Dejoy_coinvolgimento2","Dejoy_gravità2",
                                "IllusionC2","OptimismB2","Ob32",
                                "Durata","Wpt")],
             smooth = TRUE,      # If TRUE, draws loess smooths
             scale = FALSE,      # If TRUE, scales the correlation text font
             density = TRUE,     # If TRUE, adds density plots and histograms
             ellipses = FALSE,    # If TRUE, draws ellipses
             method = "pearson", # Correlation method (also "spearman" or "kendall")
             pch = 21,           # pch symbol
             lm = FALSE,         # If TRUE, plots linear fit rather than the LOESS (smoothed) fit
             cor = TRUE,         # If TRUE, reports correlations
             jiggle = TRUE,     # If TRUE, data points are jittered
             factor = 2,         # Jittering factor
             hist.col = 4,       # Histograms color
             stars = TRUE,       # If TRUE, adds significance level with stars
             ci = TRUE) 

```

# ***REGRESSIONI***

## *Percezione del rischio*

```{r, fig.width=8, fig.height=5, fig.align='center', echo=TRUE}
library(tidyr)
dsPRE_POST_ok2$Code <-  factor(dsPRE_POST_ok2$Code3)
ds_long <- reshape(dsPRE_POST_ok2, 
    varying=list(c("Percezione_rischio","Percezione_rischio2"),
                 c("Dejoy_coinvolgimento","Dejoy_coinvolgimento2"),
                 c("Dejoy_gravità","Dejoy_gravità2"),
                 c("IllusionC1","IllusionC2"),
                 c("Ob11","Ob21"),
                 c("Ob31","Ob32")),
    v.names=c("Percezione_rischio","Dejoy_coinvolgimento",
              "Dejoy_gravità","IllusionC",
              "OptimismB","PoliceStop"),
    direction="long", 
    idvar="Code3")

ds_long$Time <- factor(ds_long$time, levels = c(1,2),
                       labels = c("BASELINE","POST"))
```

### *Modello 1*

```{r, echo=TRUE}
library(lme4)
library(lmerTest)
library(car)

M1 <- lmer(Percezione_rischio ~ Test*Time + Durata + Genere +
             (1|Code), data = ds_long, REML = F)
Anova(M1, type = 3)
summary(M1)
```

#### Confidence intervals

```{r}
confint(M1, level=0.95)
```

#### Standardized coefficients

```{r}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M1)
```

```{r, fig.width=8,fig.height=5}
library(interactions)
library(jtools)
FY <- cat_plot(M1, pred = Test, modx = Time,
              interval = TRUE, plot.points = TRUE,
              point.size = 1, jitter = .5,
              y.label = "Risk perception\n", 
              x.label = "\nCondition",
              colors = c("red", "blue"),
              legend.main = "Session",
              line.thickness = .5) + theme_apa(legend.pos = "bottomleft",
                                               legend.use.title = TRUE,
                                               legend.font.size = 12,
                                               facet.title.size = 12) + 
  theme(axis.text.x = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.text.y = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.title = element_text(face = "bold")) +
  scale_y_continuous(limits=c(1,11), breaks=seq(1, 10,1))
FY
```

### *Modello 2: Gender interactions (3-way n.s.)*

```{r, echo=TRUE}
M2 <- lmer(Percezione_rischio ~ Test*Time + Test*Genere + Durata +
             Time*Genere + 
             (1|Code), data = ds_long, REML = F) # I checked the three-way and it is n.s.

Anova(M2, type = 3)
summary(M2)
```

#### Confidence intervals

```{r}
confint(M2, level=0.95)
```

#### Standardized coefficients

```{r}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M2)
```

#### Pairwise comparisons

```{r, echo=TRUE}
library(emmeans)
library(qualitema)
# x <- emm_cohend(
#   model = M2,
#   formula = pairwise ~ Time | Test,
#   adjust = "bonferroni")

```

```{r, fig.width=8,fig.height=5}
library(interactions)
library(jtools)
FY <- cat_plot(M2, pred = Test, modx = Time,
              interval = TRUE, plot.points = TRUE,
              point.size = 1, jitter = .5,
              y.label = "Risk perception\n", 
              x.label = "\nCondition",
              colors = c("red", "blue"),
              legend.main = "Session",
              line.thickness = .5) + theme_apa(legend.pos = "bottomleft",
                                               legend.use.title = TRUE,
                                               legend.font.size = 12,
                                               facet.title.size = 12) + 
  theme(axis.text.x = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.text.y = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.title = element_text(face = "bold")) +
  scale_y_continuous(limits=c(1,11), breaks=seq(1, 10,1))
FY
```

### *Modello 2b: Illusion of control*

```{r, echo=TRUE}
M3 <- lmer(Percezione_rischio ~ Test*Time + Test*Genere +
             Time*Genere + IllusionC + Durata +
             (1|Code), data = ds_long, REML = F) # I checked the three-way and it is n.s.

Anova(M3, type = 3)
summary(M3)
```

#### Confidence intervals

```{r}
confint(M3, level=0.95)
```

#### Standardized coefficients

```{r}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M3)
```

### *Modello 2c: Illusion of control + optmism*

```{r, echo=TRUE}
M4 <- lmer(Percezione_rischio ~ Test*Time + Test*Genere +
             Time*Genere + IllusionC + OptimismB + Durata +
             (1|Code), data = ds_long, REML = F) # I checked the three-way and it is n.s.

Anova(M4, type = 3)
summary(M4)
```

#### Confidence intervals

```{r}
confint(M4, level=0.95)
```

#### Standardized coefficients

```{r}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M4)
```

### *Modello 3: DEJOY SERIOUSNESS SCALE*

```{r, fig.width=8,fig.height=5, echo=TRUE}
M3 <- lmer(Percezione_rischio ~ Test*Time + Time*Dejoy_gravità +
             Test*Dejoy_gravità + Genere + Durata +
             (1|Code), data = ds_long, REML = F)
Anova(M3, type = 3)
summary(M3)
confint(M3)
stdCoef.merMod(M3)

# M3.3 <- lmer(Percezione_rischio ~ Test*Time*Dejoy_gravità +
#              Genere + Durata +
#              (1|Code), data = ds_long, REML = F)
# Anova(M3.3, type = 3)
# summary(M3.3)
# confint(M3.3)
# stdCoef.merMod(M3.3)
```

```{r, echo=TRUE}
library(interactions)
sim_slopes(M3, pred = Dejoy_gravità, modx = Time, confint = T)
```

```{r, fig.width=8,fig.height=5, echo=TRUE}
FY2 <- interact_plot(M3, pred = Dejoy_gravità, modx = Time,
              interval = TRUE, plot.points = TRUE,
              point.size = 1, jitter = .5,
              y.label = "Risk perception\n", 
              x.label = "\nDejoy seriousness",
              colors = c("red", "blue"),
              legend.main = "Session",
              line.thickness = 1) + theme_apa(legend.pos = "bottomleft",
                                               legend.use.title = TRUE,
                                               legend.font.size = 12,
                                               facet.title.size = 12) + 
  theme(axis.text.x = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.text.y = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.title = element_text(face = "bold")) +
  scale_y_continuous(limits=c(1,11), breaks=seq(1, 10,1))
FY2
```

```{r, echo=TRUE}
# write.table(dsPRE_POST_ok2, file = "ds_pre_post.csv", sep = ",",
#            na = " ", dec = ".", row.names = FALSE,
#            col.names = TRUE)
```

# ***FOLLOW UP***

```{r, echo=TRUE}
followup <- read.table(file = "followup.csv", sep = ",", dec = ".", header = TRUE)
followup$Code <- as.factor(followup$Cod_ID_min_b)

#left_join(dsPRE_POST_ok2, followup, by = "Code")
```

```{r, echo=TRUE, results='markup', fig.height=5, fig.width=8}
## Compute average score for Risk Perception in the followup
followup$Percezione_rischio_f <- rowMeans(followup[,c("RP_scale_1_b",
                                                              "RP_scale_2_b",
                                                              "RP_scale_3_b",
                                                              "RP_scale_4_b",
                                                              "RP_scale_5_b")],
                                             na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAriskb <- cronbach(cbind(followup[,c("RP_scale_1_b","RP_scale_2_b",
                                       "RP_scale_3_b","RP_scale_4_b",
                                       "RP_scale_5_b")]))
ALFAriskb

## Compute average score for Vankov scale in the followup
followup$Vankov_f <- rowMeans(followup[,c("DUI_int_b","Inst_att_b",
                                        "Per_cont_b","Self_eff_b",
                                        "Des_norm_b","Sub_norm_b",
                                        "Aff_att_b")],
                                             na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAvankovb <- cronbach(cbind(followup[,c("DUI_int_b","Inst_att_b",
                                        "Per_cont_b","Self_eff_b",
                                        "Des_norm_b","Sub_norm_b",
                                        "Aff_att_b")]))
ALFAvankovb

## Compute average score for Dejoy scale (involvement) in the followup
followup$Dejoy_coinvolgimento_f <- rowMeans(followup[,c(
  "involv_DJ_1_b","involv_DJ_2_b",
  "involv_DJ_3_b","involv_DJ_4_b",
  "involv_DJ_5_b","involv_DJ_6_b",
  "involv_DJ_7_b","involv_DJ_8_b",
  "involv_DJ_9_b","involv_DJ_10_b")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyCb <- cronbach(cbind(followup[,c(
  "involv_DJ_1_b","involv_DJ_2_b",
  "involv_DJ_3_b","involv_DJ_4_b",
  "involv_DJ_5_b","involv_DJ_6_b",
  "involv_DJ_7_b","involv_DJ_8_b",
  "involv_DJ_9_b","involv_DJ_10_b")]))
ALFAdejoyCb

## Compute average score for Dejoy scale (gravity) in the before the driving test
followup$Dejoy_gravità_f <- rowMeans(followup[,c(
  "grav_DJ_1_b","grav_DJ_2_b","grav_DJ_3_b","grav_DJ_4_b",
  "grav_DJ_5_b","grav_DJ_6_b","grav_DJ_7_b","grav_DJ_8_b",
  "grav_DJ_9_b","grav_DJ_10_b")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyGb <- cronbach(cbind(followup[,c(
  "grav_DJ_1_b","grav_DJ_2_b","grav_DJ_3_b","grav_DJ_4_b",
  "grav_DJ_5_b","grav_DJ_6_b","grav_DJ_7_b","grav_DJ_8_b",
  "grav_DJ_9_b","grav_DJ_10_b")]))
ALFAdejoyGb

## Compute average score for Dejoy scale (control) in the before the driving test
followup$Dejoy_controllo_f <- rowMeans(followup[,c(
  "contol_DJ_1_b","contol_DJ_2_b","contol_DJ_3_b","contol_DJ_4_b",
  "contol_DJ_5_b","contol_DJ_6_b","contol_DJ_7_b","contol_DJ_8_b",
  "contol_DJ_9_b","contol_DJ_10_b")],na.rm=TRUE)

## Cronbach's alfa
library(psy)
ALFAdejoyConb <- cronbach(cbind(followup[,c(
  "contol_DJ_1_b","contol_DJ_2_b","contol_DJ_3_b","contol_DJ_4_b",
  "contol_DJ_5_b","contol_DJ_6_b","contol_DJ_7_b","contol_DJ_8_b",
  "contol_DJ_9_b","contol_DJ_10_b")]))
ALFAdejoyConb

followup$Esperienza_sim <- rowMeans(followup[,c("esp_sim_1","esp_sim_2",
                                                "esp_sim_3")],na.rm=TRUE)

ALFAesp_sim <- cronbach(cbind(followup[,c(
  "esp_sim_1","esp_sim_2","esp_sim_3")]))
ALFAesp_sim


followup$Self_efficacy <- rowMeans(followup[,c("DUI_int_b","Inst_att_b",
                                                "Per_cont_b","Self_eff_b",
                                                "Des_norm_b","Sub_norm_b",
                                                "Aff_att_b")],na.rm=TRUE)

ALFAself_eff <- cronbach(cbind(followup[,c("DUI_int_b","Inst_att_b",
                                           "Sub_norm_b")]))
ALFAself_eff

followup$fr_drive <- car::recode(followup$fr_drive, "NA=0; 1=1; 5=2; 6=3; 7=4")
```

# ***STATISTICHE DESCRITTIVE***

```{r, echo=TRUE,fig.width=8, fig.height=5, echo=TRUE}
followup$Multe <- factor(followup$Multe, levels = c(1,2),
                         labels = c("Sì","No"))
followup$Patente <- factor(followup$license, levels = c(1,2),
                         labels = c("Sì","No"))
followup$Incidenti <- factor(followup$inc, levels = c(1,2),
                         labels = c("Sì","No"))
followup$DUI <- factor(followup$alc_drive, levels = c(1,2),
                         labels = c("Sì","No"))
library(table1)
T1f <- table1(~ fr_drive + Patente + Multe + Incidenti + 
               DUI + Esperienza_sim, data = followup)
T1f


followup$IllusionCf <- rowMeans(followup[,c("ill_co1_b","ill_co2_b")],
                                             na.rm=TRUE)
followup$OptimismBf <- rowMeans(followup[,c("Ob1_b","Ob2_b")],
                                             na.rm=TRUE)

ds_full2 <- left_join(followup, dsPRE_POST_ok2, by = "Code")
ds_full2 <- subset(ds_full2, Code != "va637ne")
ds_full2 <- subset(ds_full2, Test != "<NA>")
```

```{r, echo=TRUE}
ds_long_f <- reshape(ds_full2, 
    varying=list(c("Percezione_rischio","Percezione_rischio2","Percezione_rischio_f"),
                 c("Dejoy_coinvolgimento","Dejoy_coinvolgimento2","Dejoy_coinvolgimento_f"),
                 c("Dejoy_gravità","Dejoy_gravità2","Dejoy_gravità_f"),
                 c("IllusionC1","IllusionC2","IllusionCf"),
                 c("Ob11","Ob21","Ob1_b"),
                 c("Ob31","Ob32","Ob3_b")),
    v.names=c("Percezione_rischio","Dejoy_coinvolgimento",
              "Dejoy_gravità","IllusionC",
              "OptimismB","PoliceStop"),
    direction="long", 
    idvar="Code")

ds_long_f$Time <- factor(ds_long_f$time, levels = c(1,2,3),
                       labels = c("BASELINE","POST","FOLLOW_UP"))

#dsfp <- subset(ds_long_f, Time == "FOLLOW_UP")
T1p2 <- table1(~ Percezione_rischio + 
               Dejoy_coinvolgimento + 
               Dejoy_gravità +
               Dejoy_controllo +
               IllusionC + 
               OptimismB + PoliceStop | Test + Time, 
             data = ds_long_f)
T1p2

ds_long_f$game1 <- factor(ds_long_f$game1, levels = c(1,2,3,4,5),
                          labels = c("Never","Very rarely",
                                     "Only 1/2 times a week","3/4 times a week",
                                     "Very often/every day"))
ds_long_f$Fr_alcol1 <- factor(ds_long_f$Fr_alcol1, levels = c(1,2,3,4,5,6),
                          labels = c("Never","Once","2/4 times",
                                     "2/3 times a week","5/6 times a week",
                                     "6 or more times a week"))
ds_long_f$M_drink1 <- factor(ds_long_f$M_drink1, levels = c(1,2,3,4,5,6),
                          labels = c("None","1/2 drinks","3/4 drinks",
                                     "5/6 drinks","7/9 drinks",
                                     "10 or more drinks"))
ds_long_f$Fam_guida1 <- factor(ds_long_f$Fam_guida1, levels = c(1,2,3,4),
                          labels = c("Never","Sometimes","Often","Very often"))
```

# ***CORRELAZIONI***

## *Solo seconda prova - guida da sobri*

```{r, echo=TRUE, fig.width=9, fig.height=9}
library(psych)
ds_long_f$Duration <- as.numeric(ds_long_f$Duration)
post <- subset(ds_long_f, Time == "POST")
sober <- subset(post, Test == "Driving sober")
corPlot(sober[, c("Percezione_rischio", "Dejoy_coinvolgimento","Dejoy_gravità",
                    "IllusionC","OptimismB","PoliceStop",
                    "Duration1","Wpt")], stars = T,
        cex = .5, scale = F, diag = F, upper = F, show.legend = F) 
```

## *Solo seconda prova - guida da ubriachi*

```{r, echo=TRUE, fig.width=9, fig.height=9}
library(psych)
DUI <- subset(post, Test == "Driving under the influence")
corPlot(DUI[, c("Percezione_rischio", "Dejoy_coinvolgimento","Dejoy_gravità",
                    "IllusionC","OptimismB","PoliceStop",
                    "Duration1","Wpt")], stars = T,
        cex = .5, scale = F, diag = F, upper = F, show.legend = F)
```

## *Solo follow-up - guida da sobri*

```{r, echo=TRUE, fig.width=9, fig.height=9}
library(psych)
follow <- subset(ds_long_f, Time == "FOLLOW_UP")
f_sober <- subset(follow, Test == "Driving sober")
corPlot(f_sober[, c("Percezione_rischio", "Dejoy_coinvolgimento","Dejoy_gravità",
                    "IllusionC","OptimismB","PoliceStop")], stars = T,
        cex = .5, scale = F, diag = F, upper = F, show.legend = F)
```

## *Solo follow-up - guida da ubriachi*

```{r, echo=TRUE, fig.width=9, fig.height=9}
library(psych)
f_DUI <- subset(follow, Test == "Driving under the influence")
corPlot(f_DUI[, c("Percezione_rischio", "Dejoy_coinvolgimento","Dejoy_gravità",
                    "IllusionC","OptimismB","PoliceStop")], stars = T,
        cex = .5, scale = F, diag = F, upper = F, show.legend = F)
```

# ***REGRESSIONI (TUTTE LE SESSIONI)***

## *Modello 1*

```{r, echo=TRUE}
FUP<-c(-1,-1,2)
PRE_POST<-c(-1,1,0)
contrasts(ds_long_f$Time)<-cbind(FUP,PRE_POST)

M1fl <- lmer(Percezione_rischio ~ Test*Time + Durata + Genere +
             (1|Code), data = ds_long_f, REML = F)
Anova(M1fl, type = 3)
```

```{r, echo=TRUE}
summ(M1fl)
```

### Confidence intervals

```{r,echo=TRUE}
confint(M1fl)
```

### Standard coefficients

```{r, echo=TRUE}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M1fl)
```

## *Modello 2*

```{r, echo=TRUE}
M2fl <- lmer(Percezione_rischio ~ Test*Time + Test*Dejoy_gravità + Time*Dejoy_gravità + 
               Durata + Genere +
             (1|Code), data = ds_long_f, REML = F) # I checked the three-way and it is n.s.
Anova(M2fl, type = 3)
```

```{r, echo=TRUE}
summ(M2fl)
```

### Confidence intervals

```{r, echo=TRUE}
confint(M2fl)
```

### Standardized coefficients

```{r, echo=TRUE}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M2fl)
```

## *Modello 3*

```{r, echo=TRUE}
M3fl <- lmer(Percezione_rischio ~ Test*Time*Dejoy_gravità + 
               Durata + Genere +
             (1|Code), data = ds_long_f, REML = F) # I checked the three-way and it is n.s.
Anova(M3fl, type = 3)
```

```{r, echo=TRUE}
summ(M3fl)
```

### Confidence intervals

```{r, echo=TRUE}
confint(M3fl)
```

### Standardized coefficients

```{r, echo=TRUE}
stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
stdCoef.merMod(M3fl)
```

### Slope analysis: interaction Dejoy gravità x Session

Vedere ultimo modello...

## *Modello 4*

```{r, echo=TRUE}
M4fl <- lmer(Percezione_rischio ~ Test*Time*Dejoy_gravità + IllusionC + 
               + PoliceStop + OptimismBf  + Durata + Genere +
             (1|Code), data = ds_long_f, REML = F) # I checked the three-way and it is n.s.
Anova(M4fl, type = 3)
```

```{r, echo=TRUE}
summ(M4fl)
```

### Confidence intervals

```{r, echo=TRUE}
confint(M4fl)
```

### Standardized coefficients

```{r, echo=TRUE}
stdCoef.merMod(M4fl)
```

### Slope analysis: interaction Dejoy gravità x Session

```{r, echo=TRUE}
library(interactions)
sim_slopes(M4fl, pred = Dejoy_gravità, modx = Test, mod2 = Time, confint = T)
```

### Figure: interaction Condition x Session

```{r, echo=TRUE, fig.width=8,fig.height=5}
interact_plot(M4fl, pred = Dejoy_gravità, modx = Test, mod2 = Time,
              interval = TRUE, plot.points = F,
              point.size = 1, jitter = .5, 
              y.label = "Risk perception\n", 
              x.label = "\nDejoy seriousness",
              mod2.labels = c("BASELINE","POST-TEST","FOLLOW-UP"),
              colors = c("red", "blue","darkgreen"),
              legend.main = "Session",
              line.thickness = .5) + theme_apa(legend.pos = "bottomleft",
                                               legend.use.title = F,
                                               legend.font.size = 10,
                                               facet.title.size = 10) + 
  theme(axis.text.x = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.text.y = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.title = element_text(face = "bold")) +
  scale_y_continuous(limits=c(1,10), breaks=seq(1, 10,1)) +
  scale_x_continuous(limits=c(1,5), breaks=seq(1, 5,1))
```
# ***REGRESSIONE FOLLOWUP***

## *Modello 1*

```{r, echo=TRUE}
m1f <- lm(Percezione_rischio ~ Dejoy_gravità + Dejoy_coinvolgimento , follow)
library(jtools)
summ(m1f)
```

### Confidence intervals

```{r, echo=TRUE}
confint(m1f, level=0.95)
# How do I compute standardized effects for lm models?
```

### Standard coefficients

```{r, echo=TRUE}
m1fs <- lm(Percezione_rischio ~ scale(Dejoy_gravità) + scale(Dejoy_coinvolgimento), follow)
library(jtools)
summ(m1fs)
```

## *Modello 2*

```{r, echo=TRUE}
m2f <- lm(Percezione_rischio ~ 
            Dejoy_gravità*PoliceStop + esp_sim_1 + esp_sim_2 + esp_sim_3 +
            fr_drive + Genere, follow)
summ(m2f)
```

### Confidence intervals

```{r, echo=TRUE}
confint(m2f)
```

### Standard coefficients

```{r, echo=TRUE}
m2fs <- lm(Percezione_rischio ~ scale(Dejoy_gravità)*scale(PoliceStop) + 
            scale(esp_sim_1) + scale(esp_sim_2) + scale(esp_sim_3) + 
            scale(fr_drive) + Genere, follow)
summ(m2fs)
```

### Simple slopes: interaction Dejoy gravità x frequency of driving

```{r, echo=TRUE}
library(interactions)
sim_slopes(m2f, pred = Dejoy_gravità, modx = PoliceStop, confint = T)
```

### Figure: interaction Dejoy gravità x Effectiveness of simulator

```{r, echo=TRUE, fig.width=8,fig.height=5}
library(interactions)
library(jtools)
FUP1 <- interact_plot(m2f, pred = Dejoy_gravità, modx = PoliceStop,
              interval = TRUE, plot.points = F, 
              point.size = 1, jitter = .5,
              y.label = "Risk perception\n", 
              x.label = "\nDejoy seriousness",
              colors = c("red", "blue", "red"),
              modx.labels = c("Low likelihood","Average likelihood","High likelihood"),
              legend.main = "Perceived likelihood of being\nstopped by the police",
              line.thickness = 1) + theme_apa(legend.pos = "bottomright",
                                               legend.use.title = TRUE,
                                               legend.font.size = 9,
                                               facet.title.size = 8) + 
  theme(axis.text.x = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.text.y = element_text(color = "#212121", size = 10, face = "bold")) +
  theme(axis.title = element_text(face = "bold")) +
  scale_y_continuous(limits=c(1,10), breaks=seq(1, 10,1))
FUP1
```

-- ***THE END*** --