R語法檔，轉為Qmd

Author

王庭宥 M146020011

研究動機

2016年4月，蘇起先生表示，依數據顯示，中共不會打、美國不會救的狀況下，台獨支持比例接近六成；若中共不打、美國會救，支持比例提高。
本練習嘗試使用 TNSS2015 公開資料重製此發現，並進一步檢查資料背後的意義。

研究問題

  • 中共攻台與美國援台對支持台獨傾向是否有影響？
  • 對民眾來說，這兩個因素之間是否相互牽制？

研究假設

  • 假設1：若民眾認知到中共會攻打台灣，會降低支持台獨的可能。
  • 假設2：若民眾認知到美國會救台灣，會增加支持台獨的可能。
  • 假設3：中共攻台的認知具有調節美國援台對支持台獨影響的作用。

資料描述

使用 The Taiwan National Security Surveys
這個資料為 Duke University 與政大選舉研究中心合作的研究案資料，受訪對象為台灣地區（不含金門、馬祖）年滿二十歲以上的成年人。
抽樣是以「中華電信住宅部103-104年版電話號碼簿」為母體清冊，依據各縣市電話簿所刊電話數佔台灣地區所刊電話總數比例，決定各縣市抽出之電話個數。
先以等距抽樣法抽出各縣市電話樣本後，為求完整的涵蓋性，再隨機修正最後二碼及四碼，以求接觸到未登錄電話的住宅戶。
電話接通後再由訪員按照戶中抽樣的原則，抽出應受訪的對象進行訪問。
調查期間為2015年10月27日至10月30日，以及11月02日和03日，由政治大學選舉研究中心執行。
最後完成的有效問卷1,071份。
若以百分之九十五之信賴度估計，最大可能抽樣誤差為：±2.99%。

載入套件與資料

library(here)
here() starts at C:/Users/wang9/OneDrive/桌面/庭宥的R project
library(sjmisc)
library(sjPlot)
library(gmodels)
library(car)
Warning: package 'car' was built under R version 4.5.2
Loading required package: carData
Warning: package 'carData' was built under R version 4.5.2
library(VGAM)
Warning: package 'VGAM' was built under R version 4.5.2
Loading required package: stats4
Loading required package: splines
here::here()
[1] "C:/Users/wang9/OneDrive/桌面/庭宥的R project"
load(here("上課講義", "TNSS2015.rda"))

編碼

依變數：獨立傾向

12.關於臺灣和大陸的關係，有下面幾種不同的看法：
1: 儘快統一
2: 儘快宣布獨立
3: 維持現狀，以後走向統一
4: 維持現狀，以後走向獨立
5: 維持現狀，看情形再決定獨立或統一
6: 永遠維持現狀
請問您比較偏向那一種？

TNSS2015$Q12n <- rec(TNSS2015$Q12, rec = "2,4=1;else=0", as.num = TRUE)
frq(TNSS2015$Q12n, weights = TNSS2015$w)
關於臺灣和大陸的關係，有下面幾種不同的看法，請問您比較偏向那一種？ (xw) <numeric> 
# total N=1071 valid N=1071 mean=0.21 sd=0.41

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
    0 | 841 | 78.52 |   78.52 |  78.52
    1 | 230 | 21.48 |   21.48 | 100.00
 <NA> |   0 |  0.00 |    <NA> |   <NA>
frq(rec(TNSS2015$Q12, rec = "2,4,5=1;else=0", as.num = TRUE),
    weights = TNSS2015$w)
關於臺灣和大陸的關係，有下面幾種不同的看法，請問您比較偏向那一種？ (xw) <numeric> 
# total N=1071 valid N=1071 mean=0.58 sd=0.49

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
    0 | 450 | 42.02 |   42.02 |  42.02
    1 | 621 | 57.98 |   57.98 | 100.00
 <NA> |   0 |  0.00 |    <NA> |   <NA>

條件獨立1：武力犯台

13.如果臺灣宣佈獨立會引起大陸攻打臺灣，請問您贊不贊成臺灣獨立？【訪員請追問強弱程度】
01.非常不贊成
02.不贊成
03.贊成
04.非常贊成
96.看情形
97.無意見
98.不知道
95.拒答

TNSS2015$Q13n <- rec(TNSS2015$Q13, rec = "1,2=0;3,4=1; else=NA", as.num = F)
frq(TNSS2015$Q13n, weights = TNSS2015$w) # 37.3%
如果臺灣宣佈獨立會引起大陸攻打臺灣，請問您贊不贊成臺灣獨立？ (xw) <categorical> 
# total N=912 valid N=912 mean=0.37 sd=0.48

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
    0 | 572 | 62.72 |   62.72 |  62.72
    1 | 340 | 37.28 |   37.28 | 100.00
 <NA> |   0 |  0.00 |    <NA> |   <NA>

條件獨立2：不會犯台

14.那如果臺灣宣佈獨立，而大陸不會攻打臺灣，請問您贊不贊成臺灣獨立？【訪員請追問強弱程度】
01.非常不贊成
02.不贊成
03.贊成
04.非常贊成
96.看情形
97.無意見
98.不知道
95.拒答

TNSS2015$Q14n <- rec(TNSS2015$Q14, rec = "1,2=0;3,4=1; else=NA", as.num = F)
frq(TNSS2015$Q14n, weights = TNSS2015$w) # 78.95%
那如果臺灣宣佈獨立，而大陸不會攻打臺灣，請問您贊不贊成臺灣獨立？ (xw) <categorical> 
# total N=931 valid N=931 mean=0.79 sd=0.41

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
    0 | 196 | 21.05 |   21.05 |  21.05
    1 | 735 | 78.95 |   78.95 | 100.00
 <NA> |   0 |  0.00 |    <NA> |   <NA>

資料描述與解讀

見本書正文。

自變數編碼

自變數1：對大陸攻台的認知

27.如果臺灣自行宣佈獨立，請問您認為大陸會不會攻打臺灣？【訪員請追問強弱程度】
01.一定不會
02.不會
03.會
04.一定會
96.看情形
97.無意見
98.不知道
95.拒答

TNSS2015$Q27n <- rec(TNSS2015$Q27, rec = "1,2=0;3,4=1;else=NA", as.num = F)
frq(TNSS2015$Q27n, weights = TNSS2015$w)   # 60.9% 認知大陸會攻打
如果臺灣自行宣佈獨立，請問您認為大陸會不會攻打臺灣？ (xw) <categorical> 
# total N=929 valid N=929 mean=0.61 sd=0.49

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
    0 | 363 | 39.07 |   39.07 |  39.07
    1 | 566 | 60.93 |   60.93 | 100.00
 <NA> |   0 |  0.00 |    <NA> |   <NA>

自變數2：對美國援助的認知

30.如果因為臺灣宣佈獨立，大陸攻打臺灣，請問您認為美國會不會出兵幫助臺灣？【訪員請追問強弱程度】
01.一定不會
02.不會
03.會
04.一定會
96.看情形
97.無意見
98.不知道
95.拒答

TNSS2015$Q30n <- rec(TNSS2015$Q30, rec = "1,2=0;3,4=1; else=NA", as.num = F)
frq(TNSS2015$Q30n, weights = TNSS2015$w)   # 70.2% 認知美國會救
如果因為臺灣宣佈獨立，大陸攻打臺灣，請問您認為美國會不會出兵幫助臺灣？ (xw) <categorical> 
# total N=900 valid N=900 mean=0.70 sd=0.46

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
    0 | 268 | 29.78 |   29.78 |  29.78
    1 | 632 | 70.22 |   70.22 | 100.00
 <NA> |   0 |  0.00 |    <NA> |   <NA>

自變數與依變數之間的相關性檢證

獨立立場 vs 美國因素

library(gmodels)
CrossTable(TNSS2015$Q12n, TNSS2015$Q30n,
           prop.r = TRUE, prop.t = FALSE, prop.c = TRUE,
           prop.chisq = FALSE, chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Row Total |
|           N / Col Total |
|-------------------------|

 
Total Observations in Table:  916 

 
              | TNSS2015$Q30n 
TNSS2015$Q12n |         0 |         1 | Row Total | 
--------------|-----------|-----------|-----------|
            0 |       251 |       451 |       702 | 
              |     0.358 |     0.642 |     0.766 | 
              |     0.860 |     0.723 |           | 
--------------|-----------|-----------|-----------|
            1 |        41 |       173 |       214 | 
              |     0.192 |     0.808 |     0.234 | 
              |     0.140 |     0.277 |           | 
--------------|-----------|-----------|-----------|
 Column Total |       292 |       624 |       916 | 
              |     0.319 |     0.681 |           | 
--------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  20.80134     d.f. =  1     p =  5.094721e-06 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  20.04412     d.f. =  1     p =  7.567574e-06

獨立立場 vs 中共因素

CrossTable(TNSS2015$Q12n, TNSS2015$Q27n,
           prop.r = TRUE, prop.t = FALSE, prop.c = TRUE,
           prop.chisq = FALSE, chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Row Total |
|           N / Col Total |
|-------------------------|

 
Total Observations in Table:  942 

 
              | TNSS2015$Q27n 
TNSS2015$Q12n |         0 |         1 | Row Total | 
--------------|-----------|-----------|-----------|
            0 |       237 |       489 |       726 | 
              |     0.326 |     0.674 |     0.771 | 
              |     0.707 |     0.806 |           | 
--------------|-----------|-----------|-----------|
            1 |        98 |       118 |       216 | 
              |     0.454 |     0.546 |     0.229 | 
              |     0.293 |     0.194 |           | 
--------------|-----------|-----------|-----------|
 Column Total |       335 |       607 |       942 | 
              |     0.356 |     0.644 |           | 
--------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  11.76451     d.f. =  1     p =  0.0006037061 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  11.21574     d.f. =  1     p =  0.0008110655

美國因素 vs 中共因素

CrossTable(TNSS2015$Q30n, TNSS2015$Q27n,
           prop.r = TRUE, prop.t = FALSE, prop.c = TRUE,
           prop.chisq = FALSE, chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Row Total |
|           N / Col Total |
|-------------------------|

 
Total Observations in Table:  842 

 
              | TNSS2015$Q27n 
TNSS2015$Q30n |         0 |         1 | Row Total | 
--------------|-----------|-----------|-----------|
            0 |        72 |       200 |       272 | 
              |     0.265 |     0.735 |     0.323 | 
              |     0.240 |     0.369 |           | 
--------------|-----------|-----------|-----------|
            1 |       228 |       342 |       570 | 
              |     0.400 |     0.600 |     0.677 | 
              |     0.760 |     0.631 |           | 
--------------|-----------|-----------|-----------|
 Column Total |       300 |       542 |       842 | 
              |     0.356 |     0.644 |           | 
--------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  14.6958     d.f. =  1     p =  0.0001263278 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  14.11181     d.f. =  1     p =  0.0001722588

假設檢證

mod.1 <- glm(Q12n ~ Q27n + Q30n, data = TNSS2015, family = "binomial")
summary(mod.1)

Call:
glm(formula = Q12n ~ Q27n + Q30n, family = "binomial", data = TNSS2015)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.4167     0.2039  -6.947 3.74e-12 ***
Q27n1        -0.4430     0.1674  -2.646 0.008155 ** 
Q30n1         0.7226     0.1952   3.702 0.000214 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 923.20  on 841  degrees of freedom
Residual deviance: 898.54  on 839  degrees of freedom
  (229 observations deleted due to missingness)
AIC: 904.54

Number of Fisher Scoring iterations: 4
mod.2 <- update(mod.1, . ~ . + Q27n:Q30n)
summary(mod.2)

Call:
glm(formula = Q12n ~ Q27n + Q30n + Q27n:Q30n, family = "binomial", 
    data = TNSS2015)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.5126     0.3064  -4.937 7.94e-07 ***
Q27n1        -0.3027     0.3680  -0.823   0.4107    
Q30n1         0.8391     0.3369   2.491   0.0127 *  
Q27n1:Q30n1  -0.1778     0.4136  -0.430   0.6673    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 923.20  on 841  degrees of freedom
Residual deviance: 898.36  on 838  degrees of freedom
  (229 observations deleted due to missingness)
AIC: 906.36

Number of Fisher Scoring iterations: 4
library(car)
vif(mod.1)
    Q27n     Q30n 
1.010622 1.010622 
vif(mod.2)
there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif
     Q27n      Q30n Q27n:Q30n 
 4.871832  3.016957  6.173798 
exp(coef(mod.1))      # 印出每個變數的勝算
(Intercept)       Q27n1       Q30n1 
  0.2425042   0.6421274   2.0596841 
exp(confint(mod.1))   # 印出每個變數的95%信賴區間
Waiting for profiling to be done...
                2.5 %    97.5 %
(Intercept) 0.1607418 0.3580860
Q27n1       0.4626975 0.8925051
Q30n1       1.4165113 3.0493027

進一步探索

假設1a：相對於維持現狀的民眾來說，若受訪者認知到中共會攻打台灣，會降低支持台獨的可能。

假設2a：相對於維持現狀的民眾來說，若受訪者認知到美國會救台灣，會增加支持台獨的可能。

假設3a：相對於維持現狀的民眾來說，中共攻台的認知具有調節美國援台對支持台獨影響的作用（反之亦然）。

12.關於臺灣和大陸的關係，有下面幾種不同的看法：
1: 儘快統一
2: 儘快宣布獨立
3: 維持現狀，以後走向統一
4: 維持現狀，以後走向獨立
5: 維持現狀，看情形再決定獨立或統一
6: 永遠維持現狀
請問您比較偏向那一種？

library(sjPlot)
TNSS2015$Q12m <- rec(TNSS2015$Q12, rec = "1,3=2; 2,4,=1; 5,6=0; else=NA", as.num = F)
contrasts(TNSS2015$Q12m)
  1 2
0 0 0
1 1 0
2 0 1
library(VGAM)
mod.3 <- vglm(Q12m ~ Q27n + Q30n, data = TNSS2015, family = multinomial)
summary(mod.3)

Call:
vglm(formula = Q12m ~ Q27n + Q30n, family = multinomial, data = TNSS2015)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1   1.7409     0.2761   6.306 2.87e-10 ***
(Intercept):2   0.5207     0.3153   1.651  0.09868 .  
Q27n1:1        -0.2201     0.2737  -0.804  0.42128    
Q27n1:2        -0.6531     0.2970  -2.199  0.02787 *  
Q30n1:1         0.6339     0.2454   2.583  0.00978 ** 
Q30n1:2         1.2783     0.2887   4.428 9.52e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 1354.544 on 1652 degrees of freedom

Log-likelihood: -677.2722 on 1652 degrees of freedom

Number of Fisher scoring iterations: 5 


Reference group is level  3  of the response
mod.4 <- update(mod.3, . ~ . + Q27n:Q30n)
summary(mod.4)

Call:
vglm(formula = Q12m ~ Q27n + Q30n + Q27n:Q30n, family = multinomial, 
    data = TNSS2015)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1   1.6946     0.3627   4.673 2.97e-06 ***
(Intercept):2   0.3677     0.4336   0.848  0.39643    
Q27n1:1        -0.1613     0.4148  -0.389  0.69733    
Q27n1:2        -0.4367     0.5070  -0.861  0.38906    
Q30n1:1         0.7183     0.4715   1.524  0.12762    
Q30n1:2         1.4912     0.5332   2.796  0.00517 ** 
Q27n1:Q30n1:1  -0.1137     0.5516  -0.206  0.83664    
Q27n1:Q30n1:2  -0.3113     0.6345  -0.491  0.62371    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 1354.238 on 1650 degrees of freedom

Log-likelihood: -677.1188 on 1650 degrees of freedom

Number of Fisher scoring iterations: 5 


Reference group is level  3  of the response

補充盒子：用contrasts()直接從虛擬變數中調整對照組

contrasts(TNSS2015$Q12m) <- contr.treatment(levels(TNSS2015$Q12m), base = 2)