作者Archangels (Geh Zum Teufel!)
看板Statistics
标题[问题] 使用R计算effective size
时间Thu Oct 26 01:45:55 2017
如果是跟统计软体有关请重发文章,使用程式做为分类。
统计软体,如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://webptt.com/cn.aspx?n=bbs/Statistics/M.1508953561.A.327.html
1F:→ andrew43: 你可以先处理最後算power的code。 10/26 09:57
2F:→ jupit: 作业的话,建议可以去soho或part-time版找协助会比较适合 10/26 20:52
3F:推 ppm2013: R shiny 02/08 01:38