R 中比例的 95% 置信区间

机器算法验证 r 置信区间
2022-03-08 02:02:48

如何计算 95% 的区间来估计 R 城市中 SUV 的实际比例?我想计算这个数据的间隔:

vehicleType <- c("suv", "suv", "minivan", "car", "suv", "suv", "car", "car", "car", "car", "minivan", "car", "truck", "car", "car", "car", "car", "car", "car", "car", "minivan", "car", "suv", "minivan", "car", "minivan", "suv", "suv", "suv", "car", "suv", "car", "car", "suv", "truck", "truck", "minivan", "suv", "car", "truck", "suv", "suv", "car", "car", "car", "car", "suv", "car", "car", "car", "suv", "car", "car", "car", "truck", "car", "car", "suv", "suv", "minivan", "suv", "car", "car", "car", "car", "car", "minivan", "suv", "car", "car", "suv", "minivan", "car", "car", "car", "minivan", "minivan", "minivan", "car", "truck", "car", "car", "car", "suv", "suv", "suv", "car", "suv", "suv", "car", "suv", "car", "minivan", "car", "car", "car", "car", "car", "car", "car")

提前致谢。

3个回答

首先,请记住,比例的区间由下式给出:

p_hat +/- z * sqrt(p_hat * (1-p_hat)/n)

话虽如此,我们可以使用 R 来解决这样的公式:

# Set CI alpha level (1-alpha/2)*100%
alpha = 0.05

# Load Data
vehicleType = c("suv", "suv", "minivan", "car", "suv", "suv", "car", "car", "car", "car", "minivan", "car", "truck", "car", "car", "car", "car", "car", "car", "car", "minivan", "car", "suv", "minivan", "car", "minivan", "suv", "suv", "suv", "car", "suv", "car", "car", "suv", "truck", "truck", "minivan", "suv", "car", "truck", "suv", "suv", "car", "car", "car", "car", "suv", "car", "car", "car", "suv", "car", "car", "car", "truck", "car", "car", "suv", "suv", "minivan", "suv", "car", "car", "car", "car", "car", "minivan", "suv", "car", "car", "suv", "minivan", "car", "car", "car", "minivan", "minivan", "minivan", "car", "truck", "car", "car", "car", "suv", "suv", "suv", "car", "suv", "suv", "car", "suv", "car", "minivan", "car", "car", "car", "car", "car", "car", "car")

# Convert from string to factor
vehicleType = factor(vehicleType)

# Find the number of obs
n = length(vehicleType)

# Find number of obs per type
vtbreakdown = table(vehicleType)

# Get the proportion
p_hat = vtbreakdown['suv']/n

# Calculate the critical z-score
z = qnorm(1-alpha/2)

# Compute the CI
p_hat + c(-1,1)*z*sqrt(p_hat*(1-p_hat)/n)

所以,我们有:

0.1740293 0.3459707

对于p_hat

0.26

@Coatless 的方法将在大多数情况下完成工作(包括 OP 的情况)。但是,为了完整起见,我想我会添加几个其他选项。

引导方法

下面的函数n从数据向量中重新采样。对于每个重采样,它计算“成功”的比例,然后计算总体均值和 95% 置信区间

bp = function(x, lev, n = 1e3, alpha=0.05) {
  res = replicate(n, sum(sample(x, length(x), replace=TRUE) == lev)/length(x))
  return(list(mean=mean(res),
              `95% CI`=quantile(res, c(0.5*alpha,1-0.5*alpha))))
}

bp(vehicleType, "suv")

$mean
[1] 0.259628

$`95% CI`
  2.5% 97.5% 
  0.18  0.35 

binom包裹

binom包将在@Coatless 的答案中运行测试,假设错误是正态分布的。当“成功”的比例接近零或一和/或样本相对较小时,这可能导致不正确的值。binom.confintbinom包中还有其他选项可以避免这种陷阱。

在下面的输出中,asymptotic测试与@Coatless 编码的测试相同。例如,您可以使用参数获取其中一种方法的结果methods="exact"此外,binom.test()默认情况下使用精确 (Pearson-Klopper) 测试。

library(binom)

binom.confint(sum(vehicleType=="suv"), length(vehicleType))

          method  x   n      mean     lower     upper
1  agresti-coull 26 100 0.2600000 0.1836007 0.3541561
2     asymptotic 26 100 0.2600000 0.1740293 0.3459707
3          bayes 26 100 0.2623762 0.1788095 0.3485750
4        cloglog 26 100 0.2600000 0.1787357 0.3485852
5          exact 26 100 0.2600000 0.1773944 0.3573121
6          logit 26 100 0.2600000 0.1835016 0.3545416
7         probit 26 100 0.2600000 0.1818365 0.3526030
8        profile 26 100 0.2600000 0.1808127 0.3513344
9            lrt 26 100 0.2600000 0.1808329 0.3513338
10     prop.test 26 100 0.2600000 0.1797427 0.3590222
11        wilson 26 100 0.2600000 0.1840470 0.3537099

因为我们使用连续正态分布作为离散二项分布的近似值,所以应该在上述计算中添加一个校正项 (0.5/N):

在此处输入图像描述

请参阅此处了解更多详细信息