1 資料庫介紹

本研究關注臺灣家庭脈絡下之孝道價值觀差異，並需兼顧全體成年人口的代表性與推論效度，因此選用「臺灣社會變遷基本調查（2021 年家庭組）」作為分析資料。該調查由中央研究院社會學研究所執行，於2021年以分層隨機抽樣方式，透過面對面訪問蒐集臺灣地區18歲以上一般民眾之資料；其題項涵蓋家庭結構、代間關係與相關社會態度指標，能有效對應本研究之研究問題與實證分析需求。

2 資料清理

2.1 依變項（孝道責任感指標_原始）：

將題組’E1’每個題目（a～i）中的特定遺漏值代碼(93)、(97)、(98)做遺漏值處理，並將原本的 0 ~ 4 分改為 1 ~ 5 分作為計算分數。將「有效回答題數加總分數」除以「有效回答題數」，作為孝道責任感指標，越高表示孝道責任感越強。問卷中的E1題組原始題目如下： a. 對父母的養育之恩心存感激。 b. 無論父母對您如何不好，仍然善待他們。 c. 放棄個人志向，達成父母的心願。 d. 結婚後和父母 (公婆) 住在一起。 e. 奉養父母使他們生活更舒適。 f. 父母去世，不管住得多遠，都親自奔喪。 g. 為了顧及父母的面子，為父母說些好話。 h. 為了傳宗接代，至少生一個兒子。 i. 做些讓家族感到光彩的事。

2.1.1 依變項（孝道責任感指標_原始）變項敘述統計

根據敘述統計結果可得知，臺灣民眾孝道責任感（原始）的分數平均為3.868，最小值為1.000，最大值為5.000，第一四分位數為3.444，中位數為3.889，第三四分位數為4.444；而遺漏值為2。

summary(data$E1_ori)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   3.444   3.889   3.868   4.444   5.000       2

2.1.2 孝道責任感指標（原始）：統計圖表

孝道責任感指標（原始）的統計圖表以箱型圖呈現；圖表顯示整體樣本的孝道責任感集中在中高區間，且整體分數分佈偏高，且離群值多為低分者。

plot_dvo <- ggplot(data, aes(y = E1_ori)) +
  geom_boxplot(fill = "#FFCC99", colour = "#CC9966", alpha = 1,
               width = .5, na.rm = TRUE) +
  scale_x_discrete(expand = expansion(mult = c(.8, .8))) +
  labs(title = "孝道責任感指標（原始）分布箱型圖",
       subtitle = "2021 臺灣社會變遷基本調查_家庭組",
       y = "") +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.06),
    plot.subtitle = element_text(hjust = -0.06))
plot_dvo

2.1.3 因素分析

2.1.3.1 平行因素分析 (Parallel Factor Analysis, PFA)

下圖是「孝道責任感指標（原始題組）」的平行因素分析陡坡圖。藍線代表實際資料的特徵值，紅線則是隨機模擬與重抽樣的結果。從圖中可以看到，前三個因素的藍線都高於紅線，而第四個因素之後則不再超過模擬資料，顯示此題組較適合區分為三個面向。

E1_pfa <- data %>% select(a:i) %>% as.data.frame()

fa.parallel(E1_pfa,
            fa = "fa",
            fm = "pa",
            cor = "poly")

#### 探索性因素分析 (Exploratory Factor Analysis, EFA)

在進行平行因素分析後，得知該題組有3個面向，故進探索性因素分析；結果顯示第一個因素（PA1）主要由題項 a、f、e、b、g 所構成，因素負荷量約介於 0.5 至 0.9 之間；第二個因素（PA2）主要包含題項 i 與 h，負荷量約為 0.7；第三個因素（PA3）則由題項 d 與 c 所組成，負荷量約介於 0.5 至 0.7。

E1_fa <- fa(
  data %>% select(a:i),
  nfactors = 3,
  rotate = "promax",
  fm = "pa",
  cor = "poly")

fa.diagram(E1_fa)

2.1.4 合併變項

因第一個面向的題數最多，故我們取第一面向的題目進行新的變項合併。 a. 對父母的養育之恩心存感激。 f. 父母去世，不管住得多遠，都親自奔喪。 e. 奉養父母使他們生活更舒適。 b. 無論父母對您如何不好，仍然善待他們。 g. 為了顧及父母的面子，為父母說些好話。

# 根據因素分析結果，建立新變項E1_m作為分析依變項，其中包含題項a, f, e, b, g
ata <- mutate(.data = data, E1_m = a + f + e + b + g)
data <- data %>%
mutate(E1_m = rowMeans(across(c(a, f, e, b, g)), na.rm = TRUE))

2.1.4.1 依變項（孝道責任感指標）變項敘述統計

根據敘述統計結果可得知，臺灣民眾孝道責任感（合併後）的分數平均為4.463，最小值為1.000，最大值為5.000，第一四分位數為4.200，中位數為4.600，第三四分位數為5.000（最大值一樣）；而遺漏值為2。

summary(data$E1_m)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   4.200   4.600   4.463   5.000   5.000       3

2.1.4.2 ’孝道責任感指標’變項統計圖表

由孝道責任感分布箱型圖可見，受訪者之孝道責任感指標整體分布偏高，中位數落在高分區間，顯示多數民眾對孝道責任具有高度認同。分布下緣出現少數低分離群值，反映部分受訪者對孝道責任之認同程度較低，但整體變異幅度有限。

plot_dv <- ggplot(data, aes(y = E1_m)) +
  geom_boxplot(fill = "#FFCC99", colour = "#CC9966", alpha = 1,
               width = .5, na.rm = TRUE) +
  scale_x_discrete(expand = expansion(mult = c(.8, .8))) +
  labs(title = "孝道責任感指標分布箱型圖",
       subtitle = "2021 臺灣社會變遷基本調查_家庭組",
       y = "") +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.06),
    plot.subtitle = element_text(hjust = -0.06))
plot_dv

2.2 自變項（性別、年齡、職業、教育、月收入）

2.2.1 性別

建立新變項’male’，並將原變項’a1’重新編碼為”男性=1，女性=0”（女性因樣本數較多，故作為對照組），由於此性別變項無遺漏值，故不進行清理遺漏值的動作。

2.2.1.1 性別：敘述統計

gender_summary <- data %>%
  filter(!is.na(male)) %>%
  count(male, name = "n") %>%
  mutate(百分比 = round(n / sum(n) * 100, 1))

gender_summaryT <- gender_summary %>%
  rename(性別 = male) %>%
  bind_rows(tibble(性別 = "總計", n = sum(gender_summary$n), 百分比 = 100))

gender_summaryT
## # A tibble: 3 × 3
##   性別      n 百分比
##   <chr> <int>  <dbl>
## 1 女性    842   53.5
## 2 男性    731   46.5
## 3 總計   1573  100

2.2.1.2 性別：統計圖表

根據性別長條圖顯示，樣本中女性受訪者人數略多於男性，整體性別分布大致均衡。

plot_male <- ggplot(gender_summary, aes(x = male, y = n, fill = male)) +
  geom_col(width = .5) +
  labs(
    x = "", y = "",
    title = "性別之次數分配長條圖",
    subtitle = "2021 臺灣社會變遷基本調查_家庭組"
  ) +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -.06),
    plot.subtitle = element_text(hjust = -0.08), 
    legend.position = "none")

plot_male

2.2.2 年齡

首先將原始資料 ‘a2y’和 ’year_m中的遺漏值清除（997、998），並且分別建立新變項’a2yn’和’year_mn’。最後將新變項調查年份（year_mn）-受訪者出生年（a2yn）＋1（虛歲），建立成新變項’年齡（age）’。

2.2.2.1 年齡：敘述統計

根據敘述統計可得知，年齡的平均值為50.19，中位數為49.00，最小值為20.00，最大值為95.00，第一四分位數為37.00，第三四分位數為63.00，遺漏值則為16個。

summary(data$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   20.00   37.00   49.00   50.19   63.00   95.00      16

2.2.2.2 年齡：統計圖表

根據年齡箱型圖顯示，樣本之年齡中位數約落在中壯年區間，年齡分布範圍涵蓋年輕至高齡族群，顯示樣本在年齡結構上具一定多樣性。

plot_age <- ggplot(data, aes(y = age)) +
  geom_boxplot(fill = "#FFCC99", colour = "#CC9966", alpha = 1,
               width = .5, na.rm = TRUE) +
  scale_x_discrete(expand = expansion(mult = c(.8, .8))) +
  labs(title = "年齡之箱型圖",       
       subtitle = "2021 臺灣社會變遷基本調查_家庭組", y = "") +
  theme(text = element_text(family = "PingFang TC"),
        plot.title = element_text(hjust = -0.04), 
        plot.subtitle = element_text(hjust = -0.06))
plot_age

2.2.3 教育年數

建立新變項’edu_y’，將原始變項’a10’分為以下年數：；[小學的教育年數設為6年]；；；；；；；；最後再將選項其他(22)做遺漏值處理。

2.2.3.1 教育年數：敘述統計

根據敘述統計可得知，教育年數的平均值為12.76，中位數為12.00，最小值為0.00，最大值為22.00，第一四分位數為12.00，第三四分位數為16.00，遺漏值則為5個。顯示大多數受訪者都具有高中以上學歷。

summary(data$edu_y) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   12.00   12.00   12.76   16.00   22.00       5

2.2.3.2 教育年數：統計圖表

根據教育年數箱型圖顯示，樣本之教育年數中位數落在約12年左右，多數受訪者集中於高中至專科教育程度，顯示樣本教育年數分布具一定集中趨勢，並伴隨少數低教育年數之離群值。

plot_eduy <- ggplot(data, aes(y = edu_y)) +
  geom_boxplot(fill = "#FFCC99", colour = "#CC9966", alpha = 1, width = .5, na.rm = TRUE) +
  scale_x_discrete(expand = expansion(mult = c(.8, .8))) +
  labs(title = "教育年數之箱型圖",
       subtitle = "2021 臺灣社會變遷基本調查_家庭組", y = " ") +
  theme(text = element_text(family = "PingFang TC"), plot.title = element_text(hjust = -0.06),
        plot.subtitle = element_text(hjust = -0.08))
plot_eduy

2.2.4 月收入

建立新變項inc，針對原始資料k14的遺漏值代碼97、98做遺漏值處理，其他皆維持原始資料的數值。接著建立新變項”inc_n”，以月收入原始變項(inc)的組別數字作為數值，以轉換為數值變項。

2.2.4.1 月收入：敘述統計

根據敘述統計可得知，月收入組距的平均值為5.457，中位數為5.00，最小值為1.00，最大值為23.00，第一四分位數為3.00，第三四分位數為7.00，遺漏值則為47個。

summary(data$inc_n) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   3.000   5.000   5.458   7.000  23.000      47

2.2.4.2 月收入：統計圖表

下圖是月收入類別變項的次數分配長條圖；由圖表可知月收入佔比最高的區間為30,001–40,000元；整體趨勢則顯示人數隨著月收入越高而有逐漸降低的趨勢（呈現右傾）。

plot_inc1 <- ggplot(data[!is.na(data$inc),] %>% mutate(
  inc_f = factor(inc, levels = 1:23, labels = c(
    "無收入", "<10,000", "10,001~20,000", "20,001~30,000", "30,001~40,000", 
    "40,001~50,000", "50,001~60,000", "60,001~70,000", "70,001~80,000", 
    "80,001~90,000", "90,001~100,000", "100,001~110,000", "110,001~120,000", 
    "120,001~130,000", "130,001~140,000", "140,001~150,000", "150,001~160,000",
    "160,001~170,000", "170,001~180,000", "180,001~190,000", "190,001~200,000",
    "200,001~300,000", ">300,000" ), ordered = TRUE)),
  aes(x = inc_f)) + geom_bar(fill = "turquoise3", color = "grey60") +
  labs(title = "個人月收入之次數分配長條圖",
       subtitle = "2021 臺灣社會變遷基本調查_家庭組", x = " ", y = " ") +
  guides(x = guide_axis(angle = 45)) +
  theme(text = element_text(family = "PingFang TC"),
        plot.title = element_text(hjust = -0.05),
        plot.subtitle = element_text(hjust = -0.06))
plot_inc1

下圖則為月收入組距的箱型圖，該變項是將月收入類別變項的數值轉為類別變項後的。從統計圖表可得知，中位數落在5的地方，對照原本月收入組別就是「30,001–40,000元」區間，所以仍可看出整體薪資收入偏向右傾。

plot_inc2 <- ggplot(data, aes(y = inc_n)) +
  geom_boxplot(fill = "#FFCC99", colour = "#CC9966", alpha = 1,
               width = .5, na.rm = TRUE) +
  scale_x_discrete(expand = expansion(mult = c(.8, .8))) +
  labs(title = "月收入組距之箱型圖",
       subtitle = "2021 臺灣社會變遷基本調查_家庭組", y = " ") +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.06),
    plot.subtitle = element_text(hjust = -0.08))

plot_inc2

2.2.5 職業

首先建立新變項’temp’，將原始職業代碼’k8b5’除以1000再取整數。再建立新變項’occ10’，最後依據聯合國國際職業標準分類（International Standard Classification of Occupations, ISCO）分為下列類別： 若temp值為1，歸類為1＝民意代表與高階管理人員； 若temp值為2，歸類為2＝專業人員； 若temp值為3，歸類為3＝技術員及助理專業人員； 若temp值為4，歸類為4＝事務支援人員； 若temp值為5，歸類為5＝服務及銷售工作人員； 若temp值為6，歸類為6＝農、林、漁、牧業生產人員； 若temp值為7，歸類為7＝技藝有關工作人員； 若temp值為8，歸類為8＝機械設備操作及組裝人員； 若temp值為9，歸類為9＝基層技術工及勞力工； 若temp值為0，歸類為10＝軍人。

2.2.5.1 職業：敘述統計

由敘述統計可得知，在職業當中佔比最多的樣本為「服務及銷售工作人員」（21.6%），樣本數為340，故後續回歸分析時以此選項作為對照組。

occ_summary <- data %>%
  filter(!is.na(occ10)) %>%
  count(occ10) %>%
  mutate(
    百分比 = round(n / sum(n) * 100, 1)  ) %>%
  bind_rows(
    summarise(., occ10 = "總計", n = sum(n),
              百分比 = 100.0))
occ_summary
## # A tibble: 11 × 3
##    occ10                        n 百分比
##    <chr>                    <int>  <dbl>
##  1 民意代表與高階管理人員      56    3.6
##  2 專業人員                   179   11.4
##  3 技術員及助理專業人員       254   16.1
##  4 事務支援人員               184   11.7
##  5 服務及銷售工作人員         340   21.6
##  6 農、林、漁、牧業生產人員    64    4.1
##  7 技藝有關工作人員           157   10  
##  8 機械設備操作及組裝人員     140    8.9
##  9 基層技術工及勞力工         191   12.1
## 10 軍人                         8    0.5
## 11 總計                      1573  100

2.2.5.2 職業：統計圖表

職業類別以橫向次數分配長條圖呈現，其中「服務及銷售工作人員」為最多數（佔比21.6%），最少的為「軍人」（佔比0.5%）。

plot_occ <- ggplot(occ_summary %>%
                  filter(occ10 != "總計") %>%
                  mutate(
                    occ10 = factor(occ10,
                                  levels = c("民意代表與高階管理人員",
                                             "專業人員", "技術員及助理專業人員",
                                             "事務支援人員", "服務及銷售工作人員",
                                             "農、林、漁、牧業生產人員",
                                             "技藝有關工作人員",
                                             "機械設備操作及組裝人員",
                                             "基層技術工及勞力工",
                                             "軍人")),occ10 = fct_rev(occ10)),
                aes(x = occ10, y = n)) +
  geom_col(fill = "#FF6699", col = "#FF99CC", alpha = 1, width = .5) +
  labs(x = "", y = "",
    title = "職業之次數分配長條圖", subtitle = "2021 臺灣社會變遷基本調查_家庭組") +
  theme(text = element_text(family = "PingFang TC"),
        plot.title = element_text(hjust = 0),
        plot.subtitle = element_text(hjust = 0),
        axis.text.x = element_text(angle = 0, hjust = .5))+coord_flip()
plot_occ

3 迴歸分析

3.1 單迴歸分析

我們首先進行了單變項線性迴歸分析，目的在於檢視各變項單獨納入模型時對依變項（E1_m）的解釋力。結果顯示，不同社會人口學變項在解釋依變項上具有明顯的差異，據此將作為後續多元迴歸模型之變項篩選依據。 - 建議納入最終模型的變項：

  • age（年齡）
    • F-statistic = 100.8（p < 2.2e-16）
    • Multiple R² = 0.0609
    • 對 E1_m 具有穩定且顯著的解釋力
  • edu_y（教育年數）
    • F-statistic = 91.02（p < 2.2e-16）
    • Multiple R² = 0.0550
    • 與 E1_m 呈現系統性關聯
  • occ10（職業類別）
    • F-statistic = 3.52（p = 0.00025）
    • Multiple R² = 0.0199
    • 顯示職業位置對 E1_m 具有結構性影響
    • 對照組設定為樣本數最多之職業類別
  • inc_n（收入）
    • F-statistic = 5.37（p = 0.0206）
    • Multiple R² = 0.0035
    • 解釋力有限，具控制變項價值

  • 不建議納入最終模型的變項

  • male（性別）

    • F-statistic = 0.62（p = 0.43）
    • Multiple R² = 0.0004
    • 對 E1_m 之解釋力有限
#確認過"occ10"的次數分配後，得知「服務及銷售工作人員」為做多人數的組別，因此將其指定為對照組。
data$occ10 <- relevel(
  factor(data$occ10),
  ref = "服務及銷售工作人員")

vars <- c("male", "age", "edu_y", "inc_n", "occ10")

for (v in vars) {
  cat("\n========================\n")
  cat("單回歸模型：E1_m ~", v, "\n")

  fml <- as.formula(paste("E1_m ~", v))
  mod <- lm(fml, data = data)
  print(summary(mod)) 
}
## 
## ========================
## 單回歸模型：E1_m ~ male 
## 
## Call:
## lm(formula = fml, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4505 -0.2740  0.1495  0.5260  0.5495 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.47405    0.02034  219.97   <2e-16 ***
## male男性    -0.02357    0.02983   -0.79     0.43    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5895 on 1568 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.000398,   Adjusted R-squared:  -0.0002395 
## F-statistic: 0.6243 on 1 and 1568 DF,  p-value: 0.4296
## 
## 
## ========================
## 單回歸模型：E1_m ~ age 
## 
## Call:
## lm(formula = fml, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3634 -0.2586  0.1587  0.3943  0.8010 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.0258835  0.0456089   88.27   <2e-16 ***
## age         0.0086541  0.0008621   10.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5723 on 1553 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.06093,    Adjusted R-squared:  0.06033 
## F-statistic: 100.8 on 1 and 1553 DF,  p-value: < 2.2e-16
## 
## 
## ========================
## 單回歸模型：E1_m ~ edu_y 
## 
## Call:
## lm(formula = fml, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5916 -0.2838  0.1110  0.4084  0.8529 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.899389   0.048032 102.003   <2e-16 ***
## edu_y       -0.034197   0.003584  -9.541   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5737 on 1563 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.05503,    Adjusted R-squared:  0.05443 
## F-statistic: 91.02 on 1 and 1563 DF,  p-value: < 2.2e-16
## 
## 
## ========================
## 單回歸模型：E1_m ~ inc_n 
## 
## Call:
## lm(formula = fml, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4638 -0.2638  0.1462  0.5060  0.7171 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.514078   0.028103 160.626   <2e-16 ***
## inc_n       -0.010051   0.004338  -2.317   0.0206 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5902 on 1522 degrees of freedom
##   (49 observations deleted due to missingness)
## Multiple R-squared:  0.003515,   Adjusted R-squared:  0.002861 
## F-statistic: 5.369 on 1 and 1522 DF,  p-value: 0.02063
## 
## 
## ========================
## 單回歸模型：E1_m ~ occ10 
## 
## Call:
## lm(formula = fml, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5080 -0.2832  0.1538  0.4618  0.6618 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.514012   0.031783 142.025  < 2e-16 ***
## occ10民意代表與高階管理人員    0.060095   0.084412   0.712 0.476614    
## occ10專業人員                 -0.175810   0.054167  -3.246 0.001196 ** 
## occ10技術員及助理專業人員     -0.161256   0.048563  -3.321 0.000919 ***
## occ10事務支援人員             -0.067816   0.053585  -1.266 0.205849    
## occ10農、林、漁、牧業生產人員  0.154217   0.079756   1.934 0.053340 .  
## occ10技藝有關工作人員         -0.030785   0.056492  -0.545 0.585876    
## occ10機械設備操作及組裝人員   -0.012583   0.058790  -0.214 0.830546    
## occ10基層技術工及勞力工       -0.006029   0.053033  -0.114 0.909499    
## occ10軍人                     -0.089012   0.209324  -0.425 0.670724    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5852 on 1560 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.0199, Adjusted R-squared:  0.01425 
## F-statistic:  3.52 on 9 and 1560 DF,  p-value: 0.0002499

3.1.1 雙變項圖表：孝道責任感與性別

孝道責任感與性別箱型圖顯示不同性別在孝道責任感指標上的分布情形。由箱型圖可見，男性與女性的中位數與四分位距高度重疊，整體分布相近，僅在低分數處出現少數離群值。整體而言，性別在孝道責任感上的差異不明顯，顯示性別與孝道責任感之間並無顯著關聯。

ggplot(data, aes(x = male, y = E1_m, fill = male)) +
  geom_boxplot(alpha = 1, width = 0.5, na.rm = TRUE) +
  labs(
    title = "孝道責任感指標與性別之箱型圖",
    subtitle = "2021 臺灣社會變遷基本調查_家庭組",
    x = " ",
    y = " ",
    fill = "性別") +
  theme(legend.position = "none") +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.06),
    plot.subtitle = element_text(hjust = -0.06))

3.2 初步多元迴歸模型

在完成單變項迴歸分析並篩選出具解釋潛力之變項後，我們進一步建構初步的多元迴歸模型，同時納入了年齡、教育年數、收入及職業類別，檢視各變項在相互控制下對依變項（孝道責任感, E1_m） 的影響。此一步驟旨在評估各社會人口與社經變項於多元情境中的相對效果，並避免單變項分析可能產生的混淆影響。

分析結果顯示，整體模型達統計顯著（F=10.68，p<.001），顯示所納入之變項能共同解釋依變項（孝道責任感, E1_m）的變異。於控制其他變項後，年齡仍與依變項（孝道責任感, E1_m）呈現顯著正向關聯，而教育年數則呈現顯著負向關聯，顯示世代位置與教育背景對該構面具有穩定影響。相較之下，收入在多元模型中未達顯著，顯示其於單變項分析中所呈現的效果，可能部分來自與年齡及教育之重疊影響。

在職業類別方面，整體職業結構對依變項（孝道責任感, E1_m）仍具有解釋價值，但多數職業類別在控制年齡與教育後，其差異程度有所減弱，僅部分職業類別仍與對照組呈現顯著差異。此結果顯示，職業位置對依變項（孝道責任感, E1_m）的影響，部分係透過教育與年齡結構間接反映，而非完全獨立作用。

data$occ10 <- relevel(
  factor(data$occ10),
  ref = "服務及銷售工作人員")

model <- lm(E1_m ~ age + edu_y + inc_n + occ10, data = data)
summary(model)
## 
## Call:
## lm(formula = E1_m ~ age + edu_y + inc_n + occ10, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4515 -0.2467  0.1545  0.3797  0.8691 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.461212   0.113928  39.158  < 2e-16 ***
## age                            0.005845   0.001127   5.185 2.45e-07 ***
## edu_y                         -0.021356   0.005608  -3.808 0.000146 ***
## inc_n                          0.005238   0.004863   1.077 0.281547    
## occ10民意代表與高階管理人員    0.050143   0.089224   0.562 0.574208    
## occ10專業人員                 -0.084158   0.058398  -1.441 0.149761    
## occ10技術員及助理專業人員     -0.116007   0.050293  -2.307 0.021213 *  
## occ10事務支援人員             -0.016430   0.054090  -0.304 0.761358    
## occ10農、林、漁、牧業生產人員 -0.024311   0.084295  -0.288 0.773080    
## occ10技藝有關工作人員         -0.071683   0.056664  -1.265 0.206045    
## occ10機械設備操作及組裝人員   -0.037878   0.058760  -0.645 0.519266    
## occ10基層技術工及勞力工       -0.085905   0.054911  -1.564 0.117928    
## occ10軍人                     -0.059286   0.204587  -0.290 0.772021    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5706 on 1495 degrees of freedom
##   (65 observations deleted due to missingness)
## Multiple R-squared:  0.07899,    Adjusted R-squared:  0.0716 
## F-statistic: 10.68 on 12 and 1495 DF,  p-value: < 2.2e-16
stargazer(model,
  type = "text",
  title = "孝道責任感之多元線性迴歸分析",
  dep.var.labels = "孝道責任感指標",
  no.space = TRUE,
  digits = 3)
## 
## 孝道責任感之多元線性迴歸分析
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               孝道責任感指標          
## -----------------------------------------------
## age                          0.006***          
##                               (0.001)          
## edu_y                        -0.021***         
##                               (0.006)          
## inc_n                          0.005           
##                               (0.005)          
## occ10民意代表與高階管理人員               0.050           
##                               (0.089)          
## occ10專業人員                     -0.084           
##                               (0.058)          
## occ10技術員及助理專業人員              -0.116**          
##                               (0.050)          
## occ10事務支援人員                   -0.016           
##                               (0.054)          
## occ10農、林、漁、牧業生產人員             -0.024           
##                               (0.084)          
## occ10技藝有關工作人員                 -0.072           
##                               (0.057)          
## occ10機械設備操作及組裝人員              -0.038           
##                               (0.059)          
## occ10基層技術工及勞力工                -0.086           
##                               (0.055)          
## occ10軍人                       -0.059           
##                               (0.205)          
## Constant                     4.461***          
##                               (0.114)          
## -----------------------------------------------
## Observations                   1,508           
## R2                             0.079           
## Adjusted R2                    0.072           
## Residual Std. Error      0.571 (df = 1495)     
## F Statistic          10.685*** (df = 12; 1495) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

3.3 巢狀模型比較（nested model ANOVA）

為檢驗新增自變項是否具有額外解釋力，我們進行了巢狀模型比較（nested model ANOVA）。結果顯示，在控制年齡與教育年數後，加入收入變項（inc_n）並未使模型配適度顯著提升（F=0.94，p=0.334）。進一步於模型中納入職業類別（occ10）後，整體模型解釋力亦未出現顯著改善（F=1.08，p=0.373）。

此結果說明，雖然在完整模型中部分職業類別之迴歸係數達顯著水準，但整體而言，收入與職業類別作為變項組別，並未對模型解釋力帶來額外貢獻。相較之下，年齡與教育年數為本研究中較為關鍵的解釋變項。

anova(model)
## Analysis of Variance Table
## 
## Response: E1_m
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## age          1  31.83  31.826 97.7420 < 2.2e-16 ***
## edu_y        1   6.45   6.447 19.7995 9.237e-06 ***
## inc_n        1   0.30   0.305  0.9355    0.3336    
## occ10        9   3.17   0.352  1.0819    0.3729    
## Residuals 1495 486.79   0.326                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.3.1 雙變項圖表：孝道責任感與月收入

下圖為孝道責任感指標與各月收入級距之散佈情形。整體來看，各收入級距中的孝道責任感分數分布相當分散，未呈現隨收入提高而明顯上升或下降的趨勢。無論低收入或高收入族群，皆可觀察到高低不一的孝道責任感分數，顯示月收入與孝道責任感之間並未呈現明確關聯。

inc_labs <- c(
  "無收入","<10,000","10,001~20,000","20,001~30,000","30,001~40,000",
  "40,001~50,000","50,001~60,000","60,001~70,000","70,001~80,000",
  "80,001~90,000","90,001~100,000","100,001~110,000","110,001~120,000",
  "120,001~130,000","130,001~140,000","140,001~150,000","150,001~160,000",
  "160,001~170,000","170,001~180,000","180,001~190,000","190,001~200,000",
  "200,001~300,000",">300,000")
data <- data %>%
  mutate(
    inc_f = factor(inc, levels = 1:23, labels = inc_labs, ordered = TRUE),
    inc_f = fct_explicit_na(inc_f, na_level = "NA"))
pal <- c(rep("darkorange3", 23),"grey80")
names(pal) <- levels(data$inc_f)
ggplot(data, aes(x = inc_f, y = E1_m, color = inc_f)) +
  geom_jitter(width = 0.25, alpha = 0.6, size = 1.8, stroke = 0) +
  scale_color_manual(values = pal, drop = FALSE) +
  labs(
    title = "孝道責任感與各月收入級距之散布圖",
    subtitle = "2021 臺灣社會變遷基本調查_家庭組", x = " ", y = " ") +
  guides(x = guide_axis(angle = 45)) +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.07, face = "bold"),
    plot.subtitle = element_text(hjust = -0.06),
    legend.position = "none",
    axis.text.x = element_text(size = 8))

### 雙變項圖表：孝道責任感與職業 孝道責任感與職業箱型圖顯示不同職業類別在孝道責任感指標上的分布情形。整體而言，各職業類別的中位數多集中於偏高分區間，且箱體之間高度重疊，顯示不同職業間的孝道責任感差異並不明顯。雖然部分職業類別在分布範圍與離群值上略有差異，但整體未呈現明確的職業分化趨勢，顯示職業類別與孝道責任感之間的關聯程度有限。

pal_occ10 <- c(
  "民意代表與高階管理人員" = "#E15759",  "專業人員" = "#377EB8",
  "技術員及助理專業人員" = "#59A14F",  "事務支援人員" = "#9C755F",
  "服務及銷售工作人員" = "#F28E2B",  "農、林、漁、牧業生產人員" = "#EDC948",
  "技藝有關工作人員" = "#76B7B2",  "機械設備操作及組裝人員" = "#AF7AA1",
  "基層技術工及勞力工" = "#00BCD4",  "軍人" = "#FF9DA7")

ggplot(data, aes(x = fct_rev(occ10), y = E1_m, fill = occ10)) +
  geom_boxplot(na.rm = TRUE, alpha = 1, color = "grey40") +
  scale_fill_manual(values = pal_occ10) +
  labs(title = "孝道責任感與各職業類別之箱型圖",
       subtitle = "2021 臺灣社會變遷基本調查_家庭組",
       x = " ", y = " " ) +
  theme(
    legend.position = "none",
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.07, face = "bold", size = 14),
    plot.subtitle = element_text(hjust = -0.06, size = 12),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 1, margin = margin(t = 6)),
    plot.margin = margin(t = 10, r = 10, b = 0, l = 25))+coord_flip()

## 最終回歸模型 綜上分析，我們發現年齡與教育年數為影響E1_m的主要解釋變項。相較之下，收入與職業類別在控制年齡與教育年數後，未能提供顯著的額外解釋力。其最終迴歸方程式如下： \[\text{孝道責任感指標} = \beta_0 + \beta_1(\text{年齡}) + \beta_2(\text{教育年數}) + \varepsilon\]

final_model <- lm(E1_m ~ age + edu_y, data = data)
summary(final_model)
## 
## Call:
## lm(formula = E1_m ~ age + edu_y, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4667 -0.2525  0.1537  0.3917  0.7838 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.407995   0.098157  44.908  < 2e-16 ***
## age          0.006004   0.001051   5.711 1.34e-08 ***
## edu_y       -0.019490   0.004414  -4.416 1.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5691 on 1549 degrees of freedom
##   (21 observations deleted due to missingness)
## Multiple R-squared:  0.07296,    Adjusted R-squared:  0.07176 
## F-statistic: 60.95 on 2 and 1549 DF,  p-value: < 2.2e-16
stargazer(
  final_model,
  type = "text",
  title = "孝道責任感之多元線性迴歸分析",
  dep.var.labels = "孝道責任感指標",
  digits = 3)
## 
## 孝道責任感之多元線性迴歸分析
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               孝道責任感指標          
## -----------------------------------------------
## age                          0.006***          
##                               (0.001)          
##                                                
## edu_y                        -0.019***         
##                               (0.004)          
##                                                
## Constant                     4.408***          
##                               (0.098)          
##                                                
## -----------------------------------------------
## Observations                   1,552           
## R2                             0.073           
## Adjusted R2                    0.072           
## Residual Std. Error      0.569 (df = 1549)     
## F Statistic          60.953*** (df = 2; 1549)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

3.3.2 雙變項圖表：孝道責任感與年齡

圖表顯示孝道責任感指標與年齡之間的關係。從散佈情形與趨勢線可見，隨著年齡增加，孝道責任感呈現緩慢上升的趨勢。配合後續統計檢定結果可知，年齡與教育年數皆與孝道責任感達統計顯著關聯，顯示可能基於不同生命歷程與教育背景，會影響個體對孝道責任的認知與態度。

ggplot(data, aes(x = age, y = E1_m)) +
  geom_point(alpha = 0.5, color = "#66CCCC", size = 2, na.rm = TRUE) +
  geom_smooth(method = "lm", color = "red", se = TRUE, na.rm = TRUE, size = 1.4) +
  scale_x_continuous(
    breaks = seq(20, 100, by = 10),
    labels = paste0(seq(20, 100, by = 10), "歲")) +
  labs(
    title = "孝道責任感與年齡之散佈圖",
    subtitle = "2021 臺灣社會變遷基本調查_家庭組",
    x = " ",
    y = " ") +
  theme(
    text = element_text(family = "PingFang TC"),
    plot.title = element_text(hjust = -0.01, size = 14),
    plot.subtitle = element_text(hjust = -0.01))

### 雙變項圖表：孝道責任感與教育年數 該圖表顯示孝道責任感指標在不同教育程度下的分布情形。雖然教育年數本質上屬於連續變項，但為了便於比較各教育層級之整體分布情形，本研究仍以教育程度分組（原類別變項），並以箱型圖呈現。從圖中可見，不同教育程度之間的孝道責任感分布存在差異，且配合統計檢定結果顯示，教育年數與孝道責任感之間具有顯著關聯。整體而言，各教育層級內部仍可觀察到一定程度的分散，顯示個體差異仍然存在。

edu_levels <- c(0, 6, 9, 12, 14, 15, 16, 18, 22)
edu_labels <- c(
  "未受教育","小學","國中","高中（職）","專科","三專","大學","碩士","博士")

ggplot(data %>% dplyr::filter(!is.na(edu_y))
       ,aes(x = factor(edu_y,
                       levels = edu_levels,
                       labels = edu_labels),
            y = E1_m, fill = factor(edu_y))) +
  geom_boxplot(alpha = 1) +scale_fill_brewer(
    palette = "Blues") +
  labs(
    title = "孝道責任感與教育程度分布箱型圖",
    subtitle = "2021 臺灣社會變遷基本調查＿家庭組",
    x = " ", y = "孝道責任感") +
  theme(
    text = element_text(family = "PingFang TC"),
    legend.position = "none")

### 雙變項圖表：年齡與教育年數 下圖顯示不同教育程度族群的年齡分布情形。可以看出，教育程度與年齡之間呈現明顯差異，整體而言，教育程度較低者年齡偏高，而教育程度較高者則集中於較年輕的年齡區間。部分教育層級（如國中、高中／高職）內部的年齡分布較為分散，顯示不同世代同時存在；相較之下，大學以上教育程度者的年齡分布則較為集中，反映教育擴張與世代差異的影響。

ggplot(
  data %>% dplyr::filter(!is.na(edu_y)),
  aes(x = factor(edu_y,
                 levels = c(0, 6, 9, 12, 14, 15, 16, 18, 22),
                 labels = c("未受教育", "小學", "國中",
                            "高中／高職", "專科", "三專",
                            "大學", "碩士", "博士")),
      y = age, fill = factor(edu_y))) +
  geom_boxplot() +
  scale_fill_brewer(palette = "Blues") +
  labs(title = "教育程度與年齡分布箱型圖",
       subtitle = "2021 臺灣社會變遷基本調查＿家庭組",
       x = " ", y = "年齡") +
  theme(
    text = element_text(family = "PingFang TC"),
    legend.position = "none")

## 迴歸模型診斷圖 由迴歸診斷圖可知，殘差與預測值之間未呈現明顯非線性關係，顯示線性假設大致成立。Q–Q圖顯示殘差在分布尾端略有偏離常態分布，惟在樣本數較大的情況下，此偏離對迴歸估計影響有限。Scale-Location圖顯示殘差變異隨預測值略有變化，可能存在輕微異質變異，但整體程度不嚴重。此外，Residuals vs Leverage圖未顯示具有高度影響力之觀測值。整體而言，本研究之最終迴歸模型符合線性迴歸之基本假設，分析結果具有解釋上的穩定性。

par(mfrow = c(2, 2), pty = "s")
plot(final_model)

3.4 共線性檢查

我們以VIF檢驗共線性，結果顯示年齡與教育年數之VIF值皆為1.50，遠低於一般常用判準（VIF < 5），著說明模型不存在顯著的共線性問題。

vif(final_model)
##      age    edu_y 
## 1.500544 1.500544

3.5 迴歸模型圖：估計值與信賴區間

下圖呈現孝道責任感之線性迴歸模型估計結果與其信賴區間。結果顯示，在控制其他變項後，年齡對孝道責任感具有顯著正向影響，而教育年數則呈現顯著負向影響。從圖中可見，兩項變數的估計值皆未跨越零值虛線，顯示其效果具有統計顯著性。整體而言，年齡較高者傾向具有較高的孝道責任感，而教育年數較高者則相對較低。

# 取出估計值與信賴區間
out_conf <- broom::tidy(final_model, conf.int = TRUE)
# 移除截距
out_conf <- subset(out_conf, term != "(Intercept)")
out_conf$new_term <- dplyr::recode(
out_conf$term, age = "年齡", edu_y = "教育年數")
out_conf$p_sig <- ifelse(out_conf$p.value < 0.05, "p < .05", "p ≥ .05")

ggplot(out_conf,
       aes(x = estimate, y = reorder(new_term, estimate),
           xmin = conf.low, xmax = conf.high, color = p_sig)) +
  geom_pointrange(size = .5) +
  geom_errorbar(width = .5) +
  geom_vline(xintercept = 0, linetype = 4) +
  scale_color_manual(values = c("p < .05" = "tomato2",
                                "p ≥ .05" = "grey30")) +
  labs(
    x = "線性迴歸估計值", y = NULL,
    title = "孝道程度之線性迴歸模型：估計值及其信賴區間",
    subtitle = "2021 臺灣社會變遷基本調查_家庭組",
    color = "") +
  geom_pointrange(size = .8) +
  geom_errorbar(width = .4) +
  theme(legend.position = "right")

3.5.1 多變項圖表：孝道責任感與年齡、教育年數

下圖為孝道責任感、年齡與教育程度之關係分面圖，依不同教育程度分組呈現年齡與孝道責任感之關聯。整體來看，多數教育層級中，孝道責任感皆隨年齡增加而呈現緩慢上升的趨勢，但各教育層級內的斜率與分布程度略有差異。教育程度較低的族群，其孝道責任感整體水準較高且分布較集中；相較之下，高等教育族群內部差異較大，但仍可觀察到年齡與孝道責任感之正向關聯。此結果顯示，年齡對孝道責任感的影響在不同教育背景下大致一致，但教育程度仍會調整其表現型態。

edu_levels <- c(0, 6, 9, 12, 14, 15, 16, 18, 22)
edu_labels <- c("未受教育", "小學", "國中", "高中／高職", 
                "專科", "三專", "大學", "碩士", "博士")

plot_dat <- data %>%
  filter(!is.na(E1_m), !is.na(age), !is.na(edu_y)) %>%
  mutate(edu_y_f = factor(edu_y, levels = edu_levels, labels = edu_labels))

ggplot(plot_dat, aes(x = age, y = E1_m)) + geom_point(alpha = 0.25, size = 1.5) +
  geom_smooth(method = "lm", se = TRUE, color = "red", linewidth = 1.1) +
  facet_wrap(~ edu_y_f, ncol = 3) +
  scale_x_continuous(breaks = seq(20, 90, by = 20),
                     labels = paste0(seq(20, 90, by = 20), "歲")) +
  scale_y_continuous(breaks = 1:5) + coord_cartesian(ylim = c(1, 5)) +
  labs(title = "孝道責任感、年齡與教育程度之關係（分面圖）",
       subtitle = "2021 臺灣社會變遷基本調查＿家庭組",
       x = "", y = "") + theme(text = element_text(family = "PingFang TC"))

3.5.2 多變項圖表：孝道責任感與性別、教育年數

下圖顯示在不同教育程度下，男性與女性於孝道責任感指標上的分布情形。從各教育層級的箱型圖可見，孝道責任感隨教育程度而呈現差異，顯示教育背景對孝道責任感具有一定影響；相較之下，各教育層級內男性與女性的中位數與分布高度重疊，未呈現明顯差異，顯示性別本身並非影響孝道責任感的主要因素。

edu_levels <- c(0, 6, 9, 12, 14, 15, 16, 18, 22)

edu_labels <- c("未受教育", "小學", "國中",
                "高中／高職", "專科", "三專",
                "大學", "碩士", "博士")

plot_dat <- data %>%
  filter(!is.na(E1_m), !is.na(male), !is.na(edu_y)) %>%
  mutate(edu_y_f = factor(edu_y,
                          levels = edu_levels,
                          labels = edu_labels))
ggplot(plot_dat,
       aes(x = edu_y_f, y = E1_m, fill = male)) +
  geom_boxplot(
    position = position_dodge(width = 0.8),
    outlier.alpha = 0.3) +
  labs(title = "孝道責任感、性別與教育程度之關係（箱型圖）",
       subtitle = "2021 臺灣社會變遷基本調查＿家庭組",
      x = "", y = " ", fill = "") +
  theme(text = element_text(family = "PingFang TC")) +
  coord_flip()

3.5.3 多變項圖表：孝道責任感與年齡、職業

該圖表為孝道責任感、年齡與職業類別之關係分面圖。整體來看，多數職業類別中，孝道責任感皆隨年齡增加而呈現緩慢上升的趨勢，顯示年齡與孝道責任感之間具有一致的正向關聯。相較之下，不同職業類別之間的分布型態與趨勢線大致相近，未呈現明顯的職業分化，顯示職業類別本身對孝道責任感的影響有限，主要仍反映年齡效果。

ggplot(data, aes(x = age, y = E1_m)) +
  geom_point(alpha = 0.5, colour = "#6a8caf", size = 1.3) +
  geom_smooth(
    method = "lm", se = TRUE, colour = "firebrick", linewidth = 1) +
  facet_wrap(~ occ10, nrow = 2, ncol = 5) +
  scale_x_continuous(
    breaks = seq(20, 90, by = 20),
    labels = paste0(seq(20, 90, by = 20), "歲")) +
  labs(
    title = "孝道責任感、年齡與職業類別之關係（分面圖）",
    subtitle = "2021 臺灣社會變遷基本調查_家庭組",
    x = " ", y = "孝道責任感指標") +
  theme(text = element_text(family = "PingFang TC"))

# Appendix

# 安裝套件與讀取資料
library(knitr)
library(showtext)
library(dplyr)
library(magrittr)
library(haven)
library(stargazer)
library(forcats)
library(pscl)
library(lavaan)
library(psych)
library(ggplot2)
library(scales)
library(GGally)
library(car)
library(stargazer)
library(broom)
library(stargazer)

showtext_auto()
par(family = "HeitiTC")
par(family = "PingFang TC")

data <- read_dta(file = "tscs211.dta", encoding = "utf-8")
#依變項："孝道責任感指標"資料清理
data <- mutate(.data = data,
               a = case_when(
                 e1a %in% c(93, 97, 98) ~ NA,
                 e1a == 0 ~ 1,
                 e1a == 1 ~ 2,
                 e1a == 2 ~ 3,
                 e1a == 3 ~ 4,
                 e1a == 4 ~ 5 ),
               b = case_when(
                 e1b %in% c(93, 97, 98) ~ NA,
                 e1b == 0 ~ 1,
                 e1b == 1 ~ 2,
                 e1b == 2 ~ 3,
                 e1b == 3 ~ 4,
                 e1b == 4 ~ 5 ),
               c = case_when(
                 e1c %in% c(93, 97, 98) ~ NA,
                 e1c == 0 ~ 1,
                 e1c == 1 ~ 2,
                 e1c == 2 ~ 3,
                 e1c == 3 ~ 4,
                 e1c == 4 ~ 5 ),
               d = case_when(
                 e1d %in% c(93, 97, 98) ~ NA,
                 e1d == 0 ~ 1,
                 e1d == 1 ~ 2,
                 e1d == 2 ~ 3,
                 e1d == 3 ~ 4,
                 e1d == 4 ~ 5 ),
               e = case_when(
                 e1e %in% c(93, 97, 98) ~ NA,
                 e1e == 0 ~ 1,
                 e1e == 1 ~ 2,
                 e1e == 2 ~ 3,
                 e1e == 3 ~ 4,
                 e1e == 4 ~ 5 ),
               f = case_when(
                 e1f %in% c(93, 97, 98) ~ NA,
                 e1f == 0 ~ 1,
                 e1f == 1 ~ 2,
                 e1f == 2 ~ 3,
                 e1f == 3 ~ 4,
                 e1f == 4 ~ 5 ),
               g = case_when(
                 e1g %in% c(93, 97, 98) ~ NA,
                 e1g == 0 ~ 1,
                 e1g == 1 ~ 2,
                 e1g == 2 ~ 3,
                 e1g == 3 ~ 4,
                 e1g == 4 ~ 5 ),
               h = case_when(
                 e1h %in% c(93, 97, 98) ~ NA,
                 e1h == 0 ~ 1,
                 e1h == 1 ~ 2,
                 e1h == 2 ~ 3,
                 e1h == 3 ~ 4,
                 e1h == 4 ~ 5 ),
               i = case_when(
                 e1i %in% c(93, 97, 98) ~ NA,
                 e1i == 0 ~ 1,
                 e1i == 1 ~ 2,
                 e1i == 2 ~ 3,
                 e1i == 3 ~ 4,
                 e1i == 4 ~ 5 ))

# 合併為孝道責任感指標分數（原始）
items <- c("a","b","c","d","e","f","g","h","i")
ata <- mutate(.data = data, E1_ori = a + b + c + d + e + f + g + h + i)
data <- data %>%
mutate(E1_ori = rowMeans(across(c(a, b, c, d, e, f, g, h, i)), na.rm = TRUE))
#自變項："性別"資料清理
data <- mutate(.data = data, 
               male = factor(ifelse(a1 == 1, 1, 0),
                  levels = c(0, 1),
                  labels = c("女性", "男性")))
#自變項："年齡"資料清理
data <- mutate(.data = data, 
                a2yn = ifelse(a2y %in% c(997, 998), NA, a2y),
               year_mn = ifelse(year_m %in% c(997, 998), NA, year_m),
               age = year_mn - a2yn + 1)
#自變項："教育年數"資料清理
data <- mutate(.data = data, 
                edu_y = case_when(
                 a10 == 1 ~ 0, a10 == 2 ~ 0, a10 == 3 ~ 6, a10 == 4 ~ 9, a10 == 5 ~ 9,   
                 a10 %in% c(6,7,8,9,13) ~ 12, a10 %in% c(10,11,14,15) ~ 14, a10 == 12 ~ 15, 
                 a10 %in% c(16,17,18,19) ~ 16, a10 == 20 ~ 18, a10 == 21 ~ 22, a10 == 22 ~ NA ))
#自變項："月收入"資料清理
data <- mutate(.data = data, 
                inc = case_when(
                 k14 %in% c(97, 98) ~ NA,
                 k14 == 1 ~ 1,
                 k14 == 2 ~ 2,
                 k14 == 3 ~ 3,
                 k14 == 4 ~ 4,
                 k14 == 5 ~ 5,
                 k14 == 6 ~ 6,
                 k14 == 7 ~ 7,
                 k14 == 8 ~ 8,
                 k14 == 9 ~ 9,
                 k14 == 10 ~ 10,
                 k14 == 11 ~ 11,
                 k14 == 12 ~ 12,
                 k14 == 13 ~ 13,
                 k14 == 14 ~ 14,
                 k14 == 15 ~ 15,
                 k14 == 16 ~ 16,
                 k14 == 17 ~ 17,
                 k14 == 18 ~ 18,
                 k14 == 19 ~ 19,
                 k14 == 20 ~ 20,
                 k14 == 21 ~ 21,
                 k14 == 22 ~ 22,
                 k14 == 23 ~ 23 ))

data <- data %>%
  mutate(
    inc_n = ifelse(k14 %in% c(97, 98), NA, as.numeric(as.character(k14))))
#自變項："職業"資料清理
data <- data %>%
  mutate(
    temp = trunc(k8b5 / 1000),
    occ10 = case_when(
      temp == 1 ~ 1L,      temp == 2 ~ 2L,
      temp == 3 ~ 3L,      temp == 4 ~ 4L,
      temp == 5 ~ 5L,      temp == 6 ~ 6L,
      temp == 7 ~ 7L,      temp == 8 ~ 8L,
      temp == 9 ~ 9L,      temp == 0 ~ 10L,
      TRUE ~ NA
    ),
    occ10 = factor(occ10, levels = 1:10, labels = c(
      "民意代表與高階管理人員", "專業人員", "技術員及助理專業人員", "事務支援人員",
      "服務及銷售工作人員", "農、林、漁、牧業生產人員", "技藝有關工作人員",
      "機械設備操作及組裝人員", "基層技術工及勞力工", "軍人")))