0

我今天提出的问题相当多。

我想计算两个数据帧的可变年龄平均值的置信区间（99% 水平，而不是 95） ，其中：infert_controlinfert_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，加上每个子集的标准差。

4

4 回答 4

4

您可以为此使用引导程序：

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
于 2014-01-19T18:33:52.180 回答
3

这么多答案...

随机样本的平均值具有 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

此外，分位数函数（eqqt(...)和qnorm(...)）默认返回下尾，因此除非您设置，否则您的限制将被反转lower.tail=F

于 2014-01-19T18:45:57.417 回答
2

您可以轻松地手动计算置信区间：

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（对不起维基百科链接，但它有一个很好的解释并向您展示了公式）

于 2014-01-19T18:19:11.263 回答
1

...或作为小功能：

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)
于 2014-01-19T18:27:02.387 回答