#install.packages(c("psych", "irr","pwr"))
data_dir <- "/Users/haisaosmanli/Desktop/CORE/SDA/DATA"
excel_files <- list.files(path = data_dir, pattern = "\\.xlsx$", full.names = TRUE)
excel_files
## [1] "/Users/haisaosmanli/Desktop/CORE/SDA/DATA/AD.xlsx"
## [2] "/Users/haisaosmanli/Desktop/CORE/SDA/DATA/AK.xlsx"
## [3] "/Users/haisaosmanli/Desktop/CORE/SDA/DATA/CD.xlsx"
## [4] "/Users/haisaosmanli/Desktop/CORE/SDA/DATA/EM.xlsx"
## [5] "/Users/haisaosmanli/Desktop/CORE/SDA/DATA/RL.xlsx"
## [6] "/Users/haisaosmanli/Desktop/CORE/SDA/DATA/VG.xlsx"
definitions <- read_excel(excel_files[1], sheet = 1)
## New names:
## • `` -> `...1`
## • `` -> `...3`
glimpse(definitions)
## Rows: 81
## Columns: 3
## $ ...1 <chr> …
## $ `Present = 1 point; absent=0 points: Place X in each box to indicate the presence of each indicator (total “X”s at the end in yellow box).` <chr> …
## $ ...3 <chr> …
read_researcher_data <- function(file_path) {
# extract researcher initials (eg. "AD.xlsx" → "AD")
researcher_id <- toupper(str_extract(basename(file_path), "^[A-Z]+"))
# get all sheet names
sheets <- excel_sheets(file_path)
# decide which sheets to read
if (researcher_id == "EM") {
target_sheets <- sheets[2:4] # EM → sheets 2,3 and 4
} else if (researcher_id == "VG") {
target_sheets <- sheets[2:4] # VG → sheets 2, 3, and 4
} else {
target_sheets <- sheets[2] # others → only sheet 2
}
# read and combine data from selected sheets
df_list <- lapply(target_sheets, function(sheet_name) {
df <- read_excel(file_path, sheet = sheet_name) %>%
clean_names() %>%
mutate(
researcher = researcher_id,
source_file = basename(file_path),
# clean sheet name (remove "(Pilot)" and convert to date)
date_collected = suppressWarnings(
ymd(str_trim(str_remove(sheet_name, "\\(.*\\)")))
)
)
return(df)
})
# combine all sheets vertically
bind_rows(df_list)
}
# read all files and merge vertically
metadata_all <- map_dfr(excel_files, read_researcher_data)
# preview the merged data
glimpse(metadata_all)
## Rows: 1,650
## Columns: 16
## $ researcher_initials <chr> "AD", "N/A", "N/A", "N/A", "N/A", "N/A", …
## $ weather_brief_description <chr> "Slightly Rainy, Cloudy", "N/A", "N/A", "…
## $ temperature_celsius <chr> "7", "N/A", "N/A", "N/A", "N/A", "N/A", "…
## $ number_of_photos_taken <chr> "4", "N/A", "N/A", "N/A", "N/A", "N/A", "…
## $ time_start_00_00am_pm <chr> "0.375", "N/A", "N/A", "N/A", "N/A", "N/A…
## $ time_end_00_00am_pm <chr> "0.41666666666666669", "N/A", "N/A", "N/A…
## $ location <chr> "Sunalta", "Sunalta", "Sunalta", "Sunalta…
## $ indicator_type <chr> "N/A", "Physical disorder", "Physical dis…
## $ code <chr> "N/A", "001", "002", "003a", "003b", "004…
## $ indicator_description <chr> "N/A", "Pile of discarded cigarette butts…
## $ present_yes_1 <chr> "N/A", "1", "1", "1", NA, "1", "1", NA, N…
## $ frequency_tally_total <chr> "N/A", "7", "2", "14", NA, "4", "1", NA, …
## $ field_notes_where_applicable <chr> "N/A", NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ researcher <chr> "AD", "AD", "AD", "AD", "AD", "AD", "AD",…
## $ source_file <chr> "AD.xlsx", "AD.xlsx", "AD.xlsx", "AD.xlsx…
## $ date_collected <date> 2025-10-19, 2025-10-19, 2025-10-19, 2025…
metadata_all <- metadata_all %>%
mutate(present_yes_1 = suppressWarnings(as.numeric(present_yes_1))) %>% # convert to numeric if possible
filter(!is.na(present_yes_1)) # keep only numeric (non-NA) values
# drop unwanted columns
metadata_all <- metadata_all %>%
select(
-matches("researcher_initials|weather|temperature|number_of_photos|time_start|time_end|field_notes|source_file")
)
# get rid of double sum
metadata_all <- metadata_all %>%
filter(!is.na(code) & code != "N/A")
metadata_all <- metadata_all %>%
#ensure frequency_tally_total is numeric
mutate(frequency_tally_total = suppressWarnings(as.numeric(frequency_tally_total))) %>%
# keep only rows with frequency > 0 and present_yes_1 = 1
filter(!is.na(frequency_tally_total) & frequency_tally_total > 0) %>%
filter(present_yes_1 != 0)
summary_totals <- metadata_all %>%
# clean up names and formats
mutate(
date_collected = as_date(date_collected),
indicator_type = str_to_lower(str_squish(indicator_type)) # make categories consistent
) %>%
# group by researcher, date, location, and indicator type
group_by(location, date_collected, researcher, indicator_type) %>%
summarise(
total_frequency = sum(frequency_tally_total, na.rm = TRUE),
.groups = "drop"
) %>%
# reshape to wide format so each category is a column
pivot_wider(
names_from = indicator_type,
values_from = total_frequency,
values_fill = 0
) %>%
arrange(location, date_collected, researcher)
# preview results
head(summary_totals)
## # A tibble: 6 × 6
## location date_collected researcher `physical disorder` positive social orde…¹
## <chr> <date> <chr> <dbl> <dbl>
## 1 Beltline… 2025-10-17 EM 94 1
## 2 Beltline… 2025-10-17 VG 120 27
## 3 Beltline… 2025-10-19 AD 117 14
## 4 Beltline… 2025-10-19 CD 124 20
## 5 Beltline… 2025-10-19 RL 127 28
## 6 Beltline… 2025-10-19 VG 136 27
## # ℹ abbreviated name: ¹`positive social order`
## # ℹ 1 more variable: `social disorder` <dbl>
metadata_all <- metadata_all %>%
mutate(location = if_else(is.na(location) | location == "", "Sunalta", location))
df <- summary_totals
df <- df %>%
rename(
physical_disorder = `physical disorder`,
positive_social_order = `positive social order`,
social_disorder = `social disorder`
)
summary_stats <- df %>%
group_by(location) %>%
summarise(across(c(physical_disorder, positive_social_order, social_disorder),
list(mean = mean, sd = sd, median = median), .names = "{col}_{fn}"))
summary_stats
## # A tibble: 3 × 10
## location physical_disorder_mean physical_disorder_sd physical_disorder_me…¹
## <chr> <dbl> <dbl> <dbl>
## 1 Beltline S… 121. 16.4 122
## 2 Sunalta 69.9 20.6 77
## 3 <NA> 0 0 0
## # ℹ abbreviated name: ¹physical_disorder_median
## # ℹ 6 more variables: positive_social_order_mean <dbl>,
## # positive_social_order_sd <dbl>, positive_social_order_median <dbl>,
## # social_disorder_mean <dbl>, social_disorder_sd <dbl>,
## # social_disorder_median <dbl>
# convert data to long format for plotting
df_long <- df %>%
pivot_longer(
cols = c(physical_disorder, positive_social_order, social_disorder),
names_to = "indicator",
values_to = "count"
)
# boxplot: counts per indicator, colored by researcher, faceted by location
ggplot(df_long, aes(x = researcher, y = count, fill = researcher)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.1, alpha = 0.5) + # show individual points
facet_wrap(~ indicator + location, scales = "free_y") +
labs(
title = "Distribution of Observations per Indicator and Researcher",
x = "Researcher",
y = "Count"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
)
df %>%
group_by(location, date_collected) %>%
summarise(across(c(physical_disorder, positive_social_order, social_disorder), sum)) %>%
pivot_longer(cols = c(physical_disorder, positive_social_order, social_disorder),
names_to = "indicator", values_to = "total") %>%
ggplot(aes(x = date_collected, y = total, fill = indicator)) +
geom_bar(stat = "identity", position = "dodge") +
facet_wrap(~location) +
theme_minimal() +
labs(title = "Daily Indicator Counts by Location")
## `summarise()` has grouped output by 'location'. You can override using the
## `.groups` argument.
df %>%
group_by(location, date_collected) %>%
summarise(across(c(physical_disorder, positive_social_order, social_disorder), sum)) %>%
pivot_longer(cols = c(physical_disorder, positive_social_order, social_disorder),
names_to = "indicator", values_to = "total") %>%
ggplot(aes(x = as.Date(date_collected, format="%m/%d/%y"), y = total, color = indicator, group = indicator)) +
geom_line() +
geom_point() +
facet_wrap(~location) +
theme_minimal() +
labs(title = "Trend of Observations Over Time")
## `summarise()` has grouped output by 'location'. You can override using the
## `.groups` argument.
# example: for 'physical disorder'
sunalta <- df %>% filter(location == "Sunalta") %>% pull(physical_disorder)
chumir <- df %>% filter(grepl("Beltline Sheldon Chumir", location)) %>% pull(physical_disorder)
effect_size <- (mean(sunalta) - mean(chumir)) / sd(c(sunalta, chumir))
pwr_result_pd <- pwr.t.test(d = effect_size, n = length(sunalta), sig.level = 0.05, type = "two.sample", alternative = "two.sided")
pwr_result_pd
##
## Two-sample t test power calculation
##
## n = 10
## d = 1.598664
## sig.level = 0.05
## power = 0.9219632
## alternative = two.sided
##
## NOTE: n is number in *each* group
# example: for 'physical disorder'
sunalta <- df %>% filter(location == "Sunalta") %>% pull(social_disorder)
chumir <- df %>% filter(grepl("Beltline Sheldon Chumir", location)) %>% pull(social_disorder)
effect_size <- (mean(sunalta) - mean(chumir)) / sd(c(sunalta, chumir))
pwr_result_sd <- pwr.t.test(d = effect_size, n = length(sunalta), sig.level = 0.05, type = "two.sample", alternative = "two.sided")
pwr_result_sd
##
## Two-sample t test power calculation
##
## n = 10
## d = 1.238313
## sig.level = 0.05
## power = 0.7450286
## alternative = two.sided
##
## NOTE: n is number in *each* group
# current effect size
d <- 1.227779
# calculate required n per group for power = 0.8
required_n_sd <- pwr.t.test(d = d, power = 0.8, sig.level = 0.05, type = "two.sample", alternative = "two.sided")
required_n_sd
##
## Two-sample t test power calculation
##
## n = 11.45818
## d = 1.227779
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
# example: for 'physical disorder'
sunalta <- df %>% filter(location == "Sunalta") %>% pull(positive_social_order)
chumir <- df %>% filter(grepl("Beltline Sheldon Chumir", location)) %>% pull(positive_social_order)
effect_size <- (mean(sunalta) - mean(chumir)) / sd(c(sunalta, chumir))
pwr_result_pso <- pwr.t.test(d = effect_size, n = length(sunalta), sig.level = 0.05, type = "two.sample", alternative = "two.sided")
pwr_result_pso
##
## Two-sample t test power calculation
##
## n = 10
## d = 1.066858
## sig.level = 0.05
## power = 0.6168342
## alternative = two.sided
##
## NOTE: n is number in *each* group
# current effect size
d <- 0.9069098
# calculate required n per group for power = 0.8
required_n_pso <- pwr.t.test(d = d, power = 0.8, sig.level = 0.05, type = "two.sample", alternative = "two.sided")
required_n_pso
##
## Two-sample t test power calculation
##
## n = 20.09283
## d = 0.9069098
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
icc_data <- df %>%
mutate(subject = paste0(ifelse(grepl("Sunalta", location), "Sunalta", "Chumir"), " ", date_collected)) %>%
select(subject, researcher, physical_disorder, positive_social_order, social_disorder)
# View result
icc_data
## # A tibble: 26 × 5
## subject researcher physical_disorder positive_social_order social_disorder
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 Chumir 20… EM 94 1 32
## 2 Chumir 20… VG 120 27 31
## 3 Chumir 20… AD 117 14 68
## 4 Chumir 20… CD 124 20 57
## 5 Chumir 20… RL 127 28 17
## 6 Chumir 20… VG 136 27 16
## 7 Chumir 20… AK 103 23 15
## 8 Chumir 20… EM 151 33 15
## 9 Chumir 20… EM 108 33 22
## 10 Chumir 20… VG 126 0 0
## # ℹ 16 more rows
# ICC(2,1) stands for Model 2, Type 1
# It measures the consistency or absolute agreement between raters when each target (in our case a
# location) is rated by the same set of raters, and those raters are
# considered a random sample from a larger population of possible raters
manual_icc2 <- function(df, measure){
df_var <- df %>%
dplyr::select(subject, researcher, all_of(measure)) %>%
dplyr::mutate(target = as.character(subject)) %>%
dplyr::select(target, researcher, all_of(measure))
df_var[[measure]] <- as.numeric(as.character(df_var[[measure]]))
# handle duplicates by averaging
df_wide <- df_var %>%
dplyr::group_by(target, researcher) %>%
dplyr::summarise(value = mean(.data[[measure]], na.rm = TRUE), .groups = "drop") %>%
tidyr::pivot_wider(names_from = researcher, values_from = value)
mat <- as.matrix(df_wide[,-1])
mode(mat) <- "numeric"
n_targets <- nrow(mat)
n_raters <- ncol(mat)
target_mean <- rowMeans(mat, na.rm = TRUE)
rater_mean <- colMeans(mat, na.rm = TRUE)
grand_mean <- mean(mat, na.rm = TRUE)
n_i <- rowSums(!is.na(mat))
ss_targets <- sum(n_i * (target_mean - grand_mean)^2)
n_j <- colSums(!is.na(mat))
ss_raters <- sum(n_j * (rater_mean - grand_mean)^2)
ss_residual <- 0
for(i in 1:n_targets){
for(j in 1:n_raters){
if(!is.na(mat[i,j])){
ss_residual <- ss_residual + (mat[i,j] - target_mean[i] - rater_mean[j] + grand_mean)^2
}
}
}
ms_targets <- ss_targets / (n_targets - 1)
ms_raters <- ss_raters / (n_raters - 1)
ms_residual <- ss_residual / ((n_targets - 1)*(n_raters - 1))
icc <- (ms_targets - ms_residual) /
(ms_targets + (n_raters - 1)*ms_residual + (n_raters/n_targets)*(ms_raters - ms_residual))
return(list(ICC = icc,
MS_target = ms_targets,
MS_rater = ms_raters,
MS_residual = ms_residual))
}
df <- icc_data
variables <- c("physical_disorder", "positive_social_order", "social_disorder")
df[variables] <- lapply(df[variables], function(x) as.numeric(as.character(x)))
icc_results <- lapply(variables, function(v){
manual_icc2(df, v)
})
names(icc_results) <- variables
icc_results
## $physical_disorder
## $physical_disorder$ICC
## EM
## 0.5483627
##
## $physical_disorder$MS_target
## [1] 1575.15
##
## $physical_disorder$MS_rater
## [1] 769.2225
##
## $physical_disorder$MS_residual
## EM
## 118.6104
##
##
## $positive_social_order
## $positive_social_order$ICC
## EM
## 0.5755499
##
## $positive_social_order$MS_target
## [1] 106.1357
##
## $positive_social_order$MS_rater
## [1] 24.66917
##
## $positive_social_order$MS_residual
## EM
## 9.9825
##
##
## $social_disorder
## $social_disorder$ICC
## EM
## 0.36733
##
## $social_disorder$MS_target
## [1] 422.3054
##
## $social_disorder$MS_rater
## [1] 454.0067
##
## $social_disorder$MS_residual
## EM
## 55.48369
ICC(2,1) Results Interpretation:
Inter-rater reliability was assessed using ICC(2,1), which evaluates
absolute agreement
between researchers for each observational measure.
Physical disorder showed moderate reliability (ICC = 0.55),
indicating a fair-to-good level
of agreement among raters.
Positive social order demonstrated slightly higher reliability (ICC
= 0.58), suggesting
consistent scoring across observers.
Social disorder had a lower ICC (0.37), reflecting weaker agreement
and suggesting that
this construct may require clearer rating criteria or additional
rater training.
Overall, reliability ranged from fair to good across domains, with
stronger consistency
observed for positive and physical order indicators.
Social disorder had a lower ICC (0.37), reflecting weaker agreement and suggesting that