|
|
相关文章推荐
|
活泼的伤痕 · mysql中怎么存数组类型 - CSDN文库· 3 月前 · |
|
|
大力的充值卡 · 新版Edge浏览器兼容模式如何设置 - ...· 6 月前 · |
|
|
喝醉的米饭 · .Net 6 ...· 12 月前 · |
|
|
深沉的蚂蚁 · DataGrip常用快捷键 | 清澄秋爽· 1 年前 · |
evalCpp()
转换单一计算表达式
cppFunction()
转换简单的C++函数—Fibnacci例子
sourceCpp()
转换C++程序—正负交替迭代例子
sourceCpp()
转换C++源文件中的程序—正负交替迭代例子
sourceCpp()
转换C++源程序文件—卷积例子
wrap()
把C++变量返回到R中
as()
函数把R变量转换为C++类型
as()
和
wrap()
的隐含调用
//[[Rcpp::export]]
sourceCpp()
函数中直接包含C++源程序字符串
cppFunction()
函数中直接包含C++函数源程序字符串
evalCpp()
函数中直接包含C++源程序表达式字符串
depends
指定要链接的库
invisible
要求函数结果不自动显示
clone
函数
is_na
seq_along
seq_len
pmin
和
pmax
ifelse
sapply
和
lapply
sign
diff
kable()
函数制作表格
考虑连续型随机变量之间的关系。相关系数定义为 \[\rho(X, Y) = \frac{E\left[ (X - EX)(Y-EY) \right]} {\sqrt{\text{Var}(X)\text{Var}(Y)}} 又称Pearson相关系数。
\(-1 \leq \rho \leq 1\) 。 \(\rho\) 接近于 \(+1\) 表示 \(X\) 和 \(Y\) 有正向的相关； \(\rho\) 接近于 \(-1\) 表示 \(X\) 和 \(Y\) 有负向的相关。
相关系数代表的是 线性 相关性， 对于 \(X\) 和 \(Y\) 的其它相关可能反映不出来， 比如 \(X \sim \text{N}(0,1)\) ， \(Y = X^2\) ， 有 \(\rho(X, Y) = 0\) 。
给定样本 \((X_i, Y_i), i=1,2,\dots,n\) ，样本相关系数为 \[r = \frac{\sum_{i=1}^n (X_i - \bar X) (Y_i - \bar Y)} {\sqrt{\sum_{i=1}^n (X_i - \bar X)^2 \sum_{i=1}^n (Y_i - \bar Y)^2}}
用散点图和散点图矩阵直观地查看变量间的相关。
例如，线性相关的模拟数据的散点图:
图31.1: 线性相关数据二次曲线相关的模拟数据散点图:
图31.2: 二次曲线相关数据指数关系的例子:
图31.3: 指数关系相关数据对数关系的例子:
图31.4: 对数关系相关数据31.1.1 相关系数的性质
设变量
\(X\)
和
\(Y\)
的样本存放于R向量
x
和
y
中，
用
cor(x,y)
计算样本相关系数。
set.seed(1) x <- runif(30, 0, 10) xx <- seq(0, 10, length.out = 100) y <- 40 - (x-7)^2 + rnorm(30) yy <- 40 - (xx-7)^2 plot(x, y, pch=16) lines(xx, yy)
## [1] 0.8244374x <- runif(30, 0, 10) xx <- seq(0, 10, length.out = 100) y <- 40 - (x-5)^2 + rnorm(30) yy <- 40 - (xx-5)^2 plot(x, y, pch=16) lines(xx, yy)
## [1] -0.042684x <- runif(30, 0, 10) xx <- seq(0, 10, length.out = 100) y <- 40*exp(-x/2) + rnorm(30) yy <- 40*exp(-xx/2) plot(x, y, pch=16) lines(xx, yy)
## [1] -0.884111731.1.2 相关与因果
相关系数是从数据中得到的 \(x\) 和 \(y\) 两个变量关系的一种描述， 不能直接引申为因果关系。
例如，增加身高，不一定增长体重。
实际中有许多错误理解相关与因果性的例子。 例如：研究发现喝咖啡的人群比不喝咖啡的人群心脏病发病率高。 于是断言喝咖啡导致心脏病危险增加。 进一步研究更多的因素发现， 喝咖啡放糖很多的人心脏病发病率才增高了。
31.1.3 相关系数大小
31.1.4 相关系数的检验
相关系数 \(r\) 是总体相关系数 \(\rho\) 的估计。 \[H_0: \rho=0 \longleftrightarrow H_a: \rho\neq 0\] 检验统计量 \[t = \frac{r \sqrt{n-2}}{\sqrt{1 - r^2}}\] 当总体 \((X, Y)\) 服从二元正态分布且 \(H_0: \rho=0\) 成立时， \(t\) 服从 \(t(n-2)\) 分布。 设 \(t\) 统计量值为 \(t_0\) ，p值为 \[P(|t(n-2)| > |t_0|)\]
在R中用
cor.test(x,y)
计算检验。
得身高与体重的相关系数为0.88，检验的p值为
7.887e-07
，
95%置信区间为
\([0.70, 0.95]\)
。
可以用ggstatsplot包将相关系数检验进行图示：
## You can cite this package as:
## Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
## Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167ggscatterstats( data = d.class, x = height, y = weight, xlab = "身高", ylab = "体重", title = "体重与身高的相关性"## Registered S3 method overwritten by 'ggside': ## method from ## +.gg ggplot2## Pearson's product-moment correlation ## data: x and y ## t = 0.93076, df = 28, p-value = 0.3599 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## -0.1994815 0.5021658 ## sample estimates: ## cor ## 0.1732378## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
得x与y的样本相关系数为0.17, p值为0.36，不显著。
31.1.4.1 相关系数检验的功效和样本量
pwr包的
pwr.r.test()
可以计算相关系数检验的功效和样本量。
用相关系数绝对值作为效应大小，
常用值为：
比如，样本量为30， 检验“大”的效应， 计算其功效：
## Warning: package 'pwr' was built under R version 4.2.2pwr.r.test( r = 0.5, n = 30, alternative = "two.sided" ## approximate correlation power calculation (arctangh transformation) ## n = 30 ## r = 0.5 ## sig.level = 0.05 ## power = 0.8250298 ## alternative = two.sided
功效为 \(83\%\) 。
检验“小”的效应， 功效为 \(80\%\) ， 需要的样本量：
## Conventional effect size from Cohen (1982) ## test = r ## size = small ## effect.size = 0.1pwr.r.test( r = 0.1, power = 0.80, alternative = "two.sided" ## approximate correlation power calculation (arctangh transformation) ## n = 781.7516 ## r = 0.1 ## sig.level = 0.05 ## power = 0.8 ## alternative = two.sided
需要的样本量为 \(782\) 。
31.1.5 相关阵
如果
x
是一个仅包含数值型列的数据框，
则
cor(x)
计算样本相关系数矩阵；
var(x)
计算样本协方差阵。
corrgram::corrgram()
可以绘制相关系数矩阵的图形，
用颜色和阴影浓度代表相关系数正负和绝对值大小，
用饼图中阴影部分大小代表相关系数绝对值大小。
library(corrgram) corrgram( baseball[,c("Assists","Atbat","Errors","Hits","Homer","logSal", "Putouts","RBI","Runs","Walks","Years")], order=TRUE, main="Baseball data PC2/PC1 order", lower.panel=panel.shade, upper.panel=panel.pie)
其中选项
order=TRUE
重排各列的次序使得较接近的列排列在相邻位置。
ggcorrplot包也提供了类似的功能：
library(ggcorrplot) baseball |> select(all_of(c( "Assists","Atbat","Errors","Hits","Homer","logSal", "Putouts","RBI","Runs","Walks","Years"))) |> cor(use="pairwise.complete.obs") |> ggcorrplot(hc.order=TRUE)
ggstatsplot包也提供了这样的功能， 并且进行检验：
library(ggstatsplot) ggcorrmat( data = baseball[,c("Assists","Atbat","Errors","Hits","Homer","logSal", "Putouts","RBI","Runs","Walks","Years")] 图31.10: 相关阵图形（ggstatsplot）## Spearman's rank correlation rho ## data: d.class$height and d.class$weight ## S = 164.43, p-value = 2.983e-06 ## alternative hypothesis: true rho is not equal to 0 ## sample estimates: ## rho ## 0.8557611## Warning in cor.test.default(d.class$height, d.class$weight, method = ## "spearman"): Cannot compute exact p-value with ties
秩相关系数值为 \(0.8558\) 。
31.1.7 Kendall tau系数
当变量 \((X, Y)\) 正相关性很强时， 任意两个观测的 \(X\) 值的大小顺序应该与 \(Y\) 值的大小顺序相同； 如果独立， 一对观测的 \(X\) 值比较和 \(Y\) 值比较顺序相同与顺序相反的数目应该基本相同。 Kandall tau系数也是取值于 \([-1, 1]\) 区间， 用这样的思想表示两个变量的相关性和正负。
## Warning in cor.test.default(d.class$height, d.class$weight, method = "kendall"):
## Cannot compute exact p-value with tiesKendall tau系数的值为 \(0.7143\) 。
31.2 一元回归分析
31.2.1 模型
考虑两个变量 \(Y\) 与 \(X\) 的关系，希望用 \(X\) 值的变化解释 \(Y\) 值的变化。 \(X\) 称为自变量(independent variable)， \(Y\) 称为因变量(response variable)。
\[\begin{aligned} Y = a + b X + \varepsilon, \quad \varepsilon \sim \text{N}(0,\sigma^2) \end{aligned}在回归分析理论中通常假定 \(X\) 是非随机的， 因变量 \(Y\) 的随机性来自误差项 \(\varepsilon\) 。 实际应用中 \(X\) 常常也是随机的。
设观测值为 \((X_i, Y_i), i=1,2,\dots,n\) ，假设观测值满足上模型。 \[\begin{aligned} Y_i = a + b X_i + \varepsilon_i, \quad \varepsilon_i \text{ iid }\sim \text{N}(0,\sigma^2) . \end{aligned}
31.2.2 最小二乘法
直观上看，要找最优的直线 \(y=a+bx\) 使得直线与观测到的点最接近。
set.seed(1) nsamp <- 30 x <- runif(nsamp, -10, 10) y <- 20 + 0.5*x + rnorm(nsamp,0,0.5) plot(x, y) abline(lm(y ~ x), col="red", lwd=2)
\(a,b\) 的解用 最小二乘法 得到:
\[\min_{a,b} \sum_{i=1}^n \left( y_i - a - b x_i \right)^2\]
最小二乘解表达式为
\[\begin{aligned} \hat b =& \frac{\sum_i (x_i - \bar x)(y_i - \bar y)} {\sum_i (x_i - \bar x)^2} = r_{xy} \frac{S_y}{S_x} \\ \hat a =& \bar y - \hat b \bar x \end{aligned}\] 其中 \(r_{xy}\) 为 \(X\) 与 \(Y\) 的样本相关系数， \(S_x\) 与 \(S_y\) 分别为 \(X\) 和 \(Y\) 的样本标准差。
随机误差方差 \(\sigma^2\) 的估计取为
\hat\sigma^2 = \frac{1}{n-2}\sum_i (y_i - \hat a - \hat b x_i)^2标准误差: 为衡量 \(\hat a\) 和 \(\hat b\) 的估计精度， 计算 \(\text{SE}(\hat a)\) 和 \(\text{SE}(\hat b)\) 。
估计结果:
31.2.3 回归有效性
31.2.3.1 拟合优度指标
按照最小二乘估计公式，只要 \(x_1, x_2, \dots, x_n\) 不全相等， 不论 \(x, y\) 之间有没有相关关系， 都能计算出参数估计值。 得到的回归直线与观测数据的散点之间的接近程度代表了回归结果的优劣。
残差平方和 \[\text{SSE}=\sum_{i=1}^n (y_i - \hat a - \hat b x_i)^2\] 残差平方和越小， 说明回归直线与观测数据点吻合得越好。
\(y\) 是因变量， 因变量的数据的离差平方和 \[\text{SST}=\sum_{i=1}^n (y_i - \bar y)^2\] 代表了因变量的未知变动大小， 需要对这样的变动用自变量进行解释， 称SST为总平方和。
可以证明如下平方和分解公式： \[\text{SST} = \text{SSR} + \text{SSE}\] 其中SSR称为回归平方和。
令 \(\hat y_i = \hat\beta_0 + \hat\beta_1 x_i\) 称为第 \(i\) 个拟合值， 回归平方和为 \[\text{SSR}=\sum_{i=1}^n (\hat y_i - \bar y)^2 = \hat b^2 \sum_{i=1}^n (x_i - \bar x)^2
在总平方和中， 回归平方和是能够用回归斜率与自变量的变动解释的部分， 残差平方和是自变量不能解释的变动。 分解中，回归平方和越大，误差平方和越小， 拟合越好。 \[R^2 = \frac{\text{SSR}}{\text{SST}} = 1 - \frac{\text{SSE}}{\text{SST}}\] 称为回归的复相关系数，或判定系数。
\(0 \leq R^2 \leq 1\) 。判定系数越大，回归拟合越好。 \(R^2=1\) 时，观测点完全落在一条直线上面。 一元回归中 \(R^2\) 是 \(x\) 和 \(y\) 的样本相关系数的平方。 在多元回归时只能用复相关系数。
估计误差项方差 \(\sigma^2\) 为 \hat\sigma^2 = \frac{1}{n-2}\sum_{i=1}^n (y_i - \hat a - \hat b x_i)^2 = \frac{1}{n-2} \text{SSE}
\(\hat\sigma\) 是 \(\sigma\) 的估计， R的回归结果中显示为“residual standard error”（残差标准误差）。 \(\hat\sigma\) 是残差 \(e_i = y_i - \hat y_i\) 的标准差的一个粗略估计。
31.2.3.2 线性关系显著性检验
当 \(b=0\) 时，模型退化为 \(Y = a + \varepsilon\) ， \(X\) 不出现在模型中， 说明 \(Y\) 与 \(X\) 不相关。检验 H_0: b=0 \longleftrightarrow H_\text{a}: b \neq 0 \[F = \frac{\text{SSR}}{\text{SSE}/(n-2)}\] 检验p值为 \(P(F(1, n-2) > c)\) (设 \(c\) 为 \(F\) 统计量的值)。 取检验水平 \(\alpha\) , p值小于等于 \(\alpha\) 时拒绝 \(H_0\) ， 认为 \(y\) 与 \(x\) 有显著的线性相关关系， 否则认为 \(y\) 与 \(x\) 没有显著的线性相关关系。
也可以使用t统计量 \[t = \frac{\hat b}{\text{SE}(\hat b)}\] \text{SE}(\hat b) = \hat\sigma / \sqrt{\sum_{i=1}^n (x_i - \bar x)^2} 设t统计量的值为 \(c\) ， 检验p值为 \(P(|t(n-2)| > |c|)\) 。
31.2.4 R程序
设数据保存在数据框d中，变量名为y和x，用R的
lm()
函数计算回归，如:
set.seed(1) nsamp <- 30 x <- runif(nsamp, -10, 10) y <- 20 + 0.5*x + rnorm(nsamp,0,0.5) d <- data.frame(x=x, y=y) lm1 <- lm(y ~ x, data=d); lm1
结果只有回归系数。需要用
summary()
显示较详细的结果。如
结果中
\(H_0: b=0\)
的检验结果p值为
\(<2e-16\)
。
这个结果可以从输出最后一行的F统计量检验结果看到，
也可以从稀疏估计中自变量
x
所在行的t检验结果看到。
又如对d.class数据集，建立体重对身高的回归方程:
## Call: ## lm(formula = weight ~ height, data = d.class) ## Residuals: ## Min 1Q Median 3Q Max ## -17.6807 -6.0642 0.5115 9.2846 18.3698 ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -143.0269 32.2746 -4.432 0.000366 *** ## height 3.8990 0.5161 7.555 7.89e-07 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Residual standard error: 11.23 on 17 degrees of freedom ## Multiple R-squared: 0.7705, Adjusted R-squared: 0.757 ## F-statistic: 57.08 on 1 and 17 DF, p-value: 7.887e-07
身高对体重的影响是显著的(p值
7.89e-7
)。
31.2.5 回归诊断
除了系数的显著性检验以外， 对误差项是否符合回归模型假定中的独立性、方差齐性、正态性等， 可以利用回归残差进行一系列的回归诊断。 还可以计算一些异常值、强影响点的诊断。
对 \(n\) 组观测值，回归模型为 y_i = a + b x_i + \varepsilon_i \hat y_i = \hat a + \hat b x_i 残差（residual）为 \[e_i = y_i - \hat y_i\]
残差相当于模型中的随机误差， 但仅在模型假定成立时才是随机误差的合理估计。 所以残差能够反映违反模型假定的情况。
\(e_i\) 有量纲，不好比较大小。 定义标准化残差（内部学生化残差）： r_i = \frac{e_i}{S_e\sqrt{1 - \frac{1}{n} - \frac{(x_i - \bar x)^2}{\sum_{k=1}^n (x_k - \bar x)^2}}} 样本量较大时标准化残差近似服从标准正态分布， 于是取值于 \([-2,2]\) 范围外的点是可疑的离群点。
R中用
residuals()
从回归结果计算残差，
用
rstandard()
从回归结果计算标准化残差。
R中
plot(lmres)
(lmres为回归结果变量)可以做简单回归诊断。
第一幅图是残差对预测值散点图， 散点应该随机在0线上下波动，不应该有曲线模式、分散程度增大模式、 特别突出的离群点等情况。
第二幅图是残差的正态QQ图， 散点接近于直线时可以认为模型误差项的正态分布假定是合理的。
第三幅图是误差大小(标准化残差绝对值的平方根)对拟合值的图形， 可以判断方差齐性假设(方差 \(\sigma^2\) 不变)。
第四个幅是残差对杠杆量图，并叠加了Cook距离等值线。 杠杆量代表了回归自变量对结果的影响大小， 超过 \(4/n\) 的值是需要重视的。 Cook距离考察删去第 \(i\) 个观测对回归结果的影响， 如果有超出两条虚线范围之外的点， 就可能是强影响点。
显示残差与标准化残差：
标准化残差 x <- runif ( 30 , 0 , 10 ) xx <- seq ( 0 , 10 , length.out = 100 ) y <- 40 - (x -7 ) ^ 2 + rnorm ( 30 ) yy <- 40 - (xx -7 ) ^ 2 lms1 <- lm (y ~ x) opar <- par ( mfrow= c ( 1 , 2 )) plot (x, y, pch= 16 , main= "数据和真实模型" ) lines (xx, yy) plot (x, rstandard (lms1), main= "使用线性回归的残差" )x <- sort(runif(30, 0, 10)) y <- 10 + 2*x + rnorm(30)*(seq(30)/5) lms2 <- lm(y ~ x) opar <- par(mfrow=c(1,2)) plot(x, y, pch=16, main="数据和真实模型") abline(a=10, b=2) plot(x, rstandard(lms2), main="线性模型的残差")
序列自相关：
当数据随时间变化时，还需要考虑前后是否有序列自相关。 自相关残差图形例子：
ar1 <- arima.sim(list(ar=c(0.5)), n=30) plot(1:30, ar1[], type="l", xlab="time", ylab="residual", main="有序列自相关的残差图形")
图中横坐标是观测序号。
car::ncvTest()
检验方差齐性。零假设是方差齐性成立。
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.664307, Df = 1, p = 0.19702
用Durbin-Watson检验（DW检验）可以检验残差中是否有序列自相关，
零假设是没有序列自相关。
R中用
car::durbinWatsonTest()
检验。
序列自相关检验仅当数据中数据是沿时间等间隔记录时有意义。
## lag Autocorrelation D-W Statistic p-value
## 1 0.1959139 1.57293 0.24
## Alternative hypothesis: rho != 0也可以对回归残差绘制ACF图， 如果除了横坐标0之外都落在两条水平界限内则认为没有序列自相关， 如果有明显超出界限的就认为有序列自相关。
31.2.6 预测区间
设
\(X\)
取
\(x_0\)
，
\(Y\)
的预测值为
\(\hat y_0 = \hat a + \hat b x_0\)
。
置信水平为
\(1-\alpha\)
的预测区间为
\[\hat y_0 \pm \lambda \hat\sigma
\sqrt{1 + \frac{1}{n}
+ \frac{(\bar x - x_0)^2}{\sum_i (x_i - \bar x)^2}}
其含义是
\(P(y_0 \in \text{预测区间}) = 1 - \alpha\)
。
\(\lambda\)
是标准正态分布双侧
\(\alpha\)
分位数
qnorm(1 - alpha/2)
。
用
predict(lmres)
得到
\(\hat y_i, i=1,\dots,n\)
的值，
lmres
表示回归结果变量。
用
predict(lmres, interval="prediction")
同时得到预测的置信区间，
需要的话加入
level=
选项设定置信度。
回归的置信区间有两种， 上述的区间是“预测区间”(CLI)， 使得 \(P(y_0 \in \text{预测区间}) = 1 - \alpha\) ; 另一种区间称为均值的置信区间(CLM)， 使得 \(P(Ey_0 \in \text{均值置信区间}) = 1 - \alpha\) 。
31.2.7 控制
如果需要把 \(Y\) 的值控制在 \([y_l, y_u]\) 范围内，问如何控制 \(X\) 的范围， 可以求解 \(x_0\) 的范围使上面的置信区间包含在 \([y_l, y_u]\) 内。
近似地可以解不等式
\[\begin{aligned}
\hat a + \hat b x_0 - z_{1-\frac{\alpha}{2}} \hat \sigma \geq & y_l \\
\hat a + \hat b x_0 + z_{1-\frac{\alpha}{2}} \hat \sigma \leq & y_u
\end{aligned}\]
其中
\(z_{1-\frac{\alpha}{2}}\)
为标准正态分布双侧
\(\alpha\)
分位数
（用
qnorm(1-alpha/2)
计算）。