R建模范例分析

Regression analysis is widely used in statisitic, always analysis the correlation among different variables.

这是一个非常简单资料处理和建模的例子，主要是为了熟悉整个资料分析和建模的流程。其中有很多细节可以继续展开。

读取资料

大部分读取资料，读取excel档一本需要转成CSV档来读取。
这里我们以 iris.csv 这个 CSV 档案作为范例：

1 2	# 读取 iris.csv my.iris.df <- read.csv("iris.csv")

如果自己的csv档没有放在目前的目录中，可以使用绝对路径来指定：

1 2	# 使用绝对路径 my.iris.df <- read.csv("D:\\iris.csv")

双反斜线（\）在R中属于特殊的跳脱字元，所以在撰写绝对路径的时候，凡是要输入反斜线的地方，都要改为双反斜线。

读取资料完成后，为了方便查看大体量的资料可以用head（）来查看前面几行

1 2	# 查看 my.iris.df 的资料的前面20rows head(my.iris.df, 20)

将csv表格读进R之后个栏位若有空白或是特殊字元，会被自动替代为（.），它不会影响到资料，只是在变数指定的时候要使用新的名字

资料的预处理

读进去之后我们就需要更详细的了解资料的内部结构：

1 2	# 查看 my.iris.df 内部结构 str(my.iris.df)

str输出的一个row就是一个变数的栏位，标示了变数名称与类型，后面接的是实际的资料笔数。

检查完资料的类型没有什么问题之后，我们需要处理资料中残留的一些不干净的数据以及缺失值NA。

1 2	# 检查是否有 NA 的资料 my.iris.df[!complete.cases(my.iris.df),]

如果资料是完整的就会返回一个空的dataframe

分析资料与绘图

通常在开始分析资料之前，会先用 summary 看一下各个变数的基本统计量：

1 2	# 查看基本统计量 summary(my.iris.df)

当然光看数字肯定是非常不直观的，绘图会是一个查看各组数据相关性的好方法。这里我们示范使用** ggplot2 与 GGally** 套件来画图的方法。

1
2
3

# 载入 ggplot2 与 GGally 套件
library(ggplot2)
library(GGally)

使用 ggpairs 可画出很漂亮的资料分布矩阵图：

1 2	# 绘制资料分布矩阵图 ggpairs(my.iris.df)

如果只想画出简单的 XY 散布图：

# 画出 XY 散布图
require(ggplot2)
qplot(x = Petal.Length,
      y = Petal.Width,
      data = my.iris.df)

如果想要比较不同的 Species 类别的 XY 散布图，可以用color区分出不同的species：

# 画出 XY 散布图，依据 Species 上色
require(ggplot2)
qplot(x = Petal.Length,
      y = Petal.Width,
      data = my.iris.df,
      color = Species)

建立回归模型

建立回归模型可以使用 lm 这个函数（lm 代表 linear models），假设我们想要拿 Sepal.Length 作为反应变数（也就是 Y），而 Sepal.Width、Petal.Length 与 Petal.Width 作为解释变数（也就是 X），则可以执行：

1
2
3

# 建立回归模型
iris.lm <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
              data = my.iris.df)

其中第一个参数放的就是所谓的公式（formula），它用来表示回归模型的一种表示法，中间的** ~ 相当于回归公式的等号**，左边的 Sepal.Length 就是 Y，而右边放的三个变数则是 X。第二个参数 data 则是用来指定资料来源的 data frame。

建立好回归模型之后，可以使用 summary 查看模型配适的结果：

1 2	# 查看模型配适结果 summary(iris.lm)

结果为：

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
    data = my.iris.df)
Residuals:
     Min       1Q   Median       3Q      Max
-0.82816 -0.21989  0.01875  0.19709  0.84570
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.85600    0.25078   7.401 9.85e-12 ***
Sepal.Width   0.65084    0.06665   9.765  < 2e-16 ***
Petal.Length  0.70913    0.05672  12.502  < 2e-16 ***
Petal.Width  -0.55648    0.12755  -4.363 2.41e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3145 on 146 degrees of freedom
Multiple R-squared:  0.8586,    Adjusted R-squared:  0.8557
F-statistic: 295.5 on 3 and 146 DF,  p-value: < 2.2e-16

从这份输出的结果中，我们可以看出各个解释变数的系数，若将整个模型写出来就会像这样：

Sepal.Length = 0.65084 * Sepal.Width
             + 0.70913 * Petal.Length
             - 0.55648 * Petal.Width
             + 1.856

报表中的 Pr(>|t|) 就是统计上的 p-value，以这里的值来说，每一个系数都非常显著。
R-squared 的值为 0.8586，Adjusted R-squared 的值为 0.8557，表示模型的配适情况不错

模型诊断

在配适完回归模型之后，接着要进行模型诊断，通常会先画出几张常用来诊断模型的图，这里我们用 ggfortify 这个套件来画图。

1
2
3

library(ggfortify)
# 画出模型诊断用的图
autoplot(iris.lm)

在理论上回归模型配适完成后，其残差值（residual）必须符合常态性（normality）、独立性（independence）以及变异数同质性（homogeneity of variance）三种假设，所以接下来要用三种检定检查这三个条件是否符合。

残差独立性检定

独立性检定要使用到 car 这个套件：

1
2
3

# 残差独立性检定
require(car)
durbinWatsonTest(iris.lm)

1
2
3

lag Autocorrelation D-W Statistic p-value
  1     -0.03992126      2.060382   0.822
Alternative hypothesis: rho != 0

残差独立性检定的 p-value 也非常高，所以也不拒绝虚无假设，亦即残差值的独立性假设是合理的。

残差变异数同质性检定

残差的变异数同质性检定可以使用 car 套件所提供的 ncvTest（）函数：

1
2
3

# 残差变异数同质性检定
require(car)
ncvTest(iris.lm)

1
2
3

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 4.448612    Df = 1     p = 0.03492962

残差变异数同质性检定的 p-value 稍微偏小，以一般 95% 的信赖水准来说，是拒绝虚无假设的，也就是说残差的变异数没有符合同质性的假设，但是因为这个 p-value 并没有非常小，所以证据并不是非常明确。

预测

在回归模型建立好之后，就可以利用这个模型来预测新的资料，假设我们收到一些新的观测值（解释变数 X）：

1
2
3

# 新观测值
new.iris <- data.frame(Sepal.Width=3.1, Petal.Length=1.6, Petal.Width=0.3)
new.iris

若要使用回归模型预测新的反应变数，可以使用 predict 函数

1 2	# 预测资料 predict(iris.lm, new.iris)

4.841259
预测出来的 Sepal.Length 值就是 4.841259。