# 回归分析

R语言回归分析相关流程

### 数据预处理

*   中心化标准化

    ```
    data <- as.data.frame(scale(data, center = T, scale = T))    # 标准化
    data <- as.data.frame(scale(data, center = T, scale = F))    # 中心化
    ```

### 变量相关性检验

*   使用`Hmisc`

    ```
    library(Hmisc)
    rcorr(as.matrix(data), type = 'pearson')
    R <- rcorr(as.matrix(data), type = 'pearson')$r  # 相关系数阵
    P <- rcorr(as.matrix(data), type = 'pearson')$P  # 相关性检验的P值
    ```
*   使用`psych`

    ```
    library(psych)
    corr.test(data,use="complete")
    ```
*   热力图

    ```
    library(Hmisc)
    R <- rcorr(as.matrix(data), type = 'pearson')$r # 相关阵
    library(corrplot)
    corrplot(R, method="color", addCoef.col='grey', tl.col='black') # 热力图
    ```

    ![image-20221124123849236](https://euclid-picgo.oss-cn-shenzhen.aliyuncs.com/image/image-20221124123849236.png)
*   使用ggplot绘制热图

    ```
    library(Hmisc)
    R <- rcorr(as.matrix(data), type = 'pearson')$r # 相关阵
    library(reshape2)  
    melted_r <- melt(R)  # 将相关系数矩阵拉直
    library(ggplot2)
    ggplot(data = melted_r, aes(x=Var1,y=Var2,fill = value)) + # ggplot2 基本图层
      geom_raster() + # 绘制热力图
      scale_fill_gradient2(low = rgb(10,8,32,max = 255),mid = rgb(165,24,90,max = 255),
                           high = rgb(249,228,209,max = 255),limit = c(0,1),
                           name="Pearson\nCorrelation") + 
      theme( # 调整图表样式
        axis.title = element_blank(), # 删除ggplot2的坐标标签
        panel.background = element_blank(), # 删除ggplot2的灰色背景
        panel.grid.major = element_blank(), # 删除ggplot2的面板网格
        panel.border = element_blank(),
        axis.ticks = element_blank() # 删除ggplot2的轴刻度
      ) + 
      geom_text(aes(Var2,Var1,label = round(value,2)),
                color = '#228B22',size = 2) # 加上相关系数数值
    ```

    ![image-20221124124310441](https://euclid-picgo.oss-cn-shenzhen.aliyuncs.com/image/image-20221124124310441.png)

### 偏相关系数

在多要素所构成的系统中，当研究某一个要素对另一个要素的影响或相关程度时，把其他要素的影响视作常数（保持不变），即暂时不考虑其他要素影响，单独研究两个要素之间的相互关系的密切程度，所得数值结果为偏相关系数

*   计算偏相关系数并绘图

    ```
    library(corpcor)
    cor2pcor(R) # 计算得到偏相关阵
    library(corrplot)
    corrplot(cor2pcor(R), method="color", addCoef.col='grey', tl.col='black') # 热力图
    ```

    ![偏相关系数图](https://euclid-picgo.oss-cn-shenzhen.aliyuncs.com/image/202211261537897.png)

### 建立模型

*   最小二乘回归

    ```
    model <- lm(y~x1+x2+x3+x4+x5, data)     # 最小二乘回归
    #model <- lm(y~0+x1+x2+x3+x4+x5, data)   # 无截距项最小二乘回归
    fitted(model)                           # 拟合值
    model$coefficients                      # 回归系数
    model$residuals                         # 残差
    ```
*   计算杠杆值

    ```
    hii <- hatvalues(mo)                                # 杠杆值
    ```

    ![查看杠杆值(按照X大小顺序排序)](https://euclid-picgo.oss-cn-shenzhen.aliyuncs.com/image/image-20221124103158877.png)
*   计算模型残差

    ```
    SSE <- sum(mo$residuals^2)                          # 残差平方和
    sqrt(SSE/(length(data$x)-2))
    ```

### 回归模型检验

*   显著性检验

    ```
    summary(model)$r.squared         # R-square
    summary(model)$adj.r.squared     # Adjusted R-square
    sum((result - mean(result))^2)   # SSR
    sum((data$y - result)^2)         # SSE
    sum((data$y - mean(data$y))^2)   # SST, and SST = SSR + SSE
    plot(fitted(model), model$residuals, pch = 19, xlab = "fitted values",ylab = "residuals") # 残差图
    ```

    ![image-20221124124627132](https://euclid-picgo.oss-cn-shenzhen.aliyuncs.com/image/image-20221124124627132.png)

### 模型分析

*   点预测

    ```
    mo <- lm(y~x,data)                                                  # 最小二乘回归
    none_interval <- predict.lm(mo, data, interval="none")              # 计算点预测值
    ```
*   区间预测

    ```
    confidence_interval <- predict.lm(mo, data, interval="confidence")  # 计算置信区间预测值
    prediction_interval <- predict.lm(mo, data, interval="prediction")  # 计算预测区间预测值
    ```
*   绘制区间图

    ```
    # 对齐数据顺序
    data_order <- data[order(data$x),]
    confidence_interval_order <- confidence_interval[order(data$x),]
    prediction_interval_order <- prediction_interval[order(data$x),]
    ```

    开始绘图

    ```
    plot(data$x, data$y, xlab="x", ylab="y", pch = 16)  # 绘制轮廓
    polygon(c(data_order$x,rev(data_order$x)),c(prediction_interval_order[,2],rev(prediction_interval_order[,3])),col = rgb(221,234,243,max = 255),border = NA)
    polygon(c(data_order$x,rev(data_order$x)),c(confidence_interval_order[,2],rev(confidence_interval_order[,3])),col = "gray",border = NA)  # 其中polygon为绘制多边形的函数，rev为排序函数
    abline(mo, col="red", lwd=2)  
    points(data$x, data$y, pch = 16) 
    ```

    ![置信区间和预测区间](https://euclid-picgo.oss-cn-shenzhen.aliyuncs.com/image/image-20221124103249152.png)

### 预测相关

*   新值预测

    ```
    new <- data.frame(x = 3.5)  # 生成新的数据
    predict.lm(mo, new, interval="prediction")  # 点预测+预测上下界
    predict.lm(mo, new, interval="confidence")  # 点预测+置信上下界
    ```
