[問題] 使用R計算effective size
如果是跟統計軟體有關請重發文章，使用程式做為分類。
統計軟體，如SPSS, AMOS, SAS, R, STATA, Eviews，請都使用程式做為分類
請詳述問題內容，以利板友幫忙解答，過短文章依板規處置，請注意。
為避免版面混亂，請勿手動置底問題，擅用E做檔案編輯
我人在國外，這學期的統計必修課是我完全沒學過的R，英文加上統計加上新軟體，真的
很吃力啊.... 以下問題是英文的，我已經花了一整天試著理解但仍有很多疑問，主要是
我不知道我該餵給R什麼訊息讓它去跑
因為這週得了flu所以缺課，頭腦昏沉沉的理解困難，實在無法才上來請高手幫忙指點
老師給的問題如下：
Identify an experimental study with a treatment group, a control group, and
some binary outcome (some sort of success/failure outcome) and calculate the
power for that experiment with different reasonable effect sizes (say, 5%,
10%, or 15%). Check your work using a power calculator (take a screenshot of
the calculator or provide me with the information I need to check what you
found). If you are having difficulty identifying an experiment, you can use
the Minneapolis Domestic Violence Experiment; the results for that experiment
are summarized here:
https://www.policefoundation.org/publication/the-minneapolis-domestic-violence-experiment/
If you use the Minneapolis Domestic Violence Experiment, be sure to compare
the arrest treatment to both of the non-arrest treatments (separate and
mediate) combined. Also, you can use the code I sent but if you run into
difficulty with the rounding error problem, send me your code and I will help
you troubleshoot it.
他給的R code如下：
# p(fail|treatment) and p(fail|control)
pftrt <- 0.5
pfctl <- 0.6
# generate the failure variable
fail <- c(rep("yes",pftrt*treatment.table[2]),
rep("no",(1-pftrt)*treatment.table[2]),
rep("yes",pfctl*treatment.table[1]),
rep("no",(1-pfctl)*treatment.table[1]))
# create the data frame.
all <- data.frame(treatment,fail)
# crosstable of failure by treatment.
ctable <- table(all$fail,all$treatment)
ctable
# count the number of cases failing in each group.
fc <- ctable[2,1]
ft <- ctable[2,2]
c(fc,ft)
# count the total number of cases in each group.
nc <- ctable[1,1]+ctable[2,1]
nt <- ctable[1,2]+ctable[2,2]
c(nc,nt)
# calculate the failure rates in each group.
pc <- fc/nc
pt <- ft/nt
c(pc,pt)
# calculate the classical treatment effect
# for the entire population.
cte_pop <- pt-pc
cte_pop
# now we move to the sample and power assessment.
# define the number of samples we will draw.
nsamples <- 100000
# and the sample size for each of the samples is
sampsize <- 500
# need an output vector for the cte and se(cte)
scte <- vector()
se_scte <- vector()
# this is a for() loop that will run nsamples times.
for(i in 1:nsamples)
{
c <- sample(
1:popsize,
size=sampsize,
replace=T
)
strt <- all$treatment[c]
sfail <- all$fail[c]
stable <- table(sfail,strt)
sfc <- stable[2,1]
sft <- stable[2,2]
snc <- stable[1,1]+stable[2,1]
snt <- stable[1,2]+stable[2,2]
spc <- sfc/snc
spt <- sft/snt
scte[i] <- spt-spc
spooledp <- (sfc+sft)/(snc+snt)
se_scte[i] <- sqrt(
spooledp*(1-spooledp)*
(1/snc+1/snt)
)
}
# let's look at some of the output
# from the final sample.
stable
spt
spc
scte[nsamples]
se_scte[nsamples]
# here is a chart summarizing the
# sampling distribution.
par(mfrow=c(1,2))
boxplot(scte)
hist(scte,prob=T,ylim=c(0,10))
lines(density(scte),lty=1,lwd=2)
# here are some descriptive statistics for the
# sampling distribution.
mean(scte)
sd(scte)
quantile(scte,0.025)
quantile(scte,0.975)
# this is a vector of z-statistics for testing
# the hypothesis that the population cte is 0
ztest <- scte/se_scte
# Identify critical region of z-distribution
qnorm(0.025)
qnorm(0.975)
# Decision rule
decision <- ifelse(
ztest<(-1.959964) | ztest>1.959964,
"reject Ho","fail to reject Ho"
)
# summarize the decisions
table(decision)
好，我的問題首先是我實在分辨不出什麼資訊在那篇文章裡是我需要輸入給R的
我現在只知道要餵給它sample size，可是在文章裡的table 1，我分辨不出
failure/success outcomes，導致我無法給pftrt和pfctl明確的rates
我自己有試著去找另一篇比較沒那麼複雜的experiment study，但不確定是否能用
https://goo.gl/fhr1r9
第二個問題是，當我要調整effect sizes到5%，10%，15%，我是要同時調整
qnorm和quantile嗎？
文章很長，我的頭很昏（昨晚才剛退燒），如果有高手願意指點的話就太好了，感恩
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 174.194.15.13
※ 文章網址: https://www.ptt.cc/bbs/Statistics/M.1508953561.A.327.html
→
10/26 09:57,
8年前
, 1F
10/26 09:57, 1F
→
10/26 20:52,
8年前
, 2F
10/26 20:52, 2F
推
02/08 01:38, , 3F
02/08 01:38, 3F