本段程式碼載入這次需要的 R 套件，再用 read_spss()
函數讀取 PISA 2022 的資料檔，儲存為 PISA.dat。
library(haven)
library(psych)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:psych':
##
## logitPISA.dat <- read_spss("PISA_tawian2022_trimmed_lab.sav")本次分析要檢驗臺灣男生和女生的數學成績（PV1MATH）有沒有差異，Gender 中 1 是女生、2 是男生。題目沒有指定方向，所以用雙尾檢定，顯著水準 \(\alpha = .05\)。
做 t 檢定之前，要先用 Levene 檢定看兩組的變異數是否相等，再決定要用哪一種 t 檢定。
本段程式碼先把 PV1MATH
依性別分成女生（math_F）和男生（math_M）兩組，並排除缺失值。接著用
length() 算人數、sum()
算總和，自己算出兩組的平均數、變異數和標準差。
gender <- as.numeric(PISA.dat$Gender)
math <- as.numeric(PISA.dat$PV1MATH)
math_F <- math[gender == 1 & !is.na(gender) & !is.na(math)]
math_M <- math[gender == 2 & !is.na(gender) & !is.na(math)]
n_F <- length(math_F)
n_M <- length(math_M)
M_F <- sum(math_F) / n_F
M_M <- sum(math_M) / n_M
V_F <- sum((math_F - M_F)^2) / (n_F - 1)
V_M <- sum((math_M - M_M)^2) / (n_M - 1)
SD_F <- sqrt(V_F)
SD_M <- sqrt(V_M)
data.frame(Gender = c("女生", "男生"),
n = c(n_F, n_M),
M = c(M_F, M_M),
SD = c(SD_F, SD_M))## Gender n M SD
## 1 女生 2780 535.8361 102.2562
## 2 男生 2819 539.7803 116.6183本段程式碼使用 psych 套件中的 describeBy()
函數，依性別計算 PV1MATH 的描述統計（group 1 是女生、group 2
是男生），可以用來對照上面自己算的人數、平均數和標準差。
describeBy(PISA.dat$PV1MATH, PISA.dat$Gender)##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew
## X1 1 2780 535.84 102.26 538.56 536.99 108.38 227.79 821.56 593.77 -0.11
## kurtosis se
## X1 -0.4 1.94
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew
## X1 1 2819 539.78 116.62 546.54 541.26 125.62 219.06 917.38 698.32 -0.1
## kurtosis se
## X1 -0.51 2.2Levene 檢定是先算每位學生和自己那一組平均數的距離（取絕對值，記為 \(Z\)），再看兩組的平均距離差多少。F 值是組間平方和除以組內平方和（各自除以自由度）：
\[F = \frac{SS_{\text{組間}} / 1}{SS_{\text{組內}} / (N - 2)}\]
Z_F <- abs(math_F - M_F)
Z_M <- abs(math_M - M_M)
Zbar_F <- sum(Z_F) / n_F
Zbar_M <- sum(Z_M) / n_M
N <- n_F + n_M
Zbar <- (sum(Z_F) + sum(Z_M)) / N
SS_between <- n_F * (Zbar_F - Zbar)^2 + n_M * (Zbar_M - Zbar)^2
SS_within <- sum((Z_F - Zbar_F)^2) + sum((Z_M - Zbar_M)^2)
Lev_df1 <- 2 - 1
Lev_df2 <- N - 2
Lev_F <- (SS_between / Lev_df1) / (SS_within / Lev_df2)
Lev_F## [1] 59.34202F 檢定只看右尾，所以用 1 - pf() 算 p 值，自由度為 1 和
\(N - 2\)。
Lev_P <- 1 - pf(Lev_F, df1 = Lev_df1, df2 = Lev_df2)
Lev_P## [1] 1.554312e-14Lev_P 若顯示為 0，並不代表 p 值真的等於 0，而是 p
值過於小。
本段程式碼先用 as.factor() 把 Gender 轉成類別變項，再用
car 套件中的 leveneTest() 函數做 Levene
檢定，最後把自己算的和套件算的結果放在一起比較。p
值非常小的時候，最後幾位數字可能會有一點點不同，不影響結果。
PISA.dat$Gender <- as.factor(PISA.dat$Gender)
lev <- leveneTest(PISA.dat$PV1MATH ~ PISA.dat$Gender, center = "mean")
lev## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 1 59.342 1.557e-14 ***
## 5597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1data.frame(F_manual = Lev_F, F_package = lev[1, "F value"],
p_manual = Lev_P, p_package = lev[1, "Pr(>F)"])## F_manual F_package p_manual p_package
## 1 59.34202 59.34202 1.554312e-14 1.55731e-14Levene 檢定 \(p \geq .05\) 時用
Student’s t 檢定，\(p < .05\) 時用
Welch t 檢定。equal_var 是 TRUE 代表用
Student’s t 檢定。
equal_var <- Lev_P >= .05
equal_var## [1] FALSE假設兩組變異數相等時，先算合併變異數 \(s_p^2\)，再算標準誤 \(SE\) 和 t 值，自由度為 \(n_{\text{女}} + n_{\text{男}} - 2\)：
\[s_p^2 = \frac{(n_{\text{女}} - 1)s^2_{\text{女}} + (n_{\text{男}} - 1)s^2_{\text{男}}}{n_{\text{女}} + n_{\text{男}} - 2}, \qquad SE = \sqrt{s_p^2 (1/n_{\text{女}} + 1/n_{\text{男}})}, \qquad t = \frac{M_{\text{女}} - M_{\text{男}}}{SE}\]
Sp2 <- ((n_F - 1) * V_F + (n_M - 1) * V_M) / (n_F + n_M - 2)
SE_S <- sqrt(Sp2 * (1 / n_F + 1 / n_M))
t_S <- (M_F - M_M) / SE_S
df_S <- n_F + n_M - 2
t_S## [1] -1.344879df_S## [1] 5597假設兩組變異數不相等時，兩組的變異數分開算，自由度也要另外調整（公式比較長，直接寫在程式碼的
df_W）：
\[SE = \sqrt{s^2_{\text{女}} / n_{\text{女}} + s^2_{\text{男}} / n_{\text{男}}}, \qquad t = \frac{M_{\text{女}} - M_{\text{男}}}{SE}\]
SE_W <- sqrt(V_F / n_F + V_M / n_M)
t_W <- (M_F - M_M) / SE_W
df_W <- (V_F / n_F + V_M / n_M)^2 /
((V_F / n_F)^2 / (n_F - 1) + (V_M / n_M)^2 / (n_M - 1))
t_W## [1] -1.346104df_W## [1] 5521.601和上次單一樣本 t 檢定一樣，用 2 * (1 - pt(abs(t), df))
算雙尾 p 值。95% 信賴區間是兩組平均數的差，加減
qt(.975, df) 乘上 \(SE\)。
P_S <- 2 * (1 - pt(abs(t_S), df = df_S))
P_W <- 2 * (1 - pt(abs(t_W), df = df_W))
CI_S <- (M_F - M_M) + c(-1, 1) * qt(.975, df = df_S) * SE_S
CI_W <- (M_F - M_M) + c(-1, 1) * qt(.975, df = df_W) * SE_W
P_S## [1] 0.1787187P_W## [1] 0.1783242CI_S## [1] -9.693658 1.805158CI_W## [1] -9.688443 1.799944p 值若顯示為 0，並不代表 p 值真的等於 0，而是 p 值過於小。
本段程式碼使用 t.test() 函數做獨立樣本 t
檢定，var.equal = TRUE 是 Student’s t
檢定，var.equal = FALSE 是 Welch t 檢定，最後把自己算的和
t.test() 的結果放在一起比較。
tt_S <- t.test(PISA.dat$PV1MATH ~ PISA.dat$Gender, var.equal = TRUE)
tt_W <- t.test(PISA.dat$PV1MATH ~ PISA.dat$Gender, var.equal = FALSE)
tt_S##
## Two Sample t-test
##
## data: PISA.dat$PV1MATH by PISA.dat$Gender
## t = -1.3449, df = 5597, p-value = 0.1787
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -9.693658 1.805158
## sample estimates:
## mean in group 1 mean in group 2
## 535.8361 539.7803tt_W##
## Welch Two Sample t-test
##
## data: PISA.dat$PV1MATH by PISA.dat$Gender
## t = -1.3461, df = 5521.6, p-value = 0.1783
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -9.688443 1.799944
## sample estimates:
## mean in group 1 mean in group 2
## 535.8361 539.7803data.frame(test = c("Student", "Welch"),
t_manual = c(t_S, t_W),
t_package = unname(c(tt_S$statistic, tt_W$statistic)),
df_manual = c(df_S, df_W),
df_package = unname(c(tt_S$parameter, tt_W$parameter)),
p_manual = c(P_S, P_W),
p_package = c(tt_S$p.value, tt_W$p.value))## test t_manual t_package df_manual df_package p_manual p_package
## 1 Student -1.344879 -1.344879 5597.000 5597.000 0.1787187 0.1787187
## 2 Welch -1.346104 -1.346104 5521.601 5521.601 0.1783242 0.1783242因為樣本很大，就算平均數只差一點點也可能顯著，所以再算 Cohen’s \(d\) 看差異有多大：
\[d = \frac{M_{\text{女}} - M_{\text{男}}}{s_p}\]
Sp <- sqrt(Sp2)
Cohen_d <- (M_F - M_M) / Sp
Cohen_d## [1] -0.03594749本段程式碼使用 psych 套件中的 cohen.d()
函數計算效果量。cohen.d()
是用男生減女生，所以正負號會和自己算的相反，數值也會有一點點差異。
cd <- cohen.d(PISA.dat$PV1MATH, PISA.dat$Gender)
cd## Call: cohen.d(x = PISA.dat$PV1MATH, group = PISA.dat$Gender)
## Cohen d statistic of difference between two means
## lower effect upper
## [1,] -0.02 0.04 0.09
##
## Multivariate (Mahalanobis) distance between groups
## [1] 0.036
## r equivalent of difference between two means
## data
## 0.02cd$cohen.d## lower effect upper
## [1,] -0.01643992 0.03595391 0.08834502