Вычислить CI вокруг SD, полученных из заблокированного внутрисубъектного дизайна - PullRequest
1 голос
/ 26 мая 2020

У меня есть 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

Но я не знаю, подходит ли это для этого. Мы будем очень благодарны за любые корректировки курса, которые могут быть предложены.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...