Linear Mixed-Effect Model

#Shapiro to check noramlity
mydata <- read.csv("/Users/peihang/Desktop/Risk_Level.csv")
density_data <- subset(mydata, Change_Type == "Density")$Risk_Perception_Level
speed_data <- subset(mydata, Change_Type == "Speed")$Risk_Perception_Level
shapiro_density <- shapiro.test(density_data)
shapiro_speed <- shapiro.test(speed_data)
print("Shapiro-Wilk Test for Density Group:")
## [1] "Shapiro-Wilk Test for Density Group:"
print(shapiro_density)
## 
##  Shapiro-Wilk normality test
## 
## data:  density_data
## W = 0.91488, p-value = 0.0004799
print("Shapiro-Wilk Test for Speed Group:")
## [1] "Shapiro-Wilk Test for Speed Group:"
print(shapiro_speed)
## 
##  Shapiro-Wilk normality test
## 
## data:  speed_data
## W = 0.94737, p-value = 0.01173
##Data does not follow a normal distribution.

#LMM
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
str(mydata)
## 'data.frame':    120 obs. of  5 variables:
##  $ Subject_ID           : chr  "S01" "S01" "S01" "S01" ...
##  $ Change_Type          : chr  "Density" "Density" "Density" "Density" ...
##  $ Direction            : chr  "Low_High" "Low_High" "Low_High" "High_Low" ...
##  $ Probe_Time           : chr  "Pre" "Transition" "Post" "Pre" ...
##  $ Risk_Perception_Level: int  3 7 6 4 2 5 7 8 4 3 ...
mydata$Subject_ID <- as.factor(mydata$Subject_ID)
mydata$Change_Type <- as.factor(mydata$Change_Type)
mydata$Direction <- as.factor(mydata$Direction)
mydata$Probe_Time <- as.factor(mydata$Probe_Time)

str(mydata)
## 'data.frame':    120 obs. of  5 variables:
##  $ Subject_ID           : Factor w/ 10 levels "S01","S02","S03",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Change_Type          : Factor w/ 2 levels "Density","Speed": 1 1 1 1 1 1 2 2 2 2 ...
##  $ Direction            : Factor w/ 2 levels "High_Low","Low_High": 2 2 2 1 1 1 2 2 2 1 ...
##  $ Probe_Time           : Factor w/ 3 levels "Post","Pre","Transition": 2 3 1 2 3 1 2 3 1 2 ...
##  $ Risk_Perception_Level: int  3 7 6 4 2 5 7 8 4 3 ...
#Use transition as the base time to compare the difference with pre and post
mydata$Probe_Time <- relevel(mydata$Probe_Time, ref = "Transition")

density_data <- subset(mydata, Change_Type == "Density")
speed_data <- subset(mydata, Change_Type == "Speed")
#Density high to low & low to high, subject ID as random effect
density_low_high <- subset(density_data, Direction == "Low_High")
density_high_low <- subset(density_data, Direction == "High_Low")

model_density_low_high <- lmer(Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID), data = density_low_high)
summary(model_density_low_high)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID)
##    Data: density_low_high
## 
## REML criterion at convergence: 108.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17852 -0.54172  0.01747  0.56793  1.73583 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Subject_ID (Intercept) 1.919    1.385   
##  Residual               1.504    1.226   
## Number of obs: 30, groups:  Subject_ID, 10
## 
## Fixed effects:
##                Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)      4.4000     0.5850 16.5791   7.521 9.79e-07 ***
## Probe_TimePost  -1.2000     0.5484 18.0000  -2.188   0.0421 *  
## Probe_TimePre   -3.0000     0.5484 18.0000  -5.470 3.39e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prb_TmPs
## Probe_TmPst -0.469         
## Probe_TimPr -0.469  0.500
##Finding 1:
##Density from low to high, 
##From Pre to Transition, Perceived Risk Level drops significantly (-3.0, p = 3.39e-05)
##From Transition to Post, Perceived Risk Level slightly increases (+1.2, p = 0.0421)

model_density_high_low <- lmer(Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID), data = density_high_low)
summary(model_density_high_low)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID)
##    Data: density_high_low
## 
## REML criterion at convergence: 104.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0996 -0.6753 -0.1121  0.6135  2.0521 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Subject_ID (Intercept) 1.037    1.018   
##  Residual               1.515    1.231   
## Number of obs: 30, groups:  Subject_ID, 10
## 
## Fixed effects:
##                Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)      1.5000     0.5052 20.2962   2.969 0.007499 ** 
## Probe_TimePost   0.1000     0.5504 18.0000   0.182 0.857865    
## Probe_TimePre    2.5000     0.5504 18.0000   4.542 0.000253 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prb_TmPs
## Probe_TmPst -0.545         
## Probe_TimPr -0.545  0.500
##Finding 2:
##Density from high to low,
##From Pre to Transition, Perceived Risk Level significantly decreases (-2.5, p = 0.000253)
##From Transition to Post, Perceived Risk Level remains almost unchanged (+0.1, p = 0.858)

Density Change Visualization (Barplot with Error Bars)

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
density_data <- density_data %>%
  mutate(Probe_Time = factor(Probe_Time, levels = c("Pre", "Transition", "Post")))

summary_data <- density_data %>%
  mutate(Probe_Time = factor(Probe_Time, levels = c("Pre", "Transition", "Post"))) %>%
  group_by(Probe_Time, Direction) %>%
  summarise(
    mean_risk = mean(Risk_Perception_Level),
    se_risk = sd(Risk_Perception_Level) / sqrt(n())
  )
## `summarise()` has grouped output by 'Probe_Time'. You can override using the
## `.groups` argument.
ggplot(summary_data, aes(x = Probe_Time, y = mean_risk, fill = Direction)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.6) +
  geom_errorbar(aes(ymin = mean_risk - se_risk, ymax = mean_risk + se_risk), 
                width = 0.2, position = position_dodge(0.6)) +
  theme_minimal() +
  labs(title = "Perceived Risk Level by Probe Time (Density Change)",
       x = "Probe Time",
       y = "Mean Perceived Risk Level") +
  scale_fill_manual(values = c("Low_High" = "blue", "High_Low" = "red"))

Speed

mydata$Probe_Time <- relevel(mydata$Probe_Time, ref = "Transition")

speed_low_high <- subset(mydata, Change_Type == "Speed" & Direction == "Low_High")
speed_high_low <- subset(mydata, Change_Type == "Speed" & Direction == "High_Low")

model_speed_low_high <- lmer(Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID), data = speed_low_high)
summary(model_speed_low_high)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID)
##    Data: speed_low_high
## 
## REML criterion at convergence: 112.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.69246 -0.58921 -0.01409  0.47345  1.61461 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Subject_ID (Intercept) 3.315    1.821   
##  Residual               1.463    1.210   
## Number of obs: 30, groups:  Subject_ID, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     3.500e+00  6.912e-01  1.376e+01   5.064 0.000183 ***
## Probe_TimePost -1.500e+00  5.409e-01  1.800e+01  -2.773 0.012539 *  
## Probe_TimePre   1.776e-16  5.409e-01  1.800e+01   0.000 1.000000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prb_TmPs
## Probe_TmPst -0.391         
## Probe_TimPr -0.391  0.500
##Finding 3:
##From Pre to Transition, there is no significant change in perceived risk level
##From Transition to Post, percieved risk level significantly decreases (-1.5, p = 0.0125)

model_speed_high_low <- lmer(Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID), data = speed_high_low)
summary(model_speed_high_low)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Risk_Perception_Level ~ Probe_Time + (1 | Subject_ID)
##    Data: speed_high_low
## 
## REML criterion at convergence: 113.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8122 -0.3817  0.1464  0.5665  1.6511 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Subject_ID (Intercept) 1.433    1.197   
##  Residual               2.130    1.459   
## Number of obs: 30, groups:  Subject_ID, 10
## 
## Fixed effects:
##                Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)      5.1000     0.5969 20.3978   8.544 3.56e-08 ***
## Probe_TimePost  -2.0000     0.6526 18.0000  -3.065  0.00668 ** 
## Probe_TimePre   -2.5000     0.6526 18.0000  -3.831  0.00122 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prb_TmPs
## Probe_TmPst -0.547         
## Probe_TimPr -0.547  0.500
##Finding 4:
##From Pre to Transition, perceived risk level significantly decreases, indicating that participants feel less risk when transitioning to a lower speed (-2.5, p = 0.00122)
##From Transition to Post, perceived risk level continues to decrease, suggesting that participants further adapt to the lower speed over time. (-2.0, p = 0.00668)

Speed Change Visualization (Barplot with Error Bars)

library(ggplot2)
library(dplyr)

speed_data <- subset(mydata, Change_Type == "Speed") %>%
  mutate(Probe_Time = factor(Probe_Time, levels = c("Pre", "Transition", "Post")))


summary_speed <- speed_data %>%
  group_by(Probe_Time, Direction) %>%
  summarise(
    mean_risk = mean(Risk_Perception_Level),
    se_risk = sd(Risk_Perception_Level) / sqrt(n())
  )
## `summarise()` has grouped output by 'Probe_Time'. You can override using the
## `.groups` argument.
ggplot(summary_speed, aes(x = Probe_Time, y = mean_risk, fill = Direction)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.6) +  
  geom_errorbar(aes(ymin = mean_risk - se_risk, ymax = mean_risk + se_risk),
                width = 0.2, position = position_dodge(0.6)) +  
  theme_minimal() +
  labs(title = "Perceived Risk Level by Probe Time (Speed Change)",
       x = "Probe Time",
       y = "Mean Perceived Risk Level",
       fill = "Direction") +  
  scale_fill_manual(values = c("Low_High" = "blue", "High_Low" = "red"))