Расчет доверительных интервалов для двух наборов данных

Довольно много вопросов, которые я задал сегодня.

Я хотел бы рассчитать доверительный интервал (уровень 99%, а не 95) для среднего значения переменной age двух фреймов данных, infert_control и infert_patient, где:

infert_control = subset(infert$age, infert$case == 0)
infert_patient = subset(infert$age, infert$case == 1)

infert - это встроенный набор данных R, для тех, кто с ним не знаком, вот он: случай 0 обозначает пациентов контрольной группы, случай 1 - фактических.

> infert
    education age parity induced case spontaneous stratum pooled.stratum
1      0-5yrs  26      6       1    1           2       1              3
2      0-5yrs  42      1       1    1           0       2              1
3      0-5yrs  39      6       2    1           0       3              4
4      0-5yrs  34      4       2    1           0       4              2
5     6-11yrs  35      3       1    1           1       5             32
6     6-11yrs  36      4       2    1           1       6             36
7     6-11yrs  23      1       0    1           0       7              6
8     6-11yrs  32      2       0    1           0       8             22
9     6-11yrs  21      1       0    1           1       9              5
10    6-11yrs  28      2       0    1           0      10             19
11    6-11yrs  29      2       1    1           0      11             20
...
239   12+ yrs  38      6       0    0           2      74             63
240   12+ yrs  26      2       1    0           1      75             49
241   12+ yrs  31      1       1    0           0      76             45
242   12+ yrs  31      2       0    0           1      77             53
243   12+ yrs  25      1       0    0           1      78             41
244   12+ yrs  31      1       0    0           1      79             45
245   12+ yrs  34      1       0    0           0      80             47
246   12+ yrs  35      2       2    0           0      81             54
247   12+ yrs  29      1       0    0           1      82             43
248   12+ yrs  23      1       0    0           1      83             40

Как правильно решить эту проблему?

Я уже вычислил Среднее значение столбца age для infert_control и infert_patient, плюс стандартное отклонение каждого подмножества.


r dataset statistics confidence-interval
person Dimitris Sfounis    schedule 19.01.2014    source источник
comment
это вопрос статистики. Но вы можете сообщить о правильных quantiles   -  person rawr    schedule 19.01.2014
comment
Может быть, t.test сможет вам помочь?   -  person Fernando    schedule 19.01.2014
comment
@Fernando Поскольку t.test() требуется во второй части этого упражнения, я не думаю, что я должен использовать его для этой части, которая является первой. Все это - вводное задание.   -  person Dimitris Sfounis    schedule 19.01.2014
comment
Понятно, может быть, вы можете проверить? T.test, чтобы увидеть детали реализации.   -  person Fernando    schedule 19.01.2014
comment
@Fernando Я не верю, что t.test() имеет отношение к моему вопросу.   -  person Dimitris Sfounis    schedule 19.01.2014
comment
Обратите внимание, что infert_control и infert_patient являются числовыми векторами, а не кадрами данных. Также вам стоит взглянуть на ?subset. Вы не должны использовать $ внутри subset.   -  person Roland    schedule 19.01.2014

Ответы (4)


arrow_upward
2
arrow_downward

Вы можете легко рассчитать доверительный интервал вручную:

infert_control <- subset(infert$age, infert$case == 0)

# calculate needed values
m <- mean(infert_control)
s <- sd(infert_control)
n <- length(infert_control)

# calculate error for normal distribution (choose you distribution here, e.g. qt for t-distribution)
a <- 0.995 # 99% CI => 0.5% on both sides
error <- qnorm(a)*s/sqrt(n)

# calculate CI
ci_lower <- m-error
ci_upper <- m+error

См. Также http://en.wikipedia.org/wiki/Confidence_interval (извините за ссылку на википедию , но он имеет хорошее объяснение и показывает вам формулу)

person sgibb    schedule 19.01.2014
comment
Ваша помощь очень ценится. Позвольте мне изучить ваш ответ. - person Dimitris Sfounis; 19.01.2014

arrow_upward
4
arrow_downward

Вы можете использовать для этого бутстрап:

library(boot)
set.seed(42)
boot_mean <- boot(infert_control, function(x, i) mean(x[i]), R=1e4)
quantile(boot_mean$t, probs=c(0.005, 0.995))
#      0.5%    99.5% 
#  30.47273 32.58182

Или, если вы не хотите использовать библиотеку:

set.seed(42)
R <- 1e4
boot_mean <- colMeans(
                matrix(
                   sample(infert_control, R * length(infert_control), TRUE), 
                   ncol=R))
quantile(boot_mean, probs=c(0.005, 0.995))
#    0.5%    99.5% 
#30.42424 32.55152
person Roland    schedule 19.01.2014
comment
Прошу прощения, но я не знаком с библиотекой начальной загрузки. - person Dimitris Sfounis; 20.01.2014
comment
Что ж, прочтите документацию или сделайте это вручную (см. Мою правку). - person Roland; 20.01.2014

arrow_upward
3
arrow_downward

Так много ответов ...

Среднее значение случайной выборки имеет t-распределение, а не нормальное, хотя t -> N как df -> Inf.

cl <- function(data,p) {
  n  <- length(data)
  cl <- qt(p/2,n-1,lower.tail=F)*sd(data)/sqrt(n)
  m  <- mean(data)
  return(c(lower=m-cl,upper=m+cl))
}
cl.control <- cl(infert_control,0.01)
cl.control
#    lower    upper 
# 30.42493 32.55689 

cl.patient <- cl(infert_patient,0.01)
cl.patient
#    lower    upper 
# 30.00221 33.05803

aggregate(age~case,data=infert,cl,p=0.01)  # much better way...
#   case age.lower age.upper
# 1    0  30.42493  32.55689
# 2    1  30.00221  33.05803

Кроме того, функции квантилей (например, qt(...) и qnorm(...)) по умолчанию возвращают нижний хвост, поэтому ваши пределы будут отменены, если вы не установите lower.tail=F

person jlhoward    schedule 19.01.2014
comment
Это очень важно, @jhoward. Похоже, это поможет мне со второй частью упражнения. - person Dimitris Sfounis; 19.01.2014

arrow_upward
1
arrow_downward

... или как небольшая функция:

cifun <- function(data, ALPHA){
  c(mean(data) - qnorm(1-ALPHA/2) * sd(data)/sqrt(length(data)),
    mean(data) + qnorm(1-ALPHA/2) * sd(data)/sqrt(length(data)))
}

cifun(infert_control, 0.01)
person holzben    schedule 19.01.2014