У меня есть RT записали дизайн внутри Ss с двумя факторами: номер блока и время суток. Меня попросили построить SD (RT) для каждого блока с планками ошибок (я выбрал CI).
Примеры данных:
dat <- tibble(
id = c(rep('A', 18), rep('B', 18), rep('C', 18)),
block = c(rep(c(1,2,3), 18)),
ToD = c(rep(c(rep('AM', 3), rep('PM', 3)), 9)),
rt = rnorm(54, 300, 50)
)
str(dat)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 54 obs. of 4 variables:
$ id : chr "A" "A" "A" "A" ...
$ block: num 1 2 3 1 2 3 1 2 3 1 ...
$ ToD : chr "AM" "AM" "AM" "PM" ...
$ rt : num 282 274 246 316 387 ...
В процессе раскопок я натолкнулся на этот ответ { ссылка } как средство вычисления CI. Адаптированный к моему случаю, этот метод будет
ci_func = function(x) {
quantile(replicate(1e4, sd(sample(x, rep = TRUE))), c(.025, .975))
}
aggregate(rt ~ block, dat, FUN = ci_func)
block rt.2.5% rt.97.5%
1 1 27.28086 48.86919
2 2 29.66846 52.65631
3 3 35.46145 68.64686
Но он генерирует CI вокруг невзвешенных оценок SD, что является довольно большой разницей (особенно в моих «реальных» данных):
# Block-wise, unweighted
aggregate(rt ~ block, dat, sd)
block rt
1 1 39.60575
2 2 43.35119
3 3 53.70013
# Vs. block-wise, weighted
dat %>%
group_by(id, block, ToD) %>%
summarize(v = var(rt)) %>%
group_by(block) %>%
summarize(sd = sqrt(mean(v)))
block SD
<fct> <dbl>
1 1 42.8
2 2 41.9
3 3 61.8
Если бы я хотел сгенерировать CI вокруг уязвимых SD, было бы это уместно?
ci_func <- function(x) {
reps <- replicate(
1e4,
sqrt( mean( sample( x, replace = TRUE )))
)
CI <- quantile(reps, c(0.025, 0.975))
return(CI)
}
tmp <- dat %>%
group_by(id, block, ToD) %>%
summarize(V = var(rt))
CIs <- tmp %>%
group_by(block) %>%
do(as.data.frame(ci_func(.$V)))
Если мы принимаем любой результат как признак успеха, то это работает!
head(CIs)
# A tibble: 6 x 2
# Groups: block [3]
block `ci_func(.$V)`
<fct> <dbl>
1 1 31.7
2 1 51.9
3 2 20.2
4 2 57.3
5 3 43.6
6 3 78.6
Но я не знаю, подходит ли это для этого. Мы будем очень благодарны за любые корректировки курса, которые могут быть предложены.