1 + 1[1] 2Untitled
Quarto
Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.
Running Code
When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:
You can add options to executable code like this
[1] 4The echo: false option disables the printing of code (only output is displayed).
df <- read.csv("UCT_subset.csv", na.strings = c("NA", ""))df$treatment <- ifelse(
df$payment_type == "control", 0L,
ifelse(df$payment_type %in% c("monthly", "lump sum"), 1L, NA_integer_)
)
df$treatment_f <- factor(
df$treatment,
levels = c(0, 1),
labels = c("Control", "Treatment")
)
# 检查
table(df$payment_type, df$treatment_f, useNA = "ifany")
Control Treatment
control 157 0
lump sum 0 46
monthly 0 47cons_dat <- df[complete.cases(
df[, c("treatment", "cons_nondurable_ppp0", "cons_nondurable_ppp1")]
), ]
cons_dat$cons_change <- cons_dat$cons_nondurable_ppp1 - cons_dat$cons_nondurable_ppp0
# 样本量
nrow(cons_dat)[1] 250table(cons_dat$treatment_f)
Control Treatment
157 93 tapply(
cons_dat$cons_change,
cons_dat$treatment_f,
function(x) c(n = length(x), mean = mean(x), sd = sd(x))
)$Control
n mean sd
157.00000 48.96124 145.69846
$Treatment
n mean sd
93.00000 32.62212 154.23221 tapply(cons_dat$cons_nondurable_ppp0, cons_dat$treatment_f, mean) Control Treatment
98.23722 171.52908 tapply(cons_dat$cons_nondurable_ppp1, cons_dat$treatment_f, mean) Control Treatment
147.1985 204.1512 boxplot(
cons_change ~ treatment_f,
data = cons_dat,
xlab = "Group",
ylab = "Change in non-durable consumption (PPP)",
main = "Consumption change by treatment group"
)
stripchart(
cons_change ~ treatment_f,
data = cons_dat,
vertical = TRUE,
method = "jitter",
add = TRUE,
pch = 16,
col = rgb(0, 0, 0, 0.3)
)
m_ancova <- lm(
cons_nondurable_ppp1 ~ treatment + cons_nondurable_ppp0,
data = cons_dat
)
summary(m_ancova)
Call:
lm(formula = cons_nondurable_ppp1 ~ treatment + cons_nondurable_ppp0,
data = cons_dat)
Residuals:
Min 1Q Median 3Q Max
-199.37 -72.82 -6.20 46.52 349.69
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 135.54883 9.70161 13.972 < 2e-16 ***
treatment 48.26130 13.96910 3.455 0.000648 ***
cons_nondurable_ppp0 0.11859 0.05304 2.236 0.026269 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 102.5 on 247 degrees of freedom
Multiple R-squared: 0.08524, Adjusted R-squared: 0.07783
F-statistic: 11.51 on 2 and 247 DF, p-value: 1.665e-05confint(m_ancova) 2.5 % 97.5 %
(Intercept) 116.44039177 154.6572628
treatment 20.74755255 75.7750428
cons_nondurable_ppp0 0.01411096 0.2230624# 只看 treatment 这一行
coef(summary(m_ancova))["treatment", ] Estimate Std. Error t value Pr(>|t|)
48.261297671 13.969099821 3.454860964 0.000647867 confint(m_ancova)["treatment", ] 2.5 % 97.5 %
20.74755 75.77504 m_ancova$df.residual[1] 247t_cons <- t.test(
cons_dat$cons_change[cons_dat$treatment == 1],
cons_dat$cons_change[cons_dat$treatment == 0],
var.equal = FALSE
)
t_cons
Welch Two Sample t-test
data: cons_dat$cons_change[cons_dat$treatment == 1] and cons_dat$cons_change[cons_dat$treatment == 0]
t = -0.82631, df = 184.56, p-value = 0.4097
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-55.35026 22.67201
sample estimates:
mean of x mean of y
32.62212 48.96124 par(mfrow = c(1, 1)) # 重置为 1 张图
plot(m_ancova, which = 1) # 残差 vs 拟合图
plot(m_ancova, which = 2) # QQ 图
# 保存高清箱线图
png("boxplot_high_res.png", width = 800, height = 600, res = 150)
boxplot(
cons_change ~ treatment_f,
data = cons_dat,
xlab = "Group",
ylab = "Change in non-durable consumption (PPP)",
main = "Consumption change by treatment group"
)
stripchart(
cons_change ~ treatment_f,
data = cons_dat,
vertical = TRUE,
method = "jitter",
add = TRUE,
pch = 16,
col = rgb(0, 0, 0, 0.3)
)
dev.off()png
2