rm(list=ls())
library(haven)
library(ggplot2)
#for my own desktop 
pretest <- read_sav("C:/Users/liopa/My Drive (liopaper@gmail.com)/垂死病中驚坐起，懺悔仍未寫論文/2_問卷/pretest02/pretest_combine_naomit.sav")
#for laptop mac
#pretest <- read_sav("/Volumes/GoogleDrive-115666120033367351851/我的雲端硬碟/垂死病中驚坐起，懺悔仍未寫論文/2_問卷/pretest02/pretest_combine_naomit.sav")
attach (pretest)
#合格受訪者
##0為未完成訪問//1為完整完成訪問
table(qulified)
## qulified
##  1 
## 52
unit <- qulified
unit[qulified == 1] <- 1
unit[qulified == 2] <- NA
table(unit)
## unit
##  1 
## 52
length(unit)
## [1] 52
#隨機分組##0為控制組//1為實驗組
treatment <- as.integer(treatment, levels(0:1))
#政治興趣
table(interest)
## interest
##  1  2  3  4 
##  3 32 15  2
int <- interest
int[interest == 1| interest == 2] <- 1
int[interest == 3| interest == 4] <- 0
int <- factor(int, levels = c(0,1))
table(int)
## int
##  0  1 
## 17 35
length(int)
## [1] 52

#政治討論

table(discuss)
## discuss
##  1  2  3  4 
##  6 28 17  1
length(discuss)
## [1] 52
length(na.omit(discuss))
## [1] 52

#民主滿意度

table(dem_satisfied)
## dem_satisfied
##  2  3  4 
## 22 26  4
ds <- dem_satisfied
ds[dem_satisfied == 1| dem_satisfied == 2] <- 1
ds[dem_satisfied == 3| dem_satisfied == 4] <- 0
ds <- factor(ds, levels = c(0,1))
table(ds)
## ds
##  0  1 
## 30 22
length(na.omit(ds))
## [1] 52

#社會公平

table(equlity)
## equlity
##  1  2  3  4 
##  2 17 30  3
equ <- equlity
equ[equlity == 1|equlity  == 2]<- 1
equ[equlity == 3|equlity  == 4] <- 0
equ <- factor(equ, levels = c(0,1))
table(equ)
## equ
##  0  1 
## 33 19
length(na.omit(equ))
## [1] 52

#政治惡鬥

table(vic)
## vic
##  1  2  3 
## 14 30  8
vic.bi <- vic
vic.bi[vic == 1|vic  == 2] <- 1
vic.bi[vic == 3|vic  == 4] <- 0
vic.bi <- factor(vic.bi, levels = c(0,1))
table(vic.bi)
## vic.bi
##  0  1 
##  8 44
length(na.omit(vic.bi))
## [1] 52

#體制偏好

table(pref_sys)
## pref_sys
##  1  2  3 
## 26 21  5
sys <- pref_sys
sys <- factor(sys, levels = c(1:3))
table(sys)
## sys
##  1  2  3 
## 26 21  5
length(na.omit(sys))
## [1] 52

#列項實驗題組

library(list)
## 載入需要的套件：sandwich
#1
table(list1)
## list1
##  1  2  3  4  5 
##  6 23 13  9  1
newlist1 <- list1 -1
newlist1 <- as.integer(newlist1)
table(newlist1)
## newlist1
##  0  1  2  3  4 
##  6 23 13  9  1
length(na.omit(newlist1))
## [1] 52
#2
table(list2)
## list2
##  1  2  3  4 
##  9 31  9  3
newlist2 <- list2 -1
newlist2 <- as.integer(newlist2)
table(newlist2)
## newlist2
##  0  1  2  3 
##  9 31  9  3
length(na.omit(newlist2))
## [1] 52
#3
table(list3)
## list3
##  1  2  3  4 
## 21 25  5  1
newlist3 <- list3 -1
newlist3 <- as.integer(newlist3)
table(newlist3)
## newlist3
##  0  1  2  3 
## 21 25  5  1
length(na.omit(newlist3))
## [1] 52
plot(as.factor(newlist1))

plot(as.factor(newlist2))

plot(as.factor(newlist3))

treatment[is.na(list1) == TRUE] <- NA
#NO DESIGN EFFECT 
#1
test.value.affirm.1 <- ict.test(na.omit(newlist1), na.omit(treatment), J = 3, gms = TRUE, pi.table = TRUE )
print(test.value.affirm.1)
## 
## Test for List Experiment Design Effects
## 
## Estimated population proportions 
##                            est.   s.e.
## pi(Y_i(0) = 0, Z_i = 1) -0.0952 0.0903
## pi(Y_i(0) = 1, Z_i = 1) -0.0476 0.1379
## pi(Y_i(0) = 2, Z_i = 1) -0.0476 0.1086
## pi(Y_i(0) = 3, Z_i = 1)  0.0417 0.0408
## pi(Y_i(0) = 0, Z_i = 0)  0.1667 0.0761
## pi(Y_i(0) = 1, Z_i = 0)  0.5119 0.1118
## pi(Y_i(0) = 2, Z_i = 0)  0.2976 0.1211
## pi(Y_i(0) = 3, Z_i = 0)  0.1726 0.0876
## 
##  Y_i(0) is the (latent) count of 'yes' responses to the control items. Z_i is the (latent) binary response to the sensitive item.
## 
## Bonferroni-corrected p-value
## If this value is below alpha, you reject the null hypothesis of no design effect. If it is above alpha, you fail to reject the null.
## 
## Sensitive Item 1 
##        0.6543061
#2
test.value.affirm.2 <- ict.test(na.omit(newlist2), na.omit(treatment), J = 3, gms = TRUE, pi.table = TRUE )
## Error in t.test.default(y[treat == 1] <= j, y[treat == 0] <= (j - 1),  : 
##   資料是恆量
print(test.value.affirm.2)
## 
## Test for List Experiment Design Effects
## 
## Estimated population proportions 
##                           est.   s.e.
## pi(Y_i(0) = 0, Z_i = 1) 0.0119 0.1050
## pi(Y_i(0) = 1, Z_i = 1) 0.0357 0.1176
## pi(Y_i(0) = 2, Z_i = 1) 0.1250 0.0675
## pi(Y_i(0) = 3, Z_i = 1) 0.0000 0.0000
## pi(Y_i(0) = 0, Z_i = 0) 0.1667 0.0761
## pi(Y_i(0) = 1, Z_i = 0) 0.5714 0.1142
## pi(Y_i(0) = 2, Z_i = 0) 0.0893 0.1028
## pi(Y_i(0) = 3, Z_i = 0) 0.0000 0.0000
## 
##  Y_i(0) is the (latent) count of 'yes' responses to the control items. Z_i is the (latent) binary response to the sensitive item.
## 
## Bonferroni-corrected p-value
## If this value is below alpha, you reject the null hypothesis of no design effect. If it is above alpha, you fail to reject the null.
## 
## Sensitive Item 1 
##                1
#3
test.value.affirm.3 <- ict.test(na.omit(newlist3), na.omit(treatment), J = 3, gms = TRUE, pi.table = TRUE )
print(test.value.affirm.3)
## 
## Test for List Experiment Design Effects
## 
## Estimated population proportions 
##                            est.   s.e.
## pi(Y_i(0) = 0, Z_i = 1)  0.4405 0.1196
## pi(Y_i(0) = 1, Z_i = 1)  0.0179 0.0893
## pi(Y_i(0) = 2, Z_i = 1) -0.0357 0.0351
## pi(Y_i(0) = 3, Z_i = 1)  0.0000 0.0000
## pi(Y_i(0) = 0, Z_i = 0)  0.1667 0.0761
## pi(Y_i(0) = 1, Z_i = 0)  0.2679 0.1143
## pi(Y_i(0) = 2, Z_i = 0)  0.1071 0.0585
## pi(Y_i(0) = 3, Z_i = 0)  0.0357 0.0351
## 
##  Y_i(0) is the (latent) count of 'yes' responses to the control items. Z_i is the (latent) binary response to the sensitive item.
## 
## Bonferroni-corrected p-value
## If this value is below alpha, you reject the null hypothesis of no design effect. If it is above alpha, you fail to reject the null.
## 
## Sensitive Item 1 
##        0.5356749
#列項第二題為no design effect

#背書實驗題組

#還沒寫

#統獨六分類

#政黨認同體組

table(party1)
## party1
##  1  2  5  6  7  8  9 10 11 
##  3 15  1  1  5  9  1  3 14
table(party2)
## party2
##  1  2 10 
##  1  3 13
newpartyID <- c()
newpartyID[party1 == 1|party2 == 1] <- 1
newpartyID[party1 == 2|party2 == 2] <- 2
newpartyID[party1 == 3|party2 == 3] <- 3
newpartyID[party1 == 4|party2 == 4] <- 4
newpartyID[party1 == 5|party2 == 5] <- 5
newpartyID[party1 == 6|party2 == 6] <- 6
newpartyID[party1 == 7|party2 == 7] <- 7
newpartyID[party1 == 8|party2 == 8] <- 8
newpartyID[party1 == 9|party2 == 9] <- 9
newpartyID[party1 == 10|party2 == 10] <- 10
newpartyID <- factor(newpartyID)
table(newpartyID)
## newpartyID
##  1  2  5  6  7  8  9 10 
##  4 17  1  1  5  9  1 14
length(newpartyID)
## [1] 52
#E1 PartyID recode to Party_Pref
party <- c()
party[newpartyID == 1 |newpartyID == 3 |newpartyID == 4] <- 1
party[newpartyID == 2|newpartyID == 5|newpartyID == 6|newpartyID == 9] <- 2 
party[newpartyID == 7] <- 3
party[newpartyID == 8] <- 4
party[newpartyID == 10] <-5
table(party)
## party
##  1  2  3  4  5 
##  4 20  5  9 14
length(party)
## [1] 52

#人口變數(性別、年紀、教育程度、國族認同、省籍、戶籍

table(birth)
## birth
## 41 42 49 61 65 66 69 70 71 72 74 76 77 78 79 80 81 82 83 84 85 87 88 89 90 
##  1  1  1  1  1  1  2  3  1  1  1  1  2  5  5  2  2  6  5  3  1  1  1  3  1
age <- 111-birth 

table(age)
## age
## 21 22 23 24 26 27 28 29 30 31 32 33 34 35 37 39 40 41 42 45 46 50 62 69 70 
##  1  3  1  1  1  3  5  6  2  2  5  5  2  1  1  1  1  3  2  1  1  1  1  1  1

#直接問題

table(direct1)
## direct1
##  1  2  3  4 
## 10 26 14  2
direct1.bi <- c()
direct1.bi[direct1 == 1|direct1  == 2] <- 1
direct1.bi[direct1 == 3|direct1  == 4] <- 0

table(direct1.bi)
## direct1.bi
##  0  1 
## 16 36
length(na.omit(direct1.bi))
## [1] 52
table(direct2)
## direct2
##  1  2  3  4 
## 15 26  8  3
direct2.bi <- c()
direct2.bi[direct2 == 1|direct2  == 2] <- 1
direct2.bi[direct2 == 3|direct2  == 4] <- 0

table(direct2.bi)
## direct2.bi
##  0  1 
## 11 41
length(na.omit(direct1.bi))
## [1] 52
table(direct3)
## direct3
##  1  2  3  4 
## 11 26 10  5
direct3.bi <- c()
direct3.bi[direct3 == 1|direct3  == 2] <- 1
direct3.bi[direct3 == 3|direct3  == 4] <- 0

table(direct3.bi)
## direct3.bi
##  0  1 
## 15 37
length(na.omit(direct1.bi))
## [1] 52

#威權迫害經驗

table(persecution)
## persecution
##  1  2 
##  2 50
pretest.frame <- data.frame(treatment, list1, list2, list3, direct1.bi, direct2.bi, direct3.bi,edu , sex , age , int )
test_martix <- cbind(treatment, list1, list2, list3, direct1.bi, direct2.bi, direct3.bi,edu , sex , age , int)
test_data <- as.data.frame(test_martix)
list1n <- na.omit(list1)
treatn <- na.omit(treatment)
length(list1n)
## [1] 52
length(treatn)
## [1] 52
#已將一開始有遺漏值之問卷先移除，上開na.omit在未移除之情況仍有必要
#Dif in means in list.1
diff.in.means.results.1 <- ictreg( list1 ~ 1, data = pretest.frame, treat = "treatment", J=3, method = "lm")
summary(diff.in.means.results.1)
## 
## Item Count Technique Regression 
## 
## Call: ictreg(formula = list1 ~ 1, data = pretest.frame, treat = "treatment", 
##     J = 3, method = "lm")
## 
## Sensitive item 
##                 Est.    S.E.
## (Intercept) -0.14881 0.28316
## 
## Control items 
##                Est.    S.E.
## (Intercept) 2.60714 0.17638
## 
## Residual standard error: 0.986275 with 50 degrees of freedom
## 
## Number of control items J set to 3. Treatment groups were indicated by '' and '' and the control group by ''.
# calculate __ from direct question
summary(direct1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   2.000   2.154   3.000   4.000
library(ggplot2)
# combine list & direct question 1
q1.ld <- combinedListDirect(list1 ~ edu + sex + age + int , data =  test_data, treat = "treatment", direct = "direct1.bi")
summary(q1.ld)
## 
##  Combined List Estimates 
## 
## Call: combinedListDirect(formula = list1 ~ edu + sex + age + int, data = test_data, 
##     treat = "treatment", direct = "direct1.bi")
## 
##  Prevalence estimates
##                 Combined     Direct Conventional
## Estimate       0.9895642 0.69230769   -0.0831507
## Standard Error 0.1203899 0.06462831    0.2597272
## 
##  Placebo Test I
##        Beta is the conventional list experiment estimate among those who answer 'Yes' to the direct question.
##        Ho: beta = 1
##        Ha: beta != 1
##       
##         Estimate        SE            p  n
## beta -0.3409641 0.2254772 2.727254e-09 36
## 
##  Placebo Test II
##        Delta is the average effect of the receiving the treatment list on the direct question response.
##        Ho: delta = 0
##        Ha: delta != 0
##       
##          Estimate        SE         p  n
## delta 0.02644419 0.1291335 0.8377426 52
names(q1.ld)
##  [1] "comb.est"   "comb.se"    "direct.est" "direct.se"  "conv.est"  
##  [6] "conv.se"    "placebo.I"  "placebo.II" "call"       "data"      
## [11] "x"          "y"          "treat"      "direct"
q1_LD.est <- q1.ld$comb.est
q1_LD.se <- q1.ld$comb.se

q1_D.est <- q1.ld$direct.est
q1_D.se <- q1.ld$direct.se

(q1.bias.est <- q1_LD.est - q1_D.est)
##    t.nona 
## 0.2972565
(q1.bias.se <- sqrt((q1_LD.se)^2+(q1_D.se)^2))
##    t.nona 
## 0.1366402
(q1.t.ratio <- q1.bias.est/ q1.bias.se)
##   t.nona 
## 2.175468
x <- c("1. L&D", "2. Direct")
y <- c(q1_LD.est, q1_D.est)
se <- c(q1_LD.se, q1_D.se)
dat <- data.frame(x=x, y=y, se=se)
base <- ggplot(dat, aes(x, y, ymin=y-qnorm(.025)*se, ymax=y+qnorm(.025)*se)) + geom_pointrange()
base + theme(axis.text.x=element_text(angle=-30, vjust=1, hjust=0)) + xlab(NULL) + ylim(0, 1) + ylab("Estimated Proportion")
## Warning: Removed 1 rows containing missing values (geom_segment).