Tasks

這份作業分成三部分,分別是Qualtrics問卷的編輯、以R語言執行投票預測、以及將你的期末分組問卷報告撰寫成一個不超過兩段的摘要。請依照以下指示進行。

Question 1. Qualtrics問卷編輯

首先,將名為Assignment 3 Wave 1的qsf檔import進一個新的Qualtrics問卷檔中。你會看到已經替你完成的問卷格式。在此基礎上,你需要執行下列幾個步驟。

Q1a. 將共享資料匣中的名為location_DD的excel表單import進下圖的Drill-down question。並在Group 1及Group 2依序填入 縣市 及 行政區。按Preview預覽。

0.5分

Q1b. 前往Survey flow,在第一個question block下方按Add below。首先選擇Randomzier (淡紫色),接著選Group (淺藍色)。第一個group預定作為control,第二個group預定作為treatment。兩個groups下方分別依序加上Embedded data (淡綠色)和Block (灰色)。第一個group的Embedded data點下去時在Create New Field or Choose From Dropdown中鍵入treat,而一旁的Set a value now中則鍵入0;第二個group (treatment group)則在Create New Field or Choose From Dropdown中鍵入treat,而一旁的Set a value now中則鍵入1。如此,當問卷執行完後,你將有一個額外的欄位treat,標註填答者是被分入哪一個組別,便於辨識及資料整理。最後,在兩的groups各自的的Block處依序放入Control和Treatment block。完成後的樣貌會像下面這樣。

0.5分

Q1c. 前往dashboard的quotas區,按照下列圖示,依序設定quota blocks。總共需要設定六個quota blocks,題目的選擇來自第一個question block中的性別、年齡層。六個quota blocks的組名、性別年齡組合、以及佔總樣本(N = 1000)比率依序是:

男1: 男 x 18-29歲 (30%)
n: 1000 x 30% = 300
男2: 男 x 30-59歲 (15%)
n: 1000 x 15% = 150
男3: 男 x >60歲 (10%)
n: 1000 x 10% = 100
女1: 女 x 18-29歲 (15%)
n: 1000 x 15% = 150
女2: 女 x 30-59歲 (20%)
n: 1000 x 20% = 200
女3: 女 x >60歲 (10%)
n: 1000 x 10% = 100

請按下圖指示操作,不要忘記在右上角分子欄位設定各quota block的樣本數,格式為: 該block的配額為 0 / n 。 舉例說明,以第一組(男1)為例,鍵入的數字為 0 / 300。

同時,每一組quota block也請設定額滿登出的步驟,如下圖所示。

實際操作上,一個quota block額滿後,新進入的填答者即使滿足該block的條件(例如女 x 30-59歲),也無法再填寫,會被強制登出。

0.5分

Q1d. 前往Survey options處,依序如下圖操作,勾選必要之欄位。

Options \(\rightarrow\) Responses

Options \(\rightarrow\) Security 1

同一處再往下拉

Options \(\rightarrow\) Security 2

0.5分

Q1e. 產生虛擬回應(generate test responses)

請按照下圖指示生成1000則虛擬回應。

生成的速度與問卷長度和連線品質有關,請耐心等待^^;,待視窗中的灰色圈圈不再轉動後才按確認,關閉視窗(否則將無法確實產生你要求的筆數)。若因連線問題或其他因素導致未能在某一次嘗試中生成期望的回應數,則可在隔一陣子後重新產生一次(再按一次generate),補上額外的回應數量。若幾次下來總回應數 > 1000,則可在匯出資料後自行將尾端溢出的樣本列位(row)刪除。

完成後點擊進入Quotas版面,把各quota blocks的額滿進度秀給自己(以及我跟助教)看一下。顯示之數額超出配額沒有關係,就僅是資格符合被指派進來,但他們的回應不會被回收(你仍可以選擇在選單中勾選keep)。

0.5分

Q1f. 前往Data & Analysis版面查看,準備匯出資料。

選擇Excel,並分別以Export Values和Export labels兩種形式輸出。labels是指 問卷原始選項的文字內容,供你recode時參考。values則是Qualtrics自動將選項文字內容轉化為數值,例如男(1)、女(2),這種就是很標準的factor levels (而非dummy coding (0, 1))而非dummy coding (0, 1))。但有些使用者希望自行決定數值的順序或編碼規範,故不選用Export Values。又因許多人的電腦會有編碼的問題,無法顯示中文,故此處建議以Excel格式輸出。此處要求兩種格式都輸出的理由在處理下一個問題時會較明朗。

匯出後,你會在你的電腦的下載區(Downloads)看到如下的壓縮檔。請點擊開來檢視檔案。

0.5分

Q1g. 將名為Assignment 3 Wave 2的qsf檔import進一個新的Qualtrics問卷檔中。參考Q1e的步驟,生成1000則虛擬回應後以相同的方式匯出。事實上,若你不確定前面Q1b-Q1e各題怎麼做,你可以參考這份Wave 2問卷上的做法(點進去個版面參照)。兩份基本上是一樣的,Wave 2僅少了Wave 1的個資question block,因我們假設我們是對同樣一批樣本panel執行兩波調查。

0.5分

最後請在問卷上將我的email帳號列為collaborator,將問卷編輯權分享給我,方便我檢視。

Question 2. R語言分析

Q2a. Data Preparation

當你從壓縮檔中打開前面匯出的兩筆資料後,按照下方的區塊手動將資料切割。黃色那一列(row)和所有藍色的欄位(columns)全部刪除。沒錯,大膽的刪下去就對了! 只保留後面紅框標註的欄位,那是我們需要的資料。

下圖為Export labels的資料匯出格式,可以見到選項的文字內容皆被顯示。

下圖為Export values的資料匯出格式,可以見到選項的文字內容皆被轉化為數值,包括Drill-down題所收集的居住縣市(city)行政區(district)資訊。

接著,我們首先來看Wave 1的資料。我們將欄位依序命名為Gender, Age, City, District, PID, rating_2020, party_list_vote, treat.

接下來就是惱人的回應欄了。雖然Export values匯出格式替我們轉化成數值,可是縣市行政區怎麼辦? 看不懂230是什麼意思 @@ ?? 一個簡單直接的處理方式便是將Export labels的值複製貼Export values的縣市行政區欄位上。這樣就完成了。未來若當你收集上萬筆資料,此作法可能不是太實際,但本次作業希望各位能透過這種手作的方式熟習後面的原理。

另外,並請在最前方創建一個名為ID欄位,以下拉的方式依序生成1至1000,做為這些虛擬資料的ID。最後的Wave 1檔案會長成下面這個樣子。請以csv格式儲存。

也以同樣的方式產製Wave 2資料。欄位名依序為PID, rating_2024, party_list_vote, treat。並請加上ID欄位。製作好的Wave 2資料會長成下面這樣,同樣,請以csv格式儲存。Wave 2資料欄位會較少,這筆資料不需要所住縣市行政區的欄位資訊。

這樣便製作完成所需的兩個檔案了。為了下面分析題複製R code的方便,你可以選擇將此兩檔案命名為你熟記的檔案名稱,或是如範例般命名為w1_data和w2_data (這樣複製R code時也免換名稱 :p )

1分。

附帶一提的是,雖然本次未要求各位做勘誤和排除無效回應(畢竟是使用生成的虛擬資料),以下幾個匯出的資料欄位值得各位留意:

Finished
數值為1者,代表這位填答者完成了這份問卷並成功送出(submit)
Progress
數值為100者,代表這位填答者完成這份問卷(進展到最後)
Q_RecaptchaScore
偵測可疑機器人(bot)/程式填答,若數值低於0.5則將被標註(flag),建議排除此類回應。參見Qualtrics官網說明

所以到底是要測什麼?

此處我們將以這兩筆虛擬資料回答三個問題。首先,Wave 1 & 2問卷中的Treatment (將國內COVID疫情歸咎為執政黨無能)對填答者的政黨票投票選擇(我們的Y變數)是否具顯著影響?

其次,兩波問卷期間,填答者對執政者施政滿意度(參考rating這個變數)的變化是否影響受試者的政黨票投票選擇? Abramowitz (1988) 的Time-for-change模型,時至今日,還有多少解釋力?

最後,疫情期間各縣市染疫死亡人數不同,這個sub-national層次的變化是否影響填答者的政黨票投票選擇呢? 我們先以視覺化工具看看疫情期間各縣市染疫人數的分布(資料來源: CDC)

# remotes::install_github("shihjyun/twmap")    # Install twmap from GitHub
# install.packages(c("sf", "dplyr", "ggplot2")) # For spatial data and plotting

library(twmap)
library(sf)
library(dplyr)
library(ggplot2)

covid_counts <- read.csv(url("https://www.dropbox.com/scl/fi/ieva8y0nz3v1gy7ntuii1/covid_counts.csv?rlkey=rqhfftnw3weeo0ruwpc39xo5a&st=ebnskypj&dl=1"))

# Load county-level map
county_map <- tw_county   # Built-in from twmap package :contentReference[oaicite:2]{index=2}


# Merge your data with the map spatial data
map_data <- county_map %>%
  left_join(covid_counts, by = "COUNTYNAME")

# Plot a choropleth map of cases
ggplot(map_data) +
  geom_sf(aes(fill = cases), colour = "grey40", size = 0.1) +
  scale_fill_viridis_c(option = "plasma", na.value = "lightgrey", 
                       name = "累計確診數") +
  labs(title = "台灣各縣市 COVID-19 累計確診病例",
       subtitle = "資料期間:2020 年至某時點",
       caption = "資料來源:CDC") +
  theme_minimal() +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank())

# 確實存在相當大差異,但不清楚是否真會對填答者的政黨票投票選擇造成影響。

我們來實際分析一下吧。

Q2b. R programming & regression analysis 將資料輸入你的R studio或任何的介面,參考下方的步驟整併處理後。以party_list_vote作為Y變數,以Gender, Age, PID, rating,和treat作為X變數,估測一個logistic regression模型。另再以同樣的Y和X變數,納入各縣市疫情期間染疫人數的log()調整值,估測一個多層次logistic regression模型。比較兩者對填答者的政黨票投票選擇的預測能力。

## Step 1: Data wrangling 

# Load the data
w1_data <- read.csv("C://Users/hktse/Dropbox/Syllabus/Public opinion analysis/assignment3_w1.csv")
w2_data <- read.csv("C://Users/hktse/Dropbox/Syllabus/Public opinion analysis/assignment3_w2.csv")

# Check data info
dim(w1_data)
## [1] 1000    9
names(w1_data)
## [1] "ID"              "Gender"          "Age"             "city"           
## [5] "district"        "PID"             "rating_2020"     "party_list_vote"
## [9] "treat"
str(w1_data)
## 'data.frame':    1000 obs. of  9 variables:
##  $ ID             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Gender         : int  2 2 1 2 1 2 1 1 2 1 ...
##  $ Age            : int  3 3 2 2 2 2 2 3 3 2 ...
##  $ city           : chr  "台南市" "基隆市" "宜蘭縣" "連江縣" ...
##  $ district       : chr  "南區" "仁愛區" "員山鄉" "南竿鄉" ...
##  $ PID            : int  1 2 1 2 1 1 2 1 1 1 ...
##  $ rating_2020    : int  10 5 1 3 10 1 9 4 1 1 ...
##  $ party_list_vote: int  1 1 1 2 1 1 2 1 1 1 ...
##  $ treat          : int  0 1 1 0 1 0 0 1 1 0 ...
dim(w2_data)
## [1] 1000    5
names(w2_data)
## [1] "ID"              "PID"             "rating_2024"     "party_list_vote"
## [5] "treat"
str(w2_data)
## 'data.frame':    1000 obs. of  5 variables:
##  $ ID             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ PID            : int  1 1 2 1 2 2 2 2 2 2 ...
##  $ rating_2024    : int  8 8 5 8 6 6 9 5 1 6 ...
##  $ party_list_vote: int  2 2 2 1 2 1 1 1 2 1 ...
##  $ treat          : int  1 0 1 0 1 0 1 0 1 0 ...
## Step 2: Convert data "class"
# PID
w1_data$PID <- as.factor(w1_data$PID)
levels(w1_data$PID) <- c("KMT", "DPP")
w2_data$PID <- as.factor(w2_data$PID)
levels(w2_data$PID) <- c("KMT", "DPP")
 
# Gender (only for Wave 1 data)
w1_data$Gender <- as.integer(ifelse(w1_data$Gender == 2, 0, 1))

# Outcome variable
w1_data$party_list_vote <- as.factor(w1_data$party_list_vote)
levels(w1_data$party_list_vote) <- c("KMT", "DPP")

w2_data$party_list_vote <- as.factor(w2_data$party_list_vote)
levels(w2_data$party_list_vote) <- c("KMT", "DPP")
 
# Merge with covid_counts data (only for Wave 1 data), to be copied to Wave 2 data later
# Note that we are merging covid_counts to Wave 1 data, the resulting data.frame should have the same row length as Wave 1 data, so w1_data should be placed in the first argument inside left_join(), followed by covid_counts
# install.packages("dplyr")
library(dplyr)

# Remove COUNTYNAME column from covid_counts (we don't need it here)
covid_counts1 <- covid_counts[, -c(3)]

w1_data <- left_join(w1_data, covid_counts1, by = c("city"))  # Merge by city name


## Step 3: Merge personal attributes (Gender, Age, city, district) from Wave 1 data with W2's responses (since we assume these respondent attributes are fixed in the short run)
# Note that we are merging Wave 1's attributes to Wave 2, so Wave 2 should be placed on the first argument inside left_join() 
w2_data <- left_join(w2_data, w1_data[, c(1:5, 10)], by = c("ID"))

# Re-order Wave 2 data to make it resemble the order of Wave 1 data
w2_data <- w2_data[, c(1, 6:9, 2:5, 10)]

# Rename "rating_ " columns from both waves of data, which will be used as our outcome variable, they need to be the same name in order for the model estimated using Wave 1 data to be used as predictive model on Wave 2 data.

names(w1_data)[7] <- "rating"
names(w2_data)[7] <- "rating"

# The two data sets at a glance
head(w1_data)
##   ID Gender Age   city district PID rating party_list_vote treat  cases
## 1  1      0   3 台南市     南區 KMT     10             KMT     0 734368
## 2  2      0   3 基隆市   仁愛區 DPP      5             KMT     1 183765
## 3  3      1   2 宜蘭縣   員山鄉 KMT      1             KMT     1 198748
## 4  4      0   2 連江縣   南竿鄉 DPP      3             DPP     0   4077
## 5  5      1   2 金門縣   金城鎮 KMT     10             KMT     1  27035
## 6  6      0   2 彰化縣   田中鎮 KMT      1             KMT     0 480995
head(w2_data)
##   ID Gender Age   city district PID rating party_list_vote treat  cases
## 1  1      0   3 台南市     南區 KMT      8             DPP     1 734368
## 2  2      0   3 基隆市   仁愛區 KMT      8             DPP     0 183765
## 3  3      1   2 宜蘭縣   員山鄉 DPP      5             DPP     1 198748
## 4  4      0   2 連江縣   南竿鄉 KMT      8             KMT     0   4077
## 5  5      1   2 金門縣   金城鎮 DPP      6             DPP     1  27035
## 6  6      0   2 彰化縣   田中鎮 DPP      6             KMT     0 480995

現在請根據上方的兩筆資料,仿照下列的R code分別估測兩個模型(變數不顯著沒關係)。 1. 用party_list_vote作Y變數,以Gender, Age, PID, rating,和treat作為X變數,估測一個logistic regression模型。 2. 用同樣的Y和X變數,納入各縣市疫情期間染疫人數的log()調整值,估測一個多層次logistic regression模型。比較兩者對填答者的政黨票投票選擇的預測能力。

## Step 4: Estimate a fully-pooled logit model AND a random intercept model 

fp.fit <- glm(party_list_vote ~ Gender + Age + PID + rating + treat, 
                    family = binomial(link = "logit"),  # logit link function
                    data = w1_data)

# install.packages("lme4") # For multi-level / hierarchical / mixed effects models
library(lme4)

# To fit a multi-level logit model, you must specify higher level units' slope (regressors) and intercept 
# as (random_slope_predictor | grouping_variable), here "cases" are our random slopes, 
# which are grouped by "city" (the higher level unit)
# We apply log() transform on "cases" to normalize the range of covid counts
# family = binomial() means estimating a logit function

rs.fit <- glmer(
  party_list_vote ~ Gender + Age + PID + rating + treat + (log(cases) | city),
  data = w2_data,
  family = binomial()
)
## boundary (singular) fit: see help('isSingular')
# If "boundary (singular) fit: see help('isSingular')" message shows up, this means there're signs of overfitting, which is fine, just proceed as usual.


# Quickly inspect regression results
summary(fp.fit)
## 
## Call:
## glm(formula = party_list_vote ~ Gender + Age + PID + rating + 
##     treat, family = binomial(link = "logit"), data = w1_data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.551660   0.244176   2.259   0.0239 *
## Gender      -0.174371   0.127324  -1.370   0.1708  
## Age         -0.184617   0.077837  -2.372   0.0177 *
## PIDDPP       0.027198   0.127287   0.214   0.8308  
## rating      -0.008055   0.024423  -0.330   0.7415  
## treat       -0.161690   0.127252  -1.271   0.2039  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386.2  on 999  degrees of freedom
## Residual deviance: 1377.2  on 994  degrees of freedom
## AIC: 1389.2
## 
## Number of Fisher Scoring iterations: 3
summary(rs.fit)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: party_list_vote ~ Gender + Age + PID + rating + treat + (log(cases) |  
##     city)
##    Data: w2_data
## 
##      AIC      BIC   logLik deviance df.resid 
##   1394.3   1438.5   -688.2   1376.3      991 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2908 -1.0426  0.8156  0.9482  1.4734 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  city   (Intercept) 2.92076  1.7090        
##         log(cases)  0.01802  0.1343   -1.00
## Number of obs: 1000, groups:  city, 22
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.13514    0.24401   0.554    0.580
## Gender      -0.05337    0.12782  -0.418    0.676
## Age          0.01149    0.07806   0.147    0.883
## PIDDPP       0.01481    0.12796   0.116    0.908
## rating      -0.01527    0.02351  -0.650    0.516
## treat        0.11383    0.12800   0.889    0.374
## 
## Correlation of Fixed Effects:
##        (Intr) Gender Age    PIDDPP rating
## Gender -0.299                            
## Age    -0.655  0.015                     
## PIDDPP -0.244  0.028  0.000              
## rating -0.531  0.028  0.007 -0.028       
## treat  -0.258  0.007  0.019 -0.037  0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# As expected, all explanatory variables don't hold any water in both models ... Age is the only significant variable in the fully-pooled model. Treatment variable has no effect. Random effects aren't significant either.

0.5分

Q2c. 使用Wave 1建立的模型對Wave 2資料進行預測

## Step 5: Predict on Wave 2 data
fp.pred <- predict(fp.fit, w2_data, type = "response")
rs.pred <- predict(rs.fit, w2_data, type = "response", allow.new.levels = TRUE)

0.5分

Q2d. 根據課程講授的accuracy score計算方式,使用上述的預測結果對Wave 2資料的實際值(填答者在2024年大選中立委選舉政黨票投給藍綠哪一黨)作預測。請參照下列語法執行計算accuracy score,並詮釋所得之結果,對兩模型的預測準確度做比較。

# Re-lable predicted values with party names by mapping {"KMT", "DPP"} to [1, 0]
fp.pred <- ifelse(fp.pred > 0.5, 1, 0)
rs.pred <- ifelse(fp.pred > 0.5, 1, 0)

# Store predicted values into a dataframe object, along with ID and actual values

predictions <- data.frame(
  ID = w2_data$ID,
  actual = ifelse(w2_data$party_list_vote == "KMT", 1, 0),
  fp = fp.pred,
  rs = rs.pred
)

head(predictions)
##   ID actual fp rs
## 1  1      0  0  0
## 2  2      0  0  0
## 3  3      0  0  0
## 4  4      1  1  1
## 5  5      0  0  0
## 6  6      1  1  1
pred_bin <- predictions[, -c(1:2)]  # Remove ID column
actual_bin <- predictions$actual
accuracy_vals <- colMeans(pred_bin == actual_bin)
acc_df <- data.frame(Model = c("Fully-pooled", "Random slope"), Accuracy = accuracy_vals)

print(acc_df)
##           Model Accuracy
## fp Fully-pooled    0.525
## rs Random slope    0.525

1分

示範資料兩個模型的accuracy scores看來結果差不多,我們再來看看ROC-AUC曲線下面積。

Q2e. ROC-AUC計算及視覺化

請參照下列範例執行,並對產出的ROC-AUC曲線下的面積做出詮釋,比較兩模型對鑑別(a) 填答者在第二波問卷中表示自己在2024大選的政黨名單票中投給國民黨(Y = 1)及(b)投給民進黨(Y = 0)此兩種投票結果的鑑別能力。

## Calculate and plot ROC-AUC
#install.packages("pROC", "ggplot2")  # Uncomment if you haven't installed these packages
library(pROC)
library(ggplot2)

# Calculate the AUC for each model
auc_vals <- sapply(predictions[, -c(1:2)], function(x) {
  roc_obj <- roc(actual_bin, x, quiet = TRUE)
  as.numeric(auc(roc_obj))
})
auc_df <- data.frame(Model = names(auc_vals), AUC = round(auc_vals, 4))

print(auc_df)
##    Model   AUC
## fp    fp 0.522
## rs    rs 0.522
# ROC-AUC Plot
roc_list <- lapply(predictions[, -c(1:2)], function(x) roc(actual_bin, x, quiet = TRUE))
names(roc_list) <- c("Fully-pooled", "Random slope")
ggroc(roc_list) + 
  theme_minimal() + 
  labs(title = "Predicting party-list vote choices in 2024 General Election", color = "Model", x = "1 - Specificity", y = "Sensitivity")

# It doesn't seem that the more complex random slope model fares any better than plain vanilla fully-pooled model.
# Perhaps there's a lot of city/county-level heterogeneity that we should dig further into, let's see if predicted probability varies greatly across city and county?

1分

額外補充

如同上面的評論,或許對不分區政黨名單票的預測準確度在縣市層級存在極大的誤差變異量。我們來將縣市層級的平均accuracy score投放到縣市地圖上來看看。(此處不須答題,僅供你參考)

# Merge dataframe predictions with W1_data's "city" column by "ID"
pred_merged <- left_join(predictions, w1_data[, c(1, 4)], by = "ID")

pred <- pred_merged  # To save on typing
# Inspect
str(pred)
## 'data.frame':    1000 obs. of  5 variables:
##  $ ID    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ actual: num  0 0 0 1 0 1 1 1 0 1 ...
##  $ fp    : num  0 0 0 1 0 1 0 0 0 0 ...
##  $ rs    : num  0 0 0 1 0 1 0 0 0 0 ...
##  $ city  : chr  "台南市" "基隆市" "宜蘭縣" "連江縣" ...
head(pred)
##   ID actual fp rs   city
## 1  1      0  0  0 台南市
## 2  2      0  0  0 基隆市
## 3  3      0  0  0 宜蘭縣
## 4  4      1  1  1 連江縣
## 5  5      0  0  0 金門縣
## 6  6      1  1  1 彰化縣
# Normalize name function (台 & 臺 to match twmap)
normalize_name <- function(x) gsub("台", "臺", x, fixed = TRUE)

pred <- pred %>%
 mutate(COUNTYNAME = normalize_name(city))

# Compute city-level prediction accuracy for both fp and rs models
# accuracy = percentage of correct predictions

accuracy_city <- pred %>%
 group_by(COUNTYNAME) %>%
 summarise(
 fp_acc = mean(fp == actual) * 100,
 rs_acc = mean(rs == actual) * 100,
      n = n()
)

print(accuracy_city)
## # A tibble: 22 × 4
##    COUNTYNAME fp_acc rs_acc     n
##    <chr>       <dbl>  <dbl> <int>
##  1 南投縣       54.5   54.5    33
##  2 嘉義市       50     50      36
##  3 嘉義縣       52     52      50
##  4 基隆市       66.7   66.7    42
##  5 宜蘭縣       52.3   52.3    44
##  6 屏東縣       48.2   48.2    56
##  7 彰化縣       65     65      40
##  8 新北市       60.4   60.4    48
##  9 新竹市       46.2   46.2    39
## 10 新竹縣       63.8   63.8    47
## # ℹ 12 more rows
# Load county map
county_map <- tw_county %>%
  mutate(COUNTYNAME = normalize_name(COUNTYNAME))

# Merge accuracy data onto map
map_data <- county_map %>%
  left_join(accuracy_city, by = "COUNTYNAME")

# Plot fully-pooled (fp) model's accuracy map
p_fp <- ggplot(map_data) +
  geom_sf(aes(fill = fp_acc), color = "grey40") +
  scale_fill_viridis_c(option = "viridis", na.value = "lightgrey",
                       name = "FP模型\n正確率(%)") +
  labs(title = "台灣各縣市:Fully-pooled模型預測正確率",
       caption = paste0("資料筆數 = 1000 | 日期:", Sys.Date())) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.grid = element_blank())


# Plot Random Slope (rs) model's accuracy map
p_rs <- ggplot(map_data) +
  geom_sf(aes(fill = rs_acc), color = "grey40") +
  scale_fill_viridis_c(option = "plasma", na.value = "lightgrey",
                       name = "RS模型\n正確率(%)") +
  labs(title = "台灣各縣市:Random Slope模型預測正確率",
       caption = paste0("資料筆數 = 1000 | 日期:", Sys.Date())) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.grid = element_blank())


# Place the two map side-by-side using gridExtra
# install.packages("gridExtra")
library(gridExtra)

grid.arrange(p_fp, p_rs, ncol = 2)

Question 3. 你的小組期末問卷研究

重頭戲來了,請以兩個段落,約250字上下說明你的期末小組問卷研究想要做什麼? 內容必須包括下列要項:

  1. 研究主題: 你想研究什麼?
幾個主要的變數(解釋變數(x)、被解釋的現象(y)、相關控制變數、x-y之間的預期關係) 不需要完全合乎科學真理、也不需要洋洋灑灑的寫好幾頁的文獻回顧來證明x-y之間確實被學界透徹的研究過,我比較在乎你自己(和你的組員)怎麼看這件事!
  1. 如何把上述的假說或x-y關係在問卷上操作化?
什麼題型才是合適的? 要怎麼衡量? 類別尺度用multiple choices question,連續尺度用slider question? 尺度相同的一系列問題(如人格)用matrix question?
  1. 執行方法與步驟(過程)
你要怎麼執行? 誰是你研究的target group? 怎麼選取? 怎麼投放問卷? 你的受試者可能怎麼看待各題?
  1. 執行完畢後打算怎麼分析? x與y的關係適合用怎樣的模型分析?
什麼模型適合估測連續變數或二元變數? 或什麼模式適合探測不方便公開講的態度或多面向的選擇? 你不用真的去跑統計,但請解釋你的想法

以上問題沒有標準答案(你有寫助教都會給對),但請認真思考,因為這有助於你形塑期末報告的問題,使整個過程的理路更流暢。內容宜精簡不宜多,請避免以大量文獻回顧代替正文。

2分