|
“ No one konws everything, and you don't have to.” --free傻孩子本期为大家推荐的是计算数据（变量）的置信区间。置信区间是对整体样本分布的估算区间。当置信水平一定时（如95%），置信区间随检测样本量的增多而变窄（更准确）；当样本量一定时，置信区间随置信水平增加变宽。置信区间的目标是以少量的样本来评估整体对象的分布区域。那么如何使用R来实现这一功能呢？ 
当已经获得了一组个数有限的数据时，整体数据集的分布应如何评估呢？如上图所示，整体数据的分布范围在[mean-σ，mean+σ]的可能性在68%；分布范围在[mean-2σ，mean+2σ]的可能性在95%；而整体数据的分布范围在[mean-3σ，mean+3σ]的可能性是99.7%。 前面我们已经简单地了解了置信区间的概念和置信区间的作用，在不同情况下如何计算置信区间呢？条件和公式如下：参数解读：μ整体样本的平均值；σ标准差；x拔抽样样本的均值；n样本数量；S抽样样本的标准差；Z 对应置信区间的查表数据；t 对应对应置信区间的另一种情况的查表数据。第一种情况：当在已知整体数据的σ（标准差）时计算抽取样本的置信区间。这种情况一般发生的可能比较小。第二种情况：当整体样本的σ未知时。这种情况在现实生活和实验中发生的最多。即我们检测的样本数量有限，需要评估对象分布的大致范围。
第三种情况：当已知整体数据的σ，但是均值未知时，这种情况以本人目前的认知来看，发生的情况比较少。所以不做讨论，也没在接下来的函数中编辑，如果需要以后会补上。mean_ci <- function(x,sigma = NULL, conf = 0.95) {
if(is.null(sigma)) {
se = sd(x)/sqrt(length(x))
alpha = 1- conf
x_mean = mean(x)
lower = x_mean - se*qt(1-alpha/2,(length(x)-1))
upper = x_mean + se*qt(1-alpha/2,(length(x)-1))
} else {
se = sigma/sqrt(length(x))
alpha = 1- conf
x_mean = mean(x)
lower = x_mean - se*qnorm(1-alpha/2)
upper = x_mean + se*qnorm(1-alpha/2)
}
data.frame(x_mean,lower,upper)
}
# x为抽取一定数量的样本
#sigma(σ)为标准差，默认为未知
#conf为置信区间，默认为95%；可以设置为97%，也可以为99.7%如何解读这个自己编写的小函数：
主要用到了IF{ } else{ }函数；首先根据sigma判断是否给出了标准差（标准差是否已知），若未知成立（TRUE），则运行if内的代码；否者运行else内的代码。“/” 表述除以；sqrt()代表对括号内的数据开方；length()；计算数据的个数也就是n；“ - ” 表示减去；mean(); 求均值。qt()和qnorm()；分别对应σ未知和已知时的“查表表格”。 使用数据“mtcars”里的变量试验函数运行情况： data("mtcars")
x <- mtcars$mpg
x
# 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 #16.4 17.3 15.2 10.4 10.4 14.7
#32.4 30.4 33.9 21.5 15.5 15.2 13.3 19.2 27.3 26.0 30.4 15.8 #19.7 15.0 21.4
mean_ci(x,conf = 0.95)
# x_mean lower upper
# 20.09062 17.91768 22.26357
#第一个是均值，第二个和第三个分别是95%置信区间的下限和上限
mean_ci(x,conf = 0.99)
#x_mean lower upper
#20.09062 17.16706 23.01419
#若sigma已知为2
mean_ci(x,sigma = 2,conf = 0.95)
# x_mean lower upper
# 20.09062 19.39767 20.78358前面三节我们讲述了数据（单变量）的置信区间的计算方法和函数，那么模型的置信区间应该如何计算呢？模型指包含一个或多个自变量和一个因变量。如 y = ax2 +bx +c
这里就需要用到R里自带的基础函数predict了。
我们依然以mtcars数据为例。 假设已知disp与mpg满足一元二次函数，disp为因变量，mpg为自变量。 在R里两者构成的模型应该为： disp
~ I(mpg^2) + mpg fit <- lm(disp ~ I(mpg^2) + mpg, data = mtcars)
pre <- data.frame(predict(fit, mtcars, interval = "confidence", level = 0.95))
pre
#以95%为例
变量fit为根据模型预测得到的所有disp的值（与检测到的disp不一定完全一致）；lwr为95%置信区间的下限，upr为95%置信区间的上限。
绘图如下： library(tidyverse)
pre$mpg <- mtcars$mpg
pre$disp <- mtcars$disp
ggplot(data = pre)+
geom_point(aes(x=mpg, y = disp))+
geom_line(aes(x=mpg, y = fit)) +
geom_line(aes(x=mpg, y = lwr)) +
geom_line(aes(x=mpg, y = upr))
|
