Conventionally, group-level treatments are assigned in creating_session:
for g in subsession.get_groups():
g.treatment = random.choice([True, False])
However, this doesn't work when using group_by_arrival_time,
because groups are not determined until players arrive at the wait page.
(All players are in the same group initially.)
Instead, you need to assign treatments in after_all_players_arrive.
{{ extends "global/Page.html" }}
{{ block title }}Contribute{{ endblock }}
{{ block content }}
<p>
You can choose "invest" or "not invest".
</p><p>
{{C.NUM_SUCCESS}} players among those who invest will get a profit.
</p><p>
You will get
<ul><strong>{{C.PROFIT}}, when you win.</strong></ul><ul><strong>- 100 , when you lose.</strong></ul><ul><strong>0 , when you do not invest.</strong></ul>
{{ formfields }}
{{ next_button }}
{{ endblock }}
Results.html
最後は結果表示と結果に対する満足度を聞くフォームを表示するhtmlです.
{{ block title }}Results{{ endblock }}
{{ block content }}
<p>
Your profit is {{ player.profit }}.
</p><p>
Please answer the following questions.
</p>
{{ formfields }}
{{ next_button }}
{{ endblock }}
classPlayer(BasePlayer):
action = models.IntegerField(
choices=[[0, 'no'],[1,'yes (invest 100)']],
label="Do you invest? ",initial=(0)
)
outcome = models.IntegerField(
choices=[[5, 'very satisfied'],[4, 'satisfied'],[3, 'neither'],[2, 'dissatisfied'], [1, 'very dissatisfied']],
label='How do you feel about the result?',
widget=widgets.RadioSelect,
)
profit = models.IntegerField()
プレイヤーのクラスです.ここでは
action #行動
outcome #結果への評価(満足度)
profit #受け取る利得
の3属性を定義します.idのような基本属性はもともと定義されているので,呼び出すだけで使えます.
defset_profit(group: Group):#グループオブジェクトを引数にした関数
players = group.get_players()#playerを取得
group.x = sum([p.action for p in players])#投資者数の総数xを計算
xids=[]#投資者のplayeridを格納するリストfor p in players:
if p.action==1:
xids.append(p.id_in_group)#投資者のplayeridを全てリストに追加#投資者数が多い場合の抽選if C.NUM_SUCCESS < group.x:
win=random.sample(xids,C.NUM_SUCCESS)#NUM_SUCCESSだけxidsから勝者を選択for p in players:
if p.action == 0:
p.profit=0# 非投資者の利得を計算else:
if C.NUM_SUCCESS > group.x:
p.profit= 100*(C.RETURN_RATE*p.action-p.action)
else:
if p.id_in_group in win:
p.profit= 100*(C.RETURN_RATE*p.action-p.action)
else:
p.profit= -100
> summary(lm(data=data,y~d))
Call:
lm(formula = y ~ d, data = data)
Residuals:
Min 1Q Median 3Q Max
-12.3599-2.72840.03752.879111.4583
Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept)0.82420.16994.8521.42e-06***
d 11.81070.238449.549<2e-16***---
Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error:3.768 on 998 degrees of freedom
Multiple R-squared:0.711, Adjusted R-squared:0.7107F-statistic:2455 on 1 and 998 DF, p-value:<2.2e-16
介入効果の推定値は11.8107でした.真の効果は3ですから,過大推定になっています.
回帰ではなく,条件付き期待値の差をとっても同じ結果を得ます.
mean(y[d==1])-mean(y[d==0])
傾向スコアマッチング
傾向スコアはすでに計算済みですから,それをつかってマッチングデータを作り,ATEを推定します.
> matched.data <- Match(Y = data$y, Tr=data$d, X =p,estimand="ATE")> summary(matched.data)#Matched number of observations............... 500
Estimate... 3.4254
AI SE...... 0.23556T-stat..... 14.541
p.val...... <2.22e-16
Original number of observations.............. 1000
Original number of treated obs............... 508
Matched number of observations............... 1000
Matched number of observations (unweighted). 1152
matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE")
summary(matched.data)
Estimate... 7.2681
AI SE...... 0.10906T-stat..... 66.645
p.val...... <2.22e-16
Original number of observations.............. 3000
Original number of treated obs............... 1509
Matched number of observations............... 3000
Matched number of observations (unweighted). 7892
matched.data2 <- Match(Y = y, Tr=d, X =p.predict2,estimand="ATE")
summary(matched.data2)
Estimate... 3.4448
AI SE...... 0.1594T-stat..... 21.611
p.val...... <2.22e-16
Original number of observations.............. 3000
Original number of treated obs............... 1509
Matched number of observations............... 3000
Matched number of observations (unweighted). 4921
n <-3000;# sample size
tau <-3;# 真の介入効果
x1 <- runif(n,-5,5)# 共変量
x2 <- runif(n,-5,5)# 共変量
y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける
y1 <- y0+tau #潜在的結果y1
p <-1/(1+exp(-x1-x2))# make probabilities by logistic function
d <- rbinom(n,size=1,p)# treatment with bias
y <- c()for(i in1:n){# 観測値yの生成if(d[i]==1){y[i]<- y1[i]}else{y[i]<- y0[i]}}
dat <-data.frame(d,x1,x2,y)
> IPW_result
# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr><dbl><dbl><dbl><dbl>1(Intercept)-2.160.0679-31.94.56e-1922 d 7.450.095977.70.> IPW_result2
# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr><dbl><dbl><dbl><dbl>1(Intercept)-0.2690.0923-2.923.58e-32 d 4.020.13430.16.63e-174