R Notebook

library(readxl)
library(tidyverse)
library(psych)
library(ggplot2)
library(car)
library(stats)
library(viridis)

options(scipen = 999)
rm(list = ls())

library(readxl)
dataset_dfe <- read_excel("~/Desktop/QRMSEMINARS/vi. assignment2/dataAssignment2.xlsx")
View(dataset_dfe)

class(dataset_dfe)

[1] "tbl_df"     "tbl"        "data.frame"

str(dataset_dfe)

tibble [11,673 × 10] (S3: tbl_df/tbl/data.frame)
 $ laName        : chr [1:11673] "Barking and Dagenham" "Barking and Dagenham" "Barking and Dagenham" "Barking and Dagenham" ...
 $ timePeriod    : num [1:11673] 202425 202425 202425 202425 202425 ...
 $ regionName    : chr [1:11673] "London" "London" "London" "London" ...
 $ sex           : chr [1:11673] "Girls" "Boys" "Boys" "Boys" ...
 $ ethnicityMajor: chr [1:11673] "Mixed / Multiple ethnic groups" "Black / African / Caribbean / Black British" "White" "Asian / Asian British" ...
 $ ethnicityMinor: chr [1:11673] "White and Asian" "Caribbean" "Any other White background" "Any other Asian background" ...
 $ fsmStatus     : chr [1:11673] "FSM eligible" "Not known to be FSM eligible" "FSM eligible" "Not known to be FSM eligible" ...
 $ yearGroup     : chr [1:11673] "Year 1" "Year 1" "Year 1" "Year 1" ...
 $ percentTarget : num [1:11673] 100 85 72 76 82 79 NA 100 NA 68 ...
 $ percentFSM    : num [1:11673] 23.3 23.3 23.3 23.3 23.3 ...

dataset_dfe$laName <- as.factor(dataset_dfe$laName)
dataset_dfe$timePeriod <- as.factor(dataset_dfe$timePeriod)
dataset_dfe$regionName <- as.factor(dataset_dfe$regionName)
dataset_dfe$sex <- as.factor(dataset_dfe$sex)
dataset_dfe$ethnicityMajor <- as.factor(dataset_dfe$ethnicityMajor)
dataset_dfe$ethnicityMinor <- as.factor(dataset_dfe$ethnicityMinor)
dataset_dfe$fsmStatus <- as.factor(dataset_dfe$fsmStatus)
dataset_dfe$yearGroup <- as.factor(dataset_dfe$yearGroup)

str(dataset_dfe)

tibble [11,673 × 10] (S3: tbl_df/tbl/data.frame)
 $ laName        : Factor w/ 153 levels "Barking and Dagenham",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ timePeriod    : Factor w/ 1 level "202425": 1 1 1 1 1 1 1 1 1 1 ...
 $ regionName    : Factor w/ 9 levels "East Midlands",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ sex           : Factor w/ 2 levels "Boys","Girls": 2 1 1 1 2 1 1 2 2 1 ...
 $ ethnicityMajor: Factor w/ 5 levels "Asian / Asian British",..: 3 2 5 1 2 1 1 3 5 2 ...
 $ ethnicityMinor: Factor w/ 22 levels "African","All Asian / Asian British",..: 20 12 10 6 1 2 13 20 18 1 ...
 $ fsmStatus     : Factor w/ 2 levels "FSM eligible",..: 1 2 1 2 1 2 1 2 1 1 ...
 $ yearGroup     : Factor w/ 1 level "Year 1": 1 1 1 1 1 1 1 1 1 1 ...
 $ percentTarget : num [1:11673] 100 85 72 76 82 79 NA 100 NA 68 ...
 $ percentFSM    : num [1:11673] 23.3 23.3 23.3 23.3 23.3 ...

#Part 1 - Task 1

data_dfe <- dataset_dfe %>%
  filter(!is.na(percentTarget))

overall_stats <- data_dfe %>%
  summarise(
    Median = median(percentTarget),
    Mean = mean(percentTarget),
    SD = sd(percentTarget),
  )

print(overall_stats)

library(ggplot2)

mean_val <- mean(data_dfe$percentTarget, na.rm = TRUE)
median_val <- median(data_dfe$percentTarget, na.rm = TRUE)

ggplot(data_dfe, aes(x = percentTarget)) +
  geom_density(fill = "#56B4E9", alpha = 0.8, color = "black") +
  geom_vline(aes(xintercept = mean_val), color = "red", linetype = "dashed", size = 1) +
  geom_vline(aes(xintercept = median_val), color = "black", linetype = "solid", size = 1) +
  annotate("text", x = mean_val - 5, y = 0.01, 
           label = paste("Mean:", round(mean_val, 1)), color = "red", angle = 0, vjust = 1) +
  annotate("text", x = median_val + 5, y = 0.01, 
           label = paste("Median:", round(median_val, 1)), color = "black", angle = 0, vjust = -0.5) +
  
 theme_minimal() +
  labs(
    title = "Distribution of Phonics Attainment Scores in England",
    x = "percentTarget",
    y = "Density"
  ) 

#Part 1 - Task 2

sex_stats <- data_dfe %>%
  group_by(sex) %>%
  summarise(
    Median = median(percentTarget),
    Mean = mean(percentTarget),
    SD = sd(percentTarget)
  )

print(sex_stats)

boxplot_sex <- ggplot(data_dfe, aes(x = sex, y = percentTarget, fill = sex)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red", outlier.shape = 1) +

  theme_minimal() +
  labs(
    title = "Phonics Attainment by Sex",
    x = "Sex",
    y = "percentTarget"
  )

print(boxplot_sex)

library(tidyverse)
library(car)

#DV needs to be numeric
is.numeric(data_dfe$percentTarget)

[1] TRUE

boys_dfe <- data_dfe$percentTarget[data_dfe$sex == "Boys"]

print("Shapiro-Wilk Test: Boys")

[1] "Shapiro-Wilk Test: Boys"

if (length(boys_dfe) > 5000) {
  shapiro.test(sample(boys_dfe, 5000))
} else {
  shapiro.test(boys_dfe)
}

    Shapiro-Wilk normality test

data:  boys_dfe
W = 0.9608, p-value < 0.00000000000000022

girls_dfe <- data_dfe$percentTarget[data_dfe$sex == "Girls"]

print("Shapiro-Wilk Text: Girls")

[1] "Shapiro-Wilk Text: Girls"

if (length(girls_dfe) > 5000) {
  shapiro.test(sample(girls_dfe, 5000))
} else {
  shapiro.test(girls_dfe)
}

    Shapiro-Wilk normality test

data:  girls_dfe
W = 0.93387, p-value < 0.00000000000000022

leveneTest(percentTarget ~ sex, data = data_dfe)

Levene's Test for Homogeneity of Variance (center = median)
        Df F value                Pr(>F)    
group    1  110.63 < 0.00000000000000022 ***
      8243                                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

wilcox_result_sex <- wilcox.test(data_dfe$percentTarget ~ data_dfe$sex, alternative = "two.sided")
print(wilcox_result_sex)

    Wilcoxon rank sum test with continuity correction

data:  data_dfe$percentTarget by data_dfe$sex
W = 5105979, p-value < 0.00000000000000022
alternative hypothesis: true location shift is not equal to 0

#Part 1 Task 3

fsm_stats <- data_dfe %>%
  group_by(fsmStatus) %>%
  summarise(
    Median = median(percentTarget),
    Mean = mean(percentTarget),
    SD = sd(percentTarget)
  )

print(fsm_stats)

boxplot_fsm <- ggplot(data_dfe, aes(x = fsmStatus, y = percentTarget, fill = fsmStatus)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red", outlier.shape = 1) +

  theme_minimal() +
  labs(
    title = "Phonics Attainment by FSM Status",
    x = "FSM Status",
    y = "percentTarget"
  )

print(boxplot_fsm)

fsm_eligible <- data_dfe$percentTarget[data_dfe$fsmStatus == "FSM eligible"]

print("Shapiro-Wilk Test: FSM eligible")

[1] "Shapiro-Wilk Test: FSM eligible"

if (length(fsm_eligible) > 5000) {
  shapiro.test(sample(fsm_eligible, 5000))
} else {
  shapiro.test(fsm_eligible)
}

    Shapiro-Wilk normality test

data:  fsm_eligible
W = 0.98263, p-value < 0.00000000000000022

fsm_notKnown <- data_dfe$percentTarget[data_dfe$fsmStatus == "Not known to be FSM eligible"]

print("Shapiro-Wilk Text: Not known to be FSM eligible")

[1] "Shapiro-Wilk Text: Not known to be FSM eligible"

if (length(fsm_notKnown) > 5000) {
  shapiro.test(sample(fsm_notKnown, 5000))
} else {
  shapiro.test(fsm_notKnown)
}

    Shapiro-Wilk normality test

data:  fsm_notKnown
W = 0.93218, p-value < 0.00000000000000022

leveneTest(percentTarget ~ fsmStatus, data = data_dfe)

Levene's Test for Homogeneity of Variance (center = median)
        Df F value                Pr(>F)    
group    1  413.11 < 0.00000000000000022 ***
      8243                                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

wilcox_result_fsm <- wilcox.test(data_dfe$percentTarget ~ data_dfe$fsmStatus, alternative = "two.sided")
print(wilcox_result_fsm)

    Wilcoxon rank sum test with continuity correction

data:  data_dfe$percentTarget by data_dfe$fsmStatus
W = 3500711, p-value < 0.00000000000000022
alternative hypothesis: true location shift is not equal to 0

#Part 1 - Task 4

region_stats <- data_dfe %>%
  group_by(regionName) %>%
  summarise(
    Median = median(percentTarget),
    Mean = mean(percentTarget),
    SD = sd(percentTarget),
  ) %>%
arrange(desc(Median))
print(region_stats)

scatterplot <- ggplot(data = data_dfe, 
                      aes(x = percentFSM, y = percentTarget)) +
  geom_point(color = "#56B4E9") + 
  geom_smooth(method = lm, color = "red") +    
  ylim(0, 100) +
  theme_minimal() +
  labs(
    x = "percentFSM",
    y = "percentTarget"
  )

print(scatterplot)

model <- lm(percentTarget ~ percentFSM, data = data_dfe)
summary(model)

Call:
lm(formula = percentTarget ~ percentFSM, data = data_dfe)

Residuals:
    Min      1Q  Median      3Q     Max 
-77.685  -7.696   2.353   9.647  23.626 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 75.75804    0.49200 153.978 < 0.0000000000000002 ***
percentFSM   0.08403    0.02187   3.841             0.000123 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.55 on 8243 degrees of freedom
Multiple R-squared:  0.001787,  Adjusted R-squared:  0.001666 
F-statistic: 14.76 on 1 and 8243 DF,  p-value: 0.0001232

plot(model, which = 2)

#Part 2 - Task 1 

# the three Local Authorities are: Richmond Upon Thames, West Sussex, and Manchester 

la_summary <- data_dfe %>%
  group_by(laName) %>%
  summarise(
    Median_FSM = median(percentFSM, na.rm = TRUE )
    ) %>%
  arrange(Median_FSM)

print(la_summary)

library(dplyr)

fsm_rankings <- data_dfe %>%
  group_by(laName) %>%
  summarise(
    Median_FSM = median(percentFSM, na.rm = TRUE)) %>%
  arrange(desc(Median_FSM)) 

print(head(fsm_rankings, 5))

threeLas_comparison <- data_dfe %>%
  filter(laName %in% c("Richmond upon Thames", "West Sussex", "Manchester")) %>%
  group_by(laName) %>%
  summarise(
    Median_FSM = median(percentFSM, na.rm = TRUE),      
    Median_Score = median(percentTarget, na.rm = TRUE)
  )

print(threeLas_comparison)

#Part 2 - Task 2

library(tidyverse)
library(car)

three_las <- c("Richmond upon Thames", "West Sussex", "Manchester")
data_three_las <- data_dfe %>% 
  filter(laName %in% three_las)

data_three_las %>%
  group_by(laName) %>%
  summarise(p_value = shapiro.test(percentTarget)$p.value) %>%
  print()

leveneTest(percentTarget ~ laName, data = data_three_las)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value   Pr(>F)   
group   2  4.9269 0.008167 **
      197                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

kruskal_result <- kruskal.test(percentTarget ~ laName, data = data_three_las)
print(kruskal_result)

    Kruskal-Wallis rank sum test

data:  percentTarget by laName
Kruskal-Wallis chi-squared = 44.369, df = 2, p-value = 0.0000000002319

library(FSA)

dunn_result <- dunnTest(percentTarget ~ laName, 
                        data = data_three_las, 
                        method = "bonferroni")

print(dunn_result)

#Part 2 - 3 

library(dplyr)

data_three_las %>%
  group_by(laName) %>%
  summarise(
    Median_Score = median(percentTarget, na.rm = TRUE), 
    Mean_Score = mean(percentTarget, na.rm = TRUE)
  ) %>%
  arrange(Median_Score) 

NA

ggplot(data_three_las, aes(x = reorder (laName, percentTarget, median),
                           y = percentTarget,
                           fill = laName)) + 
         geom_boxplot(alpha = 0.7, show.legend = FALSE) + 
         theme_minimal() +
         labs(title = "Comparison of Three LAs",
              x = "Local Authority",
              y = "Phonics Attainment (%)") 

#part 2 - 4

data_WS <- data_dfe %>% 
  filter(laName == "West Sussex")

sex_stats_ws <- data_WS %>% 
  group_by(sex) %>% 
  summarise(Median = median(percentTarget, na.rm=TRUE), Count = n())
print(sex_stats_ws)

data_WS %>%
  group_by(sex) %>%
  summarise(p_value = shapiro.test(percentTarget)$p.value) %>%
  print()

leveneTest(percentTarget ~ sex, data = data_WS)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.8376 0.3633
      69               

result_sex <- t.test(percentTarget ~ sex, data = data_WS, var.equal = TRUE)
print(result_sex)

    Two Sample t-test

data:  percentTarget by sex
t = -2.136, df = 69, p-value = 0.03623
alternative hypothesis: true difference in means between group Boys and group Girls is not equal to 0
95 percent confidence interval:
 -13.5009383  -0.4609665
sample estimates:
 mean in group Boys mean in group Girls 
           66.68571            73.66667 

fsm_stats_ws <- data_WS %>% 
  group_by(fsmStatus) %>% 
  summarise(Median = median(percentTarget, na.rm=TRUE), Count = n())
print(fsm_stats_ws)

data_WS %>%
  group_by(fsmStatus) %>%
  summarise(p_value = shapiro.test(percentTarget)$p.value) %>%
  print()

leveneTest(percentTarget ~ fsmStatus, data = data_WS)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.7892 0.3774
      69               

print(wilcox.test(percentTarget ~ fsmStatus, data = data_WS))

    Wilcoxon rank sum test with continuity correction

data:  percentTarget by fsmStatus
W = 209, p-value = 0.000002897
alternative hypothesis: true location shift is not equal to 0

ethnicity_stats_ws <- data_WS %>% 
  group_by(ethnicityMajor) %>% 
  summarise(Median = median(percentTarget, na.rm=TRUE), Count = n()) %>% 
  arrange(desc(Median))
print(ethnicity_stats_ws)

data_WS %>%
  group_by(ethnicityMajor) %>%
  summarise(p_value = shapiro.test(percentTarget)$p.value) %>%
  print()

library(car)
leveneTest(percentTarget ~ ethnicityMajor, data = data_WS)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  4   1.749 0.1498
      66               

print(kruskal.test(percentTarget ~ ethnicityMajor, data = data_WS))

    Kruskal-Wallis rank sum test

data:  percentTarget by ethnicityMajor
Kruskal-Wallis chi-squared = 7.1045, df = 4, p-value = 0.1305