load student survey
library(haven)
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
student <- read_sav(
"/Users/xutianyu/Downloads/PISA2025/CY09_MS_STU_PUF.sav",
encoding = "latin1"
)
rename the variables
agency <- agency %>%
rename(
country = CNT,
math_score = PV1MATH,
science_score = PV1SCIE,
plan_study = ST432Q02DA,
identify_gaps = ST432Q03DA,
self_question = ST432Q04DA,
adapt_strategy = ST432Q06DA,
ask_classmates = ST433Q01DA,
ask_teachers = ST433Q02DA,
ask_general = ST433Q04DA,
ask_after_trying = ST433Q06DA
)
ncol(agency)
## [1] 12
find NA proportion
library(tibble)
agency_na_prop <- agency %>%
group_by(country) %>%
summarise(across(where(is.numeric), ~ mean(is.na(.x))))
#View(agency_na_prop)
use NA proportion = 0.5 as a threshold, make a table for all the
countries
high_na_vars_per_country <- agency_na_prop %>%
mutate(across(-country, ~ .x > 0.5)) %>%
mutate(Variables_Over_50pct_NA = rowSums(select(., -country))) %>%
select(country, Variables_Over_50pct_NA)
#View(high_na_vars_per_country)
print(as_tibble(high_na_vars_per_country), n = Inf)
## # A tibble: 90 × 2
## country Variables_Over_50pct_NA
## <chr+lbl> <dbl>
## 1 ALB [Albania] 0
## 2 ARE [United Arab Emirates] 0
## 3 ARG [Argentina] 0
## 4 ARM [Armenia] 8
## 5 AUS [Australia] 0
## 6 AUT [Austria] 0
## 7 AZE [Azerbaijan] 0
## 8 BEL [Belgium] 0
## 9 BGR [Bulgaria] 0
## 10 BRA [Brazil] 0
## 11 BRN [Brunei Darussalam] 0
## 12 CAN [Canada] 0
## 13 CHE [Switzerland] 0
## 14 CHL [Chile] 0
## 15 COL [Colombia] 0
## 16 CRI [Costa Rica] 0
## 17 CZE [Czechia] 0
## 18 DEU [Germany] 0
## 19 DNK [Denmark] 0
## 20 DOM [Dominican Republic] 0
## 21 ECU [Ecuador] 0
## 22 ESP [Spain] 0
## 23 EST [Estonia] 0
## 24 FIN [Finland] 0
## 25 FRA [France] 0
## 26 GBR [United Kingdom] 0
## 27 GEO [Georgia] 0
## 28 GRC [Greece] 0
## 29 GTM [Guatemala] 0
## 30 HKG [Hong Kong (China)] 0
## 31 HRV [Croatia] 0
## 32 HUN [Hungary] 0
## 33 IDN [Indonesia] 0
## 34 IRL [Ireland] 0
## 35 ISL [Iceland] 0
## 36 ISR [Israel] 0
## 37 ITA [Italy] 0
## 38 JOR [Jordan] 0
## 39 JPN [Japan] 0
## 40 KAZ [Kazakhstan] 0
## 41 KEN [Kenya] 0
## 42 KGZ [Kyrgyzstan] 0
## 43 KHM [Cambodia] 8
## 44 KOR [Korea] 0
## 45 KSV [Kosovo] 0
## 46 LBN [Lebanon] 8
## 47 LTU [Lithuania] 0
## 48 LUX [Luxembourg] 0
## 49 LVA [Latvia] 0
## 50 MAC [Macao (China)] 0
## 51 MAR [Morocco] 0
## 52 MDA [Moldova] 0
## 53 MEX [Mexico] 0
## 54 MKD [North Macedonia] 0
## 55 MLT [Malta] 0
## 56 MNE [Montenegro] 0
## 57 MNG [Mongolia] 0
## 58 MUS [Mauritius] 8
## 59 MYS [Malaysia] 0
## 60 NLD [Netherlands] 0
## 61 NOR [Norway] 0
## 62 NZL [New Zealand] 0
## 63 PER [Peru] 0
## 64 PHL [Philippines] 0
## 65 POL [Poland] 0
## 66 PRT [Portugal] 0
## 67 PRY [Paraguay] 0
## 68 PSE [Palestinian Authority] 0
## 69 QAT [Qatar] 0
## 70 QCI [B-S-J-Z (China)] 0
## 71 QKI [Kurdistan Region (Iraq)] 8
## 72 QTJ [Dushanbe (Tajikistan)] 8
## 73 QUA [Ukrainian regions (17 of 27)] 0
## 74 ROU [Romania] 0
## 75 RWA [Rwanda] 8
## 76 SAU [Saudi Arabia] 0
## 77 SGP [Singapore] 0
## 78 SLV [El Salvador] 0
## 79 SRB [Serbia] 0
## 80 SVK [Slovak Republic] 0
## 81 SVN [Slovenia] 0
## 82 SWE [Sweden] 0
## 83 TAP [Chinese Taipei] 0
## 84 THA [Thailand] 0
## 85 TUR [Türkiye] 0
## 86 URY [Uruguay] 0
## 87 USA [United States] 4
## 88 UZB [Uzbekistan] 1
## 89 VNM [Viet Nam] 0
## 90 ZMB [Zambia] 8
make a table for countries NA proportion > 0.5
countries_with_high_na <- high_na_vars_per_country %>%
filter(Variables_Over_50pct_NA > 0)
print(as_tibble(countries_with_high_na), n = Inf)
## # A tibble: 10 × 2
## country Variables_Over_50pct_NA
## <chr+lbl> <dbl>
## 1 ARM [Armenia] 8
## 2 KHM [Cambodia] 8
## 3 LBN [Lebanon] 8
## 4 MUS [Mauritius] 8
## 5 QKI [Kurdistan Region (Iraq)] 8
## 6 QTJ [Dushanbe (Tajikistan)] 8
## 7 RWA [Rwanda] 8
## 8 USA [United States] 4
## 9 UZB [Uzbekistan] 1
## 10 ZMB [Zambia] 8
look into exact variables over the 0.5 threshold
agency_prop <- agency %>%
group_by(country) %>%
summarise(across(where(is.numeric), ~ mean(is.na(.x))))
t(agency_prop %>% filter(country == "ARM"))
## [,1]
## country "ARM"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "0.6203704"
## identify_gaps "0.6457966"
## self_question "0.6497648"
## adapt_strategy "0.6537331"
## ask_classmates "0.5793651"
## ask_teachers "0.5927396"
## ask_general "0.6015579"
## ask_after_trying "0.6050852"
t(agency_prop %>% filter(country == "KHM"))
## [,1]
## country "KHM"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "LBN"))
## [,1]
## country "LBN"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "0.5456516"
## identify_gaps "0.5428068"
## self_question "0.5464644"
## adapt_strategy "0.5539149"
## ask_classmates "0.5149011"
## ask_teachers "0.5239772"
## ask_general "0.5284476"
## ask_after_trying "0.5302086"
t(agency_prop %>% filter(country == "MUS"))
## [,1]
## country "MUS"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "QKI"))
## [,1]
## country "QKI"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "QTJ"))
## [,1]
## country "QTJ"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "RWA"))
## [,1]
## country "RWA"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "ZMB"))
## [,1]
## country "ZMB"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "USA"))
## [,1]
## country "USA"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "0.5018386"
## identify_gaps "0.4990266"
## self_question "0.5009734"
## adapt_strategy "0.5005408"
## ask_classmates "0.4970798"
## ask_teachers "0.4994592"
## ask_general "0.5003245"
## ask_after_trying "0.4983777"
t(agency_prop %>% filter(country == "UZB"))
## [,1]
## country "UZB"
## OECD "0"
## math_score "1"
## science_score "0"
## plan_study "0.1500231"
## identify_gaps "0.1526815"
## self_question "0.1573047"
## adapt_strategy "0.1637772"
## ask_classmates "0.1488673"
## ask_teachers "0.1508322"
## ask_general "0.1552242"
## ask_after_trying "0.1575358"
change the threshold to 0.55
library(tibble)
agency_na_prop1 <- agency %>%
group_by(country) %>%
summarise(across(where(is.numeric), ~ mean(is.na(.x))))
high_na_vars_per_country1 <- agency_na_prop1 %>%
mutate(across(-country, ~ .x > 0.55)) %>%
mutate(Variables_Over_55pct_NA1 = rowSums(select(., -country))) %>%
select(country, Variables_Over_55pct_NA1)
#View(high_na_vars_per_country)
print(as_tibble(high_na_vars_per_country1), n = Inf)
## # A tibble: 90 × 2
## country Variables_Over_55pct_NA1
## <chr+lbl> <dbl>
## 1 ALB [Albania] 0
## 2 ARE [United Arab Emirates] 0
## 3 ARG [Argentina] 0
## 4 ARM [Armenia] 8
## 5 AUS [Australia] 0
## 6 AUT [Austria] 0
## 7 AZE [Azerbaijan] 0
## 8 BEL [Belgium] 0
## 9 BGR [Bulgaria] 0
## 10 BRA [Brazil] 0
## 11 BRN [Brunei Darussalam] 0
## 12 CAN [Canada] 0
## 13 CHE [Switzerland] 0
## 14 CHL [Chile] 0
## 15 COL [Colombia] 0
## 16 CRI [Costa Rica] 0
## 17 CZE [Czechia] 0
## 18 DEU [Germany] 0
## 19 DNK [Denmark] 0
## 20 DOM [Dominican Republic] 0
## 21 ECU [Ecuador] 0
## 22 ESP [Spain] 0
## 23 EST [Estonia] 0
## 24 FIN [Finland] 0
## 25 FRA [France] 0
## 26 GBR [United Kingdom] 0
## 27 GEO [Georgia] 0
## 28 GRC [Greece] 0
## 29 GTM [Guatemala] 0
## 30 HKG [Hong Kong (China)] 0
## 31 HRV [Croatia] 0
## 32 HUN [Hungary] 0
## 33 IDN [Indonesia] 0
## 34 IRL [Ireland] 0
## 35 ISL [Iceland] 0
## 36 ISR [Israel] 0
## 37 ITA [Italy] 0
## 38 JOR [Jordan] 0
## 39 JPN [Japan] 0
## 40 KAZ [Kazakhstan] 0
## 41 KEN [Kenya] 0
## 42 KGZ [Kyrgyzstan] 0
## 43 KHM [Cambodia] 8
## 44 KOR [Korea] 0
## 45 KSV [Kosovo] 0
## 46 LBN [Lebanon] 1
## 47 LTU [Lithuania] 0
## 48 LUX [Luxembourg] 0
## 49 LVA [Latvia] 0
## 50 MAC [Macao (China)] 0
## 51 MAR [Morocco] 0
## 52 MDA [Moldova] 0
## 53 MEX [Mexico] 0
## 54 MKD [North Macedonia] 0
## 55 MLT [Malta] 0
## 56 MNE [Montenegro] 0
## 57 MNG [Mongolia] 0
## 58 MUS [Mauritius] 8
## 59 MYS [Malaysia] 0
## 60 NLD [Netherlands] 0
## 61 NOR [Norway] 0
## 62 NZL [New Zealand] 0
## 63 PER [Peru] 0
## 64 PHL [Philippines] 0
## 65 POL [Poland] 0
## 66 PRT [Portugal] 0
## 67 PRY [Paraguay] 0
## 68 PSE [Palestinian Authority] 0
## 69 QAT [Qatar] 0
## 70 QCI [B-S-J-Z (China)] 0
## 71 QKI [Kurdistan Region (Iraq)] 8
## 72 QTJ [Dushanbe (Tajikistan)] 8
## 73 QUA [Ukrainian regions (17 of 27)] 0
## 74 ROU [Romania] 0
## 75 RWA [Rwanda] 8
## 76 SAU [Saudi Arabia] 0
## 77 SGP [Singapore] 0
## 78 SLV [El Salvador] 0
## 79 SRB [Serbia] 0
## 80 SVK [Slovak Republic] 0
## 81 SVN [Slovenia] 0
## 82 SWE [Sweden] 0
## 83 TAP [Chinese Taipei] 0
## 84 THA [Thailand] 0
## 85 TUR [Türkiye] 0
## 86 URY [Uruguay] 0
## 87 USA [United States] 0
## 88 UZB [Uzbekistan] 1
## 89 VNM [Viet Nam] 0
## 90 ZMB [Zambia] 8
make a table for countries NA proportion > 0.55
countries_with_high_na1 <- high_na_vars_per_country1 %>%
filter(Variables_Over_55pct_NA1 > 0)
print(as_tibble(countries_with_high_na1), n = Inf)
## # A tibble: 9 × 2
## country Variables_Over_55pct_NA1
## <chr+lbl> <dbl>
## 1 ARM [Armenia] 8
## 2 KHM [Cambodia] 8
## 3 LBN [Lebanon] 1
## 4 MUS [Mauritius] 8
## 5 QKI [Kurdistan Region (Iraq)] 8
## 6 QTJ [Dushanbe (Tajikistan)] 8
## 7 RWA [Rwanda] 8
## 8 UZB [Uzbekistan] 1
## 9 ZMB [Zambia] 8
t(agency_prop %>% filter(country == "ARM"))
## [,1]
## country "ARM"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "0.6203704"
## identify_gaps "0.6457966"
## self_question "0.6497648"
## adapt_strategy "0.6537331"
## ask_classmates "0.5793651"
## ask_teachers "0.5927396"
## ask_general "0.6015579"
## ask_after_trying "0.6050852"
t(agency_prop %>% filter(country == "KHM"))
## [,1]
## country "KHM"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "LBN"))
## [,1]
## country "LBN"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "0.5456516"
## identify_gaps "0.5428068"
## self_question "0.5464644"
## adapt_strategy "0.5539149"
## ask_classmates "0.5149011"
## ask_teachers "0.5239772"
## ask_general "0.5284476"
## ask_after_trying "0.5302086"
t(agency_prop %>% filter(country == "MUS"))
## [,1]
## country "MUS"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "QKI"))
## [,1]
## country "QKI"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "QTJ"))
## [,1]
## country "QTJ"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "RWA"))
## [,1]
## country "RWA"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "ZMB"))
## [,1]
## country "ZMB"
## OECD "0"
## math_score "0"
## science_score "0"
## plan_study "1"
## identify_gaps "1"
## self_question "1"
## adapt_strategy "1"
## ask_classmates "1"
## ask_teachers "1"
## ask_general "1"
## ask_after_trying "1"
t(agency_prop %>% filter(country == "UZB"))
## [,1]
## country "UZB"
## OECD "0"
## math_score "1"
## science_score "0"
## plan_study "0.1500231"
## identify_gaps "0.1526815"
## self_question "0.1573047"
## adapt_strategy "0.1637772"
## ask_classmates "0.1488673"
## ask_teachers "0.1508322"
## ask_general "0.1552242"
## ask_after_trying "0.1575358"
delete rows about these countries
countries_to_exclude <- c("ARM", "KHM", "LBN", "MUS", "QKI", "QTJ", "RWA", "ZMB", "UZB")
agency1 <- agency %>%
filter(!country %in% countries_to_exclude)
now no NA proportion over the threshold of 0.5
agency_na_prop2 <- agency1 %>%
group_by(country) %>%
summarise(across(where(is.numeric), ~ mean(is.na(.x))))
#View(agency_na_prop1)
now overall, 81 countries
high_na_vars_per_country2 <- agency_na_prop2 %>%
mutate(across(-country, ~ .x > 0.55)) %>%
mutate(Variables_Over_55pct_NA2 = rowSums(select(., -country))) %>%
select(country, Variables_Over_55pct_NA2)
#View(high_na_vars_per_country1)
print(as_tibble(high_na_vars_per_country2), n = Inf)
## # A tibble: 81 × 2
## country Variables_Over_55pct_NA2
## <chr+lbl> <dbl>
## 1 ALB [Albania] 0
## 2 ARE [United Arab Emirates] 0
## 3 ARG [Argentina] 0
## 4 AUS [Australia] 0
## 5 AUT [Austria] 0
## 6 AZE [Azerbaijan] 0
## 7 BEL [Belgium] 0
## 8 BGR [Bulgaria] 0
## 9 BRA [Brazil] 0
## 10 BRN [Brunei Darussalam] 0
## 11 CAN [Canada] 0
## 12 CHE [Switzerland] 0
## 13 CHL [Chile] 0
## 14 COL [Colombia] 0
## 15 CRI [Costa Rica] 0
## 16 CZE [Czechia] 0
## 17 DEU [Germany] 0
## 18 DNK [Denmark] 0
## 19 DOM [Dominican Republic] 0
## 20 ECU [Ecuador] 0
## 21 ESP [Spain] 0
## 22 EST [Estonia] 0
## 23 FIN [Finland] 0
## 24 FRA [France] 0
## 25 GBR [United Kingdom] 0
## 26 GEO [Georgia] 0
## 27 GRC [Greece] 0
## 28 GTM [Guatemala] 0
## 29 HKG [Hong Kong (China)] 0
## 30 HRV [Croatia] 0
## 31 HUN [Hungary] 0
## 32 IDN [Indonesia] 0
## 33 IRL [Ireland] 0
## 34 ISL [Iceland] 0
## 35 ISR [Israel] 0
## 36 ITA [Italy] 0
## 37 JOR [Jordan] 0
## 38 JPN [Japan] 0
## 39 KAZ [Kazakhstan] 0
## 40 KEN [Kenya] 0
## 41 KGZ [Kyrgyzstan] 0
## 42 KOR [Korea] 0
## 43 KSV [Kosovo] 0
## 44 LTU [Lithuania] 0
## 45 LUX [Luxembourg] 0
## 46 LVA [Latvia] 0
## 47 MAC [Macao (China)] 0
## 48 MAR [Morocco] 0
## 49 MDA [Moldova] 0
## 50 MEX [Mexico] 0
## 51 MKD [North Macedonia] 0
## 52 MLT [Malta] 0
## 53 MNE [Montenegro] 0
## 54 MNG [Mongolia] 0
## 55 MYS [Malaysia] 0
## 56 NLD [Netherlands] 0
## 57 NOR [Norway] 0
## 58 NZL [New Zealand] 0
## 59 PER [Peru] 0
## 60 PHL [Philippines] 0
## 61 POL [Poland] 0
## 62 PRT [Portugal] 0
## 63 PRY [Paraguay] 0
## 64 PSE [Palestinian Authority] 0
## 65 QAT [Qatar] 0
## 66 QCI [B-S-J-Z (China)] 0
## 67 QUA [Ukrainian regions (17 of 27)] 0
## 68 ROU [Romania] 0
## 69 SAU [Saudi Arabia] 0
## 70 SGP [Singapore] 0
## 71 SLV [El Salvador] 0
## 72 SRB [Serbia] 0
## 73 SVK [Slovak Republic] 0
## 74 SVN [Slovenia] 0
## 75 SWE [Sweden] 0
## 76 TAP [Chinese Taipei] 0
## 77 THA [Thailand] 0
## 78 TUR [Türkiye] 0
## 79 URY [Uruguay] 0
## 80 USA [United States] 0
## 81 VNM [Viet Nam] 0
now delete rows with cells with missing values
agency1 <- na.omit(agency1)
anyNA(agency1)
## [1] FALSE
glimpse(agency1)
## Rows: 549,099
## Columns: 12
## $ country <chr+lbl> "ALB", "ALB", "ALB", "ALB", "ALB", "ALB", "ALB", …
## $ OECD <dbl+lbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ math_score <dbl+lbl> 328.7014, 384.2887, 438.3557, 344.1591, 376.3203,…
## $ science_score <dbl+lbl> 289.6784, 416.2245, 488.5965, 326.0128, 348.6663,…
## $ plan_study <dbl+lbl> 3, 3, 3, 4, 2, 1, 3, 2, 3, 3, 3, 2, 2, 2, 3, 3, 1…
## $ identify_gaps <dbl+lbl> 3, 3, 3, 3, 2, 2, 3, 2, 3, 3, 3, 2, 1, 2, 2, 2, 2…
## $ self_question <dbl+lbl> 3, 3, 3, 3, 2, 2, 3, 2, 3, 3, 3, 1, 2, 2, 3, 2, 1…
## $ adapt_strategy <dbl+lbl> 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 1, 3, 2, 3, 4, 1…
## $ ask_classmates <dbl+lbl> 2, 2, 2, 1, 1, 1, 1, 2, 3, 2, 2, 2, 2, 2, 2, 3, 1…
## $ ask_teachers <dbl+lbl> 2, 2, 2, 1, 1, 1, 1, 2, 3, 2, 2, 2, 2, 2, 1, 3, 1…
## $ ask_general <dbl+lbl> 2, 2, 2, 1, 1, 1, 1, 2, 3, 2, 2, 2, 1, 3, 2, 2, 2…
## $ ask_after_trying <dbl+lbl> 2, 2, 2, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 2, 1, 3, 2…
summary statistics table
numeric variables
library(dplyr)
library(tidyr)
numeric_table <- agency1 %>%
select(OECD, math_score, science_score) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
group_by(Variable) %>%
summarise(
Mean = round(mean(Value, na.rm = TRUE), 2),
Median = round(median(Value, na.rm = TRUE), 2),
SD = round(sd(Value, na.rm = TRUE), 2),
Min = round(min(Value, na.rm = TRUE), 2),
Max = round(max(Value, na.rm = TRUE), 2)
)
print(numeric_table)
## # A tibble: 3 × 6
## Variable Mean Median SD Min Max
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 OECD 0.45 0 0.5 0 1
## 2 math_score 451. 444. 104. 85.9 941.
## 3 science_score 468. 464. 105. 86.0 893.
write.csv(numeric_table, "~/Downloads/numeric_summary_statistics.csv", row.names = FALSE)
categorical variables-> treat as numeric variables
agency1 <- agency1 %>%
mutate(across(c(plan_study, identify_gaps, self_question, adapt_strategy,
ask_classmates, ask_teachers, ask_general, ask_after_trying),
as.numeric))
numeric_summary_table <- agency1 %>%
select(plan_study, identify_gaps, self_question, adapt_strategy,
ask_classmates, ask_teachers, ask_general, ask_after_trying) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
group_by(Variable) %>%
summarise(
Mean = round(mean(Value, na.rm = TRUE), 2),
Median = round(median(Value, na.rm = TRUE), 2),
SD = round(sd(Value, na.rm = TRUE), 2),
Min = round(min(Value, na.rm = TRUE), 2),
Max = round(max(Value, na.rm = TRUE), 2)
)
print(numeric_summary_table, n = Inf)
## # A tibble: 8 × 6
## Variable Mean Median SD Min Max
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 adapt_strategy 2.52 2 0.9 1 4
## 2 ask_after_trying 1.99 2 0.79 1 4
## 3 ask_classmates 1.96 2 0.76 1 4
## 4 ask_general 2.01 2 0.75 1 4
## 5 ask_teachers 2.02 2 0.76 1 4
## 6 identify_gaps 2.33 2 0.89 1 4
## 7 plan_study 2.42 2 0.91 1 4
## 8 self_question 2.57 3 0.9 1 4
write.csv(numeric_summary_table, "~/Downloads/numeric_summary_statistics_agency.csv", row.names = FALSE)
total row: 549,089
nrow(agency1)
## [1] 549099
EDA
Plot 1: perceived agency for goal setting and help seeking: Lollipop
plot
library(dplyr)
library(tidyr)
library(ggplot2)
country_means <- agency1 %>%
group_by(country) %>%
summarise(
plan_study = mean(plan_study, na.rm = TRUE),
identify_gaps = mean(identify_gaps, na.rm = TRUE),
self_question = mean(self_question, na.rm = TRUE),
adapt_strategy = mean(adapt_strategy, na.rm = TRUE),
ask_classmates = mean(ask_classmates, na.rm = TRUE),
ask_teachers = mean(ask_teachers, na.rm = TRUE),
ask_general = mean(ask_general, na.rm = TRUE),
ask_after_trying = mean(ask_after_trying, na.rm = TRUE)
)
# 2. OECD average
oecd_country_means <- agency1 %>%
filter(OECD == 1) %>%
group_by(country) %>%
summarise(
plan_study = mean(plan_study, na.rm = TRUE),
identify_gaps = mean(identify_gaps, na.rm = TRUE),
self_question = mean(self_question, na.rm = TRUE),
adapt_strategy = mean(adapt_strategy, na.rm = TRUE),
ask_classmates = mean(ask_classmates, na.rm = TRUE),
ask_teachers = mean(ask_teachers, na.rm = TRUE),
ask_general = mean(ask_general, na.rm = TRUE),
ask_after_trying = mean(ask_after_trying, na.rm = TRUE)
)
oecd_means <- oecd_country_means %>%
summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE))) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "OECD_Mean")
# long format
long_means <- country_means %>%
pivot_longer(
cols = -country,
names_to = "Variable",
values_to = "Mean_Value"
)
# mark the countries with highest and lowest value for each variable
extremes <- long_means %>%
group_by(Variable) %>%
filter(Mean_Value == max(Mean_Value, na.rm = TRUE) |
Mean_Value == min(Mean_Value, na.rm = TRUE)) %>%
mutate(Type = if_else(Mean_Value == max(Mean_Value, na.rm = TRUE), "Highest", "Lowest"))
var_levels <- c("plan_study", "identify_gaps", "self_question", "adapt_strategy",
"ask_classmates", "ask_teachers", "ask_general", "ask_after_trying")
long_means$Variable <- factor(long_means$Variable, levels = var_levels)
extremes$Variable <- factor(extremes$Variable, levels = var_levels)
oecd_means$Variable <- factor(oecd_means$Variable, levels = var_levels)
# make the y-axis start from 1
baseline_val <- 1
p <- ggplot() +
geom_segment(data = long_means,
aes(x = Variable, xend = Variable, y = baseline_val, yend = Mean_Value, group = country),
color = "grey75", linewidth = 0.4, alpha = 0.4) +
geom_segment(data = extremes,
aes(x = Variable, xend = Variable, y = baseline_val, yend = Mean_Value, group = country),
color = "#2b83ba", linewidth = 0.8) +
geom_point(data = data.frame(Variable = factor(var_levels, levels = var_levels)),
aes(x = Variable, y = baseline_val),
color = "#1f78b4", shape = 16, size = 3.5) +
geom_point(data = oecd_means,
aes(x = Variable, y = OECD_Mean),
color = "#1f78b4", shape = 18, size = 5) +
geom_text(data = oecd_means,
aes(x = Variable, y = OECD_Mean, label = "OECD average"),
hjust = -0.2, vjust = 0.5, size = 3, fontface = "bold", color = "#1f78b4") +
geom_point(data = extremes,
aes(x = Variable, y = Mean_Value, color = Type),
shape = 15, size = 3) +
geom_text(data = extremes,
aes(x = Variable, y = Mean_Value, label = country, color = Type),
hjust = -0.15, vjust = 0.5, size = 3, show.legend = FALSE) +
scale_color_manual(values = c("Highest" = "#d7191c", "Lowest" = "#2b83ba")) +
scale_y_continuous(limits = c(1, NA), breaks = seq(1, 4, by = 1), expand = expansion(mult = c(0.05, 0.22))) +
labs(
title = "Goal Setting & Help Seeking: Perceived Agency",
x = NULL,
y = "Mean index",
color = "Country Extremes"
) +
theme_minimal(base_size = 11) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 25, hjust = 1, face = "bold", size = 9),
axis.text.y = element_text(size = 9),
axis.line.y = element_line(color = "black"),
legend.position = "bottom"
)
based on this Lollipop graph, change the country names into full
formats
country_name_map <- c(
"KEN" = "Kenya",
"CZE" = "Czech Republic",
"ALB" = "Albania",
"UPN" = "Japan",
"JPN" = "Japan",
"GRC" = "Greece",
"BRN" = "Brunei",
"POL" = "Poland",
"QCI" = "B-S-J-Z (China)",
"SGP" = "Singapore"
)
extremes <- extremes %>%
mutate(country = recode(country, !!!country_name_map))
p <- ggplot() +
geom_segment(data = long_means,
aes(x = Variable, xend = Variable, y = baseline_val, yend = Mean_Value, group = country),
color = "grey75", linewidth = 0.4, alpha = 0.4) +
geom_segment(data = extremes,
aes(x = Variable, xend = Variable, y = baseline_val, yend = Mean_Value, group = country),
color = "#2b83ba", linewidth = 0.8) +
geom_point(data = data.frame(Variable = factor(var_levels, levels = var_levels)),
aes(x = Variable, y = baseline_val),
color = "#1f78b4", shape = 16, size = 3.5) +
geom_point(data = oecd_means,
aes(x = Variable, y = OECD_Mean),
color = "#238b45", shape = 18, size = 5) +
geom_text(data = oecd_means,
aes(x = Variable, y = OECD_Mean, label = "OECD average"),
hjust = -0.2, vjust = 0.5, size = 3, fontface = "bold", color = "#238b45") +
geom_point(data = extremes,
aes(x = Variable, y = Mean_Value, color = Type),
shape = 15, size = 3) +
geom_text(data = extremes,
aes(x = Variable, y = Mean_Value, label = country, color = Type),
hjust = -0.15, vjust = 0.5, size = 3, show.legend = FALSE) +
scale_color_manual(values = c("Highest" = "#d7191c", "Lowest" = "#2b83ba")) +
scale_y_continuous(limits = c(1, NA), breaks = seq(1, 4, by = 1), expand = expansion(mult = c(0.05, 0.28))) +
labs(
title = "Goal Setting & Help Seeking: Perceived Agency",
x = NULL,
y = "Mean index",
color = "Country Extremes"
) +
theme_minimal(base_size = 11) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 25, hjust = 1, face = "bold", size = 9),
axis.text.y = element_text(size = 9),
axis.line.y = element_line(color = "black"),
legend.position = "bottom"
)
print(p)
ggsave(filename = "~/Downloads/lollipop_plot.png", plot = p, width = 14, height = 7, dpi = 300)
Plot 2: goal setting: dumbbell plot
country_name_map <- c(
"ALB" = "Albania",
"ARE" = "United Arab Emirates",
"ARG" = "Argentina",
"AUS" = "Australia",
"AUT" = "Austria",
"AZE" = "Azerbaijan",
"BEL" = "Belgium",
"BGR" = "Bulgaria",
"BRA" = "Brazil",
"BRN" = "Brunei Darussalam",
"CAN" = "Canada",
"CHE" = "Switzerland",
"CHL" = "Chile",
"COL" = "Colombia",
"CRI" = "Costa Rica",
"CZE" = "Czechia",
"DEU" = "Germany",
"DNK" = "Denmark",
"DOM" = "Dominican Republic",
"ECU" = "Ecuador",
"ESP" = "Spain",
"EST" = "Estonia",
"FIN" = "Finland",
"FRA" = "France",
"GBR" = "United Kingdom",
"GEO" = "Georgia",
"GRC" = "Greece",
"GTM" = "Guatemala",
"HKG" = "Hong Kong (China)",
"HRV" = "Croatia",
"HUN" = "Hungary",
"IDN" = "Indonesia",
"IRL" = "Ireland",
"ISL" = "Iceland",
"ISR" = "Israel",
"ITA" = "Italy",
"JOR" = "Jordan",
"JPN" = "Japan",
"KAZ" = "Kazakhstan",
"KEN" = "Kenya",
"KGZ" = "Kyrgyzstan",
"KOR" = "Korea",
"KSV" = "Kosovo",
"LTU" = "Lithuania",
"LUX" = "Luxembourg",
"LVA" = "Latvia",
"MAC" = "Macao (China)",
"MAR" = "Morocco",
"MDA" = "Moldova",
"MEX" = "Mexico",
"MKD" = "North Macedonia",
"MLT" = "Malta",
"MNE" = "Montenegro",
"MNG" = "Mongolia",
"MYS" = "Malaysia",
"NLD" = "Netherlands",
"NOR" = "Norway",
"NZL" = "New Zealand",
"PER" = "Peru",
"PHL" = "Philippines",
"POL" = "Poland",
"PRT" = "Portugal",
"PRY" = "Paraguay",
"PSE" = "Palestinian Authority",
"QAT" = "Qatar",
"QCI" = "B-S-J-Z (China)",
"QCY" = "Cyprus",
"QUA" = "Ukrainian regions (17 of 27)",
"ROU" = "Romania",
"SAU" = "Saudi Arabia",
"SGP" = "Singapore",
"SLV" = "El Salvador",
"SRB" = "Serbia",
"SVK" = "Slovak Republic",
"SVN" = "Slovenia",
"SWE" = "Sweden",
"TAP" = "Chinese Taipei",
"THA" = "Thailand",
"TUR" = "Türkiye",
"URY" = "Uruguay",
"USA" = "United States",
"VNM" = "Viet Nam"
)
# create the goal setting variable
agency1 <- agency1 %>%
mutate(goal_setting = plan_study + identify_gaps + self_question + adapt_strategy)
# figure out the bottom quarter and top quarter
country_quarter_means <- agency1 %>%
group_by(country) %>%
mutate(
goal_quartile = ntile(goal_setting, 4)
) %>%
filter(goal_quartile %in% c(1, 4)) %>%
group_by(country, goal_quartile) %>%
summarise(
mean_val = mean(goal_setting, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_wider(
names_from = goal_quartile,
values_from = mean_val,
names_prefix = "q_"
) %>%
rename(Bottom_Quarter = q_1, Top_Quarter = q_4) %>%
mutate(country = recode(country, !!!country_name_map))
# OECD average
oecd_subset <- agency1 %>%
filter(OECD == 1) %>%
mutate(goal_quartile = ntile(goal_setting, 4)) %>%
filter(goal_quartile %in% c(1, 4)) %>%
group_by(goal_quartile) %>%
summarise(
mean_val = mean(goal_setting, na.rm = TRUE)
)
oecd_row <- data.frame(
country = "OECD average",
Bottom_Quarter = oecd_subset$mean_val[oecd_subset$goal_quartile == 1],
Top_Quarter = oecd_subset$mean_val[oecd_subset$goal_quartile == 4]
)
plot_wide <- bind_rows(country_quarter_means, oecd_row)
# order countries by Top quarter ranking
plot_wide <- plot_wide %>%
mutate(country = factor(country, levels = country[order(Top_Quarter, decreasing = FALSE)]))
non_oecd_levels <- setdiff(levels(plot_wide$country), "OECD average")
bottom_5_countries <- head(non_oecd_levels, 5)
top_5_countries <- tail(non_oecd_levels, 5)
plot_wide <- plot_wide %>%
mutate(
row_group = case_when(
country == "OECD average" ~ "OECD",
country %in% top_5_countries ~ "Top 5",
country %in% bottom_5_countries ~ "Bottom 5",
TRUE ~ "Other"
)
)
plot_long <- plot_wide %>%
pivot_longer(
cols = c(Bottom_Quarter, Top_Quarter),
names_to = "Quartile",
values_to = "Score"
)
# add the shaded rows for OECD average (light blue), top 5, and bottom 5 (light green)
p <- ggplot(plot_wide, aes(y = country)) +
geom_rect(data = filter(plot_wide, row_group == "Bottom 5"),
aes(xmin = -Inf, xmax = Inf, ymin = as.numeric(country) - 0.5, ymax = as.numeric(country) + 0.5),
fill = "#e5f5e0", alpha = 0.6, inherit.aes = FALSE) +
geom_rect(data = filter(plot_wide, row_group == "Top 5"),
aes(xmin = -Inf, xmax = Inf, ymin = as.numeric(country) - 0.5, ymax = as.numeric(country) + 0.5),
fill = "#e5f5e0", alpha = 0.6, inherit.aes = FALSE) +
geom_rect(data = filter(plot_wide, country == "OECD average"),
aes(xmin = -Inf, xmax = Inf, ymin = as.numeric(country) - 0.5, ymax = as.numeric(country) + 0.5),
fill = "#c6dbef", alpha = 0.6, inherit.aes = FALSE) +
geom_segment(aes(x = Bottom_Quarter, xend = Top_Quarter, y = country, yend = country),
color = "grey60", linewidth = 0.8) +
geom_point(data = plot_long, aes(x = Score, color = Quartile, shape = Quartile), size = 3) +
scale_color_manual(
values = c("Bottom_Quarter" = "#e6550d", "Top_Quarter" = "#3182bd"),
labels = c("Bottom quarter of students' goal setting", "Top quarter of students' goal setting")
) +
scale_shape_manual(
values = c("Bottom_Quarter" = 16, "Top_Quarter" = 18),
labels = c("Bottom quarter of students' goal setting", "Top quarter of students' goal setting")
) +
labs(
title = "Country-Level Student Goal-Setting Index by Quartile",
x = "Goal-Setting Index Score",
y = NULL,
color = NULL,
shape = NULL
) +
theme_minimal(base_size = 11) +
theme(
panel.grid.major.y = element_line(color = "grey90"),
panel.grid.minor = element_blank(),
axis.text.y = element_text(
size = 8.5,
face = "bold",
color = ifelse(levels(plot_wide$country) == "OECD average", "#08519c", "black")
),
legend.position = "bottom"
)
print(p)
ggsave(filename = "~/Downloads/dumbbell_plot_shaded.png", plot = p, width = 9, height = 11, dpi = 300)
Plot 3: Help seeking: dot plots for 4 variables
country_means <- agency1 %>%
group_by(country) %>%
summarise(
`Ask Classmates` = mean(ask_classmates, na.rm = TRUE),
`Ask Teachers` = mean(ask_teachers, na.rm = TRUE),
`Ask General` = mean(ask_general, na.rm = TRUE),
`Ask After Trying` = mean(ask_after_trying, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(country = recode(country, !!!country_name_map))
# OECD average
oecd_means <- agency1 %>%
filter(OECD == 1) %>%
summarise(
`Ask Classmates` = mean(ask_classmates, na.rm = TRUE),
`Ask Teachers` = mean(ask_teachers, na.rm = TRUE),
`Ask General` = mean(ask_general, na.rm = TRUE),
`Ask After Trying` = mean(ask_after_trying, na.rm = TRUE)
) %>%
mutate(country = "OECD average")
plot_data <- bind_rows(country_means, oecd_means)
plot_data <- plot_data %>%
mutate(
overall_mean = (`Ask Classmates` + `Ask Teachers` + `Ask General` + `Ask After Trying`) / 4,
country = factor(country, levels = country[order(overall_mean, decreasing = FALSE)])
)
plot_long <- plot_data %>%
pivot_longer(
cols = c(`Ask Classmates`, `Ask Teachers`, `Ask General`, `Ask After Trying`),
names_to = "Help_Seeking_Type",
values_to = "Mean_Score"
)
variable_highlights <- plot_long %>%
filter(country != "OECD average") %>%
group_by(Help_Seeking_Type) %>%
arrange(Mean_Score) %>%
mutate(
rank_val = row_number(),
max_rank = max(rank_val),
highlight = case_when(
rank_val <= 3 ~ "Bottom 3",
rank_val >= max_rank - 2 ~ "Top 3",
TRUE ~ "Other"
)
) %>%
ungroup()
oecd_data <- plot_long %>%
filter(country == "OECD average")
# 4. for each plot, mark top/bottom 3 shading countries (light green) and OECD shading (light blue)
p <- ggplot(plot_long, aes(x = Mean_Score, y = country, color = Help_Seeking_Type)) +
geom_rect(data = filter(variable_highlights, highlight == "Bottom 3"),
aes(xmin = -Inf, xmax = Inf, ymin = as.numeric(country) - 0.5, ymax = as.numeric(country) + 0.5),
fill = "#e5f5e0", alpha = 0.6, inherit.aes = FALSE) +
geom_rect(data = filter(variable_highlights, highlight == "Top 3"),
aes(xmin = -Inf, xmax = Inf, ymin = as.numeric(country) - 0.5, ymax = as.numeric(country) + 0.5),
fill = "#e5f5e0", alpha = 0.6, inherit.aes = FALSE) +
geom_rect(data = oecd_data,
aes(xmin = -Inf, xmax = Inf, ymin = as.numeric(country) - 0.5, ymax = as.numeric(country) + 0.5),
fill = "#c6dbef", alpha = 0.6, inherit.aes = FALSE) +
geom_point(size = 2.5) +
facet_wrap(~ Help_Seeking_Type, scales = "free_x", nrow = 1) +
scale_color_brewer(palette = "Dark2") +
labs(
title = "Country_Level Student Help-Seeking Profiles",
x = "Mean Index Score (Ordered by Variable Value)",
y = NULL,
color = NULL
) +
theme_minimal(base_size = 11) +
theme(
panel.grid.major.y = element_line(color = "grey90"),
panel.grid.minor = element_blank(),
axis.text.y = element_text(
size = 7.5,
face = "bold",
color = ifelse(levels(plot_data$country) == "OECD average", "#08519c", "black")
),
strip.text = element_text(face = "bold", size = 10),
legend.position = "none"
)
print(p)
ggsave(filename = "~/Downloads/dumbbell_help_shaded.png", plot = p, width = 14, height = 11, dpi = 300)
Plot 4: scatterplot of math and science score
library(dplyr)
library(ggplot2)
library(ggrepel)
country_scores <- agency1 %>%
group_by(country) %>%
summarise(
math_score = mean(math_score, na.rm = TRUE),
science_score = mean(science_score, na.rm = TRUE),
.groups = "drop"
) %>%
filter(!is.na(math_score) & !is.na(science_score)) %>%
mutate(country = recode(country, !!!country_name_map))
# OECD average
oecd_scores <- agency1 %>%
filter(OECD == 1) %>%
summarise(
math_score = mean(math_score, na.rm = TRUE),
science_score = mean(science_score, na.rm = TRUE)
)
oecd_math <- oecd_scores$math_score
oecd_science <- oecd_scores$science_score
# mark countries in the first quadrant, whose math and science score are higher than OECD mean
country_scores <- country_scores %>%
mutate(
quadrant = case_when(
math_score >= oecd_math & science_score >= oecd_science ~ "Above OECD (Both)",
TRUE ~ "Other"
)
)
p <- ggplot(country_scores, aes(x = math_score, y = science_score, label = country)) +
geom_vline(xintercept = oecd_math, color = "#3182bd", linetype = "dashed", linewidth = 0.8) +
geom_hline(yintercept = oecd_science, color = "#3182bd", linetype = "dashed", linewidth = 0.8) +
annotate("text", x = oecd_math, y = max(country_scores$science_score, na.rm = TRUE),
label = paste0("OECD average (", round(oecd_math, 1), ")"), color = "#08519c", fontface = "bold", size = 2,
angle = 90, vjust = -0.5, hjust = 0.9) +
annotate("text", x = max(country_scores$math_score, na.rm = TRUE), y = oecd_science,
label = paste0("OECD average (", round(oecd_science, 1), ")"), color = "#08519c", fontface = "bold", size = 2,
vjust = -0.8, hjust = 0.7) +
geom_point(aes(color = quadrant), size = 2.5, shape = 18) +
scale_color_manual(values = c("Above OECD (Both)" = "#e6550d", "Other" = "#1f78b4")) +
geom_text_repel(size = 3, box.padding = 0.3, point.padding = 0.3, max.overlaps = Inf) +
labs(
title = "Country-Level Mean Performance in Math and Science",
x = "Mean Math Score",
y = "Mean Science Score"
) +
theme_minimal(base_size = 11) +
theme(
panel.grid.minor = element_blank(),
axis.text = element_text(size = 9),
plot.title = element_text(face = "bold", size = 12),
legend.position = "none"
)
print(p)
ggsave(filename = "/Users/xutianyu/Downloads/country_math_science_scatter.png", plot = p, width = 12, height = 20, dpi = 300)
Plot 5: HCA heatmaps: uncentered correlation + average linkage
load the package
library(devtools)
## Loading required package: usethis
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.25.3
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(circlize)
## ========================================
## circlize version 0.4.18
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================
library(hopach)
## Loading required package: cluster
## Loading required package: Biobase
## Loading required package: BiocGenerics
## Loading required package: generics
##
## Attaching package: 'generics'
## The following object is masked from 'package:dplyr':
##
## explain
## The following objects are masked from 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, union
##
## Attaching package: 'BiocGenerics'
## The following object is masked from 'package:dplyr':
##
## combine
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
## unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'hopach'
## The following object is masked from 'package:generics':
##
## prune
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("hopach", update = FALSE, ask = FALSE)
## Bioconductor version 3.22 (BiocManager 1.30.27), R 4.5.3 (2026-03-11)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'hopach'
selected_vars <- c("plan_study", "identify_gaps", "self_question", "adapt_strategy",
"ask_classmates", "ask_teachers", "ask_general", "ask_after_trying")
country_data <- agency1 %>%
group_by(country) %>%
summarise(
across(all_of(selected_vars), ~ mean(.x, na.rm = TRUE)),
math_score = mean(math_score, na.rm = TRUE),
science_score = mean(science_score, na.rm = TRUE)
)
# change into full country names
country_names_map <- country_name_map
country_data$Country_Name <- unname(country_names_map[country_data$country])
country_data$Country_Name[is.na(country_data$Country_Name)] <- country_data$country[is.na(country_data$Country_Name)]
# scale the dataset
data_matrix <- as.matrix(country_data %>% select(all_of(selected_vars)))
rownames(data_matrix) <- country_data$Country_Name
data.scale.subset <- scale(data_matrix)
# distance measures: uncentered correlation
uncenter.dist <- function(m) {
as.dist(as.matrix(hopach::distancematrix(m, d = "cosangle")))
}
# agglomeration: average linkage
row.clus <- hclust(uncenter.dist(data.scale.subset), method = "ave")
col.clus <- hclust(uncenter.dist(t(data.scale.subset)), method = "ave")
# add annotation of math and science score on the right side
right_anno <- rowAnnotation(
Math = country_data$math_score,
Science = country_data$science_score,
col = list(
Math = colorRamp2(c(min(country_data$math_score, na.rm = TRUE),
mean(country_data$math_score, na.rm = TRUE),
max(country_data$math_score, na.rm = TRUE)),
c("#edf8e5", "#a1d99b", "#238b45")),
Science = colorRamp2(c(min(country_data$science_score, na.rm = TRUE),
mean(country_data$science_score, na.rm = TRUE),
max(country_data$science_score, na.rm = TRUE)),
c("#f7fbff", "#a6bddb", "#2b8cbe"))
),
simple_anno_size = unit(0.35, "cm"),
annotation_name_gp = gpar(fontsize = 8)
)
png("~/Downloads/agency_heatmap.png", width = 12, height = 12, units = "in", res = 300)
Heatmap(
data.scale.subset,
name = "agency heatmap",
cluster_rows = row.clus,
cluster_columns = col.clus,
row_split = 4,
column_split = 2,
top_annotation = HeatmapAnnotation(
foo = anno_block(
gp = gpar(fill = c("#E41A1C", "#377EB8")),
labels = c("Goal Setting", "Help Seeking"),
labels_gp = gpar(col = "white", fontsize = 10)
)
),
left_annotation = rowAnnotation(
foo = anno_block(
gp = gpar(fill = 2:5),
labels = c("Low Goal / High Help", "Low Goal / Low Help", "Moderate Goal / Moderate Help", "High Goal / Low Help"),
labels_gp = gpar(col = "white", fontsize = 8)
)
),
right_annotation = right_anno,
show_row_names = TRUE,
row_names_side = "right",
row_names_gp = gpar(fontsize = 7),
column_names_gp = gpar(fontsize = 9),
column_title = "Hierarchical Clustering Analysis of perceived agency PISA 2025",
column_title_gp = gpar(fontsize = 12, fontface = "bold")
)