R 语言统计分析

前言

R 中的统计分析通过使用许多内置函数来执行。这些函数大多数是 R 基础包的一部分。这些函数将 R 向量作为输入和参数，并给出结果。

平均值

通过求出数据集的和再除以求和数的总量得到平均值

函数 mean() 用于在 R 语言中计算平均值。

语法

用于计算 R 中的平均值的基本语法是：

mean(x, trim = 0, na.rm = FALSE, ...)

以下是所使用的参数的说明：

• x 是输入向量。
• trim 用于从排序向量的两端丢弃一些观察结果。
• na.rm 用于从输入向量中删除缺失值。

# Create a vector. 
x <- c(12, 7, 3, 4.2, 18, 2, 54, -21, 8, -5)

# Find Mean.
result.mean <- mean(x)
print(result.mean)

# 当我们执行上面的代码，它产生以下结果：

[1] 8.22

应用 trim 选项

当提供 trim 参数时，向量中的值被排序，然后从计算平均值中减去所需的观察值。

当 trim = 0.3 时，来自每端的 3 个值将从计算中减去以找到均值。

在这种情况下，排序的向量是 (-21, -5, 2, 3, 4.2, 7, 8, 12, 18, 54)，并且从用于计算平均值的向量中移除的值是左边的 (-21, -5, 2) 从和右边的 (12, 18, 54)。

# Create a vector.
x <- c(12, 7, 3, 4.2, 18, 2, 54, -21, 8, -5)

# Find Mean.
result.mean <- mean(x, trim = 0.3)
print(result.mean)

# 当我们执行上面的代码，它产生以下结果：

[1] 5.55

应用 NA 选项

如果有缺失值，则平均函数返回 NA。

要从计算中删除缺少的值，请使用 na.rm = TRUE。这意味着去除 NA 值。

# Create a vector. 
x <- c(12, 7, 3, 4.2, 18, 2, 54, -21, 8, -5, NA)

# Find mean.
result.mean <- mean(x)
print(result.mean)

# Find mean dropping NA values.
result.mean <- mean(x, na.rm = TRUE)
print(result.mean)

# 当我们执行上面的代码，它产生以下结果：

[1] NA
[1] 8.22

中位数

数据系列中的最中间值称为中位数。在 R 语言中使用 median() 函数来计算此值。

语法

计算 R 语言中位数的基本语法是：

median(x, na.rm = FALSE)

以下是所使用的参数的说明：

• x 是输入向量。
• na.rm 用于从输入向量中删除缺失值。

# Create the vector.
x <- c(12, 7, 3, 4.2, 18, 2, 54, -21, 8, -5)

# Find the median.
median.result <- median(x)
print(median.result)

# 当我们执行上面的代码，它产生以下结果：

[1] 5.6

众数

众数是一组数据中出现次数最多的值，不同于平均值和中位数，众数可以同时包含数字和字符数据。R 没有标准的内置函数来计算众数。因此，我们将创建一个用户自定义函数来计算 R 语言中的数据集的众数。该函数将向量作为输入，并将众数值作为输出。

# Create the function.
getmode <- function(v) {
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}

# Create the vector with numbers.
v <- c(2, 1, 2, 3, 1, 2, 3, 4, 1, 5, 5, 3, 2, 3)

# Calculate the mode using the user function.
result <- getmode(v)
print(result)

# Create the vector with characters.
charv <- c("o", "it", "the", "it", "it")

# Calculate the mode using the user function.
result <- getmode(charv)
print(result)

# 当我们执行上面的代码，它产生以下结果：

[1] 2
[1] "it"

线性回归

回归分析是一种非常广泛使用的统计工具，用于建立两个变量之间的关系模型。这些变量之一称为预测变量，其值通过实验收集。另一个变量称为响应变量，其值从预测变量派生。

在线性回归中，这两个变量通过方程相关，其中这两个变量的指数（幂）为 1，数学上，线性关系表示当绘制为曲线图时的直线。任何变量的指数不等于 1 的非线性关系将创建一条曲线。

线性回归的一般数学方程为：

y = ax + b

以下是所使用的参数的说明：

• y 是响应变量。
• x 是预测变量。
• a，b 被称为系数常数。

建立回归的步骤

回归的简单例子是当人的身高已知时预测人的体重。为了做到这一点，我们需要有一个人的身高和体重之间的关系。

创建关系的步骤是：

1. 进行收集高度和相应重量的观测值的样本的实验。
2. 使用 R 语言中的 lm() 函数创建关系模型。
3. 从创建的模型中找到系数，并使用这些创建数学方程
4. 获得关系模型的摘要以了解预测中的平均误差。也称为残差。
5. 为了预测新人的体重，使用 R 中的 predict() 函数。

输入数据

下面是代表观察的样本数据：

# Values of height
151, 174, 138, 186, 128, 136, 179, 163, 152, 131

# Values of weight.
63, 81, 56, 91, 47, 57, 76, 72, 62, 48

lm()函数

此函数创建预测变量和响应变量之间的关系模型。

语法

线性回归中 lm() 函数的基本语法是：

lm(formula, data)

以下是所使用的参数的说明：

• formula 是表示 x 和 y 之间的关系的符号。
• data 是应用公式的向量。

创建关系模型并获取系数

x <- c(151, 174, 138, 186, 128, 136, 179, 163, 152, 131)
y <- c(63, 81, 56, 91, 47, 57, 76, 72, 62, 48)

# Apply the lm() function.
relation <- lm(y ~ x)

print(relation)

# 当我们执行上面的代码，它产生以下结果：

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x
   -38.4551       0.6746

获取相关的摘要

x <- c(151, 174, 138, 186, 128, 136, 179, 163, 152, 131)
y <- c(63, 81, 56, 91, 47, 57, 76, 72, 62, 48)

# Apply the lm() function.
relation <- lm(y ~ x)

print(summary(relation))

# 当我们执行上面的代码，它产生以下结果：

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3002 -1.6629  0.0412  1.8944  3.9775 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -38.45509    8.04901  -4.778  0.00139 ** 
x             0.67461    0.05191  12.997 1.16e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.253 on 8 degrees of freedom
Multiple R-squared:  0.9548,	Adjusted R-squared:  0.9491 
F-statistic: 168.9 on 1 and 8 DF,  p-value: 1.164e-06

predict() 函数

语法

线性回归中的 predict() 的基本语法是：

predict(object, newdata)

以下是所使用的参数的说明：

• object 是已使用 lm() 函数创建的公式。
• newdata 是包含预测变量的新值的向量。

预测新人的体重

# The predictor vector.
x <- c(151, 174, 138, 186, 128, 136, 179, 163, 152, 131)

# The resposne vector.
y <- c(63, 81, 56, 91, 47, 57, 76, 72, 62, 48)

# Apply the lm() function.
relation <- lm(y ~ x)

# Find weight of a person with height 170.
a <- data.frame(x = 170)
result <- predict(relation, a)
print(result)

# 当我们执行上面的代码，它产生以下结果：

       1 
76.22869 

以图形方式可视化回归

# Create the predictor and response variable.
x <- c(151, 174, 138, 186, 128, 136, 179, 163, 152, 131)
y <- c(63, 81, 56, 91, 47, 57, 76, 72, 62, 48)
relation <- lm(y ~ x)

# Give the chart file a name.
png(file = "linear_regression.png")

# Plot the chart.
plot(y, x, col = "blue", main = "Height & Weight Regression", abline(lm(x ~ y)), cex = 1.3, pch = 16, xlab = "Weight in Kg", ylab = "Height in cm")

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下结果：

多元回归

多元回归是线性回归到两个以上变量之间的关系的延伸。在简单线性关系中，我们有一个预测变量和一个响应变量，但在多元回归中，我们有多个预测变量和一个响应变量。

多元回归的一般数学方程为：

y = a + b1x1 + b2x2 + ... + bnxn

以下是所使用的参数的说明：

• y 是响应变量。
• a, b1, b2, ..., bn 是系数。
• x1, x2, ..., xn 是预测变量。

我们使用 R 语言中的 lm() 函数创建回归模型。模型使用输入数据确定系数的值。接下来，我们可以使用这些系数来预测给定的一组预测变量的响应变量的值。

lm() 函数函数创建预测变量和响应变量之间的关系模型。

语法

lm() 函数在多元回归中的基本语法是：

lm(y ~ x1 + x2 + x3 + ..., data)

以下是所使用的参数的说明：

• y ~ x1 + x2 + x3... 是表示响应变量和预测变量之间的关系的符号。
• data 是应用公式的向量。

实例

输入数据

考虑在 R 语言环境中可用的数据集 mtcars 。它给出了每加仑里程 mpg，气缸排量 disp，马力 hp，汽车重量 wt 和一些其他参数的不同汽车模型之间的比较。

模型的目标是建立 mpg 作为响应变量与 disp，hp 和 wt 作为预测变量之间的关系。为此，我们从 mtcars 数据集中创建这些变量的子集。

input <- mtcars[, c("mpg", "disp", "hp", "wt")]
print(head(input))

# 当我们执行上面的代码，它产生以下结果：

                   mpg disp  hp    wt
Mazda RX4         21.0  160 110 2.620
Mazda RX4 Wag     21.0  160 110 2.875
Datsun 710        22.8  108  93 2.320
Hornet 4 Drive    21.4  258 110 3.215
Hornet Sportabout 18.7  360 175 3.440
Valiant           18.1  225 105 3.460

创建关系模型并获取系数

input <- mtcars[, c("mpg", "disp", "hp", "wt")]

# Create the relationship model.
model <- lm(mpg ~ disp + hp + wt, data = input)

# Show the model.
print(model)

# Get the Intercept and coefficients as vector elements.
cat("# # # # The Coefficient Values # # # #", "")

a <- coef(model)[1]
print(a)

Xdisp <- coef(model)[2]
Xhp <- coef(model)[3]
Xwt <- coef(model)[4]

print(Xdisp)
print(Xhp)
print(Xwt)

# 当我们执行上面的代码，它产生以下结果：

Call:
lm(formula = mpg ~ disp + hp + wt, data = input)

Coefficients:
(Intercept)         disp           hp           wt  
  37.105505    -0.000937    -0.031157    -3.800891 

# # # # The Coefficient Values # # # #
(Intercept) 
   37.10551 

         disp 
-0.0009370091
         hp 
-0.03115655 
       wt 
-3.800891 

创建回归模型的方程

基于上述截距和系数值，我们创建了数学方程。

Y = a + Xdisp*x1 + Xhp*x2 + Xwt*x3
or
Y = 37.105505 + (-0.000937)*x1 + (-0.031157)*x2 + (-3.800891)*x3

应用方程预测新值

当提供一组新的位移，马力和重量值时，我们可以使用上面创建的回归方程来预测里程数。 对于 disp = 221，hp = 102 和 wt = 2.91 的汽车，预测里程为：

Y = 37.105505 + (-0.000937)*221 + (-0.031157)*102 + (-3.800891)*2.91 = 22.65982

使用 predict() 函数预测

x = data.frame(disp = 221, hp = 102, wt = 2.91)
result = predict(model, x)
print(result)

# 当我们执行上面的代码，它产生以下结果：

       1 
22.65987

逻辑回归

逻辑回归是回归模型，其中响应变量（因变量）具有诸如 True/False 或 0/1 的分类值。它实际上基于将其与预测变量相关的数学方程测量二元响应的概率作为响应变量的值。

逻辑回归的一般数学方程为：

y = 1/(1 + e^-(a + b1x1 + b2x2 + b3x3 + ...))

以下是所使用的参数的说明：

• y 是响应变量。
• x 是预测变量。
• a 和 b 是作为数字常数的系数。

用于创建回归模型的函数是 glm() 函数。

语法

逻辑回归中 glm() 函数的基本语法是：

glm(formula, data, family)

以下是所使用的参数的说明：

• formula 是表示变量之间的关系的符号。
• data 是给出这些变量的值的数据集。
• family 是 R 语言对象来指定模型的细节。它的值是二项逻辑回归。

实例

内置数据集 mtcars 描述具有各种发动机规格的、不同型号的汽车。在 mtcars 数据集中，传输模式（自动或手动）由 am 列描述，它是一个二进制值 0 或 1。我们可以在列 am 和其他 3 列（hp，wt 和 cyl）之间创建逻辑回归模型。

# Select some columns form mtcars.
input <- mtcars[, c("am", "cyl", "hp", "wt")]

print(head(input))

# 当我们执行上面的代码，它产生以下结果：

                  am   cyl  hp    wt
Mazda RX4          1   6    110   2.620
Mazda RX4 Wag      1   6    110   2.875
Datsun 710         1   4     93   2.320
Hornet 4 Drive     0   6    110   3.215
Hornet Sportabout  0   8    175   3.440
Valiant            0   6    105   3.460

创建回归模型

我们使用 glm() 函数创建回归模型，并得到其摘要进行分析。

input <- mtcars[, c("am", "cyl", "hp", "wt")]

am.data = glm(formula = am ~ cyl + hp + wt, data = input, family = binomial)

print(summary(am.data))

# 当我们执行上面的代码，它产生以下结果：

Call:
glm(formula = am ~ cyl + hp + wt, family = binomial, data = input)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.17272  -0.14907  -0.01464   0.14116   1.27641  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) 19.70288    8.11637   2.428   0.0152 *
cyl          0.48760    1.07162   0.455   0.6491  
hp           0.03259    0.01886   1.728   0.0840 .
wt          -9.14947    4.15332  -2.203   0.0276 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.2297  on 31  degrees of freedom
Residual deviance:  9.8415  on 28  degrees of freedom
AIC: 17.841

Number of Fisher Scoring iterations: 8

结论

在 summary 中，对于变量 cyl 和 hp ，最后一列中的 P 值大于 0.05，我们认为它们对变量 am 的值有贡献是无关紧要的。只有重量 wt 影响该回归模型中的 am 值。

标准分布

在来自独立源的数据的随机集合中，通常观察到数据的分布是正常的。这意味着，在绘制水平轴上的变量值和垂直轴上的值的计数的图形时，我们得到钟形曲线。曲线的中心表示数据集的平均值。在图中，50% 的值位于平均值的左侧，另外 50% 位于图表的右侧。这在统计学中被称为正态分布。

R 语言有四个内置函数来产生正态分布。它们描述如下：

语法

dnorm(x, mean, sd)
pnorm(x, mean, sd)
qnorm(p, mean, sd)
rnorm(n, mean, sd)

以下是所使用的参数的说明：

• x 是数字的向量。
• P 是概率的向量。
• n 是观察的数量（样本大小）。
• mean 是样本数据的平均值。它的默认值为 0。
• sd 是标准偏差。它的默认值为 1。

dnorm()

该函数给出给定平均值和标准偏差在每个点的概率分布的高度。

# Create a sequence of numbers between -10 and 10 incrementing by 0.1.
x <- seq(-10, 10, by = .1)

# Choose the mean as 2.5 and standard deviation as 0.5.
y <- dnorm(x, mean = 2.5, sd = 0.5)

# Give the chart file a name.
png(file = "dnorm.png")

plot(x, y)

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下结果：

pnorm()

该函数给出正态分布随机数的概率小于给定数的值。它也被称为“累积分布函数”。

# Create a sequence of numbers between -10 and 10 incrementing by 0.2.
x <- seq(-10, 10, by = .2)
 
# Choose the mean as 2.5 and standard deviation as 2. 
y <- pnorm(x, mean = 2.5, sd = 2)

# Give the chart file a name.
png(file = "pnorm.png")

# Plot the graph.
plot(x, y)

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下结果：

qnorm()

该函数采用概率值，并给出累积值与概率值匹配的数字。

# Create a sequence of probability values incrementing by 0.02.
x <- seq(0, 1, by = 0.02)

# Choose the mean as 2 and standard deviation as 3.
y <- qnorm(x, mean = 2, sd = 1)

# Give the chart file a name.
png(file = "qnorm.png")

# Plot the graph.
plot(x, y)

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下结果：

rnorm()

此函数用于生成分布正常的随机数。它将样本大小作为输入，并生成许多随机数。我们绘制一个直方图来显示生成的数字的分布。

# Create a sample of 50 numbers which are normally distributed.
y <- rnorm(50)

# Give the chart file a name.
png(file = "rnorm.png")

# Plot the histogram for this sample.
hist(y, main = "Normal DIstribution")

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下结果：

二项分布

二项分布模型处理在一系列实验中仅发现两个可能结果的事件的成功概率。例如，掷硬币总是给出正或反。在二项分布期间估计在 10 次重复抛掷硬币中精确找到 3 个正的概率。

R 语言有四个内置函数来生成二项分布。它们描述如下：

dbinom(x, size, prob)
pbinom(x, size, prob)
qbinom(p, size, prob)
rbinom(n, size, prob)

以下是所使用的参数的说明：

• x 是数字的向量。
• P 是概率向量。
• n 是观察的数量。
• size 是试验的数量。
• prob 是每个试验成功的概率。

dbinom()

该函数给出每个点的概率密度分布。

# Create a sample of 50 numbers which are incremented by 1.
x <- seq(0, 50, by = 1)

# Create the binomial distribution.
y <- dbinom(x, 50, 0.5)

# Give the chart file a name.
png(file = "dbinom.png")

# Plot the graph for this sample.
plot(x, y)

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下结果：

pbinom()

此函数给出事件的累积概率。它是表示概率的单个值。

# Probability of getting 26 or less heads from a 51 tosses of a coin.
x <- pbinom(26, 51, 0.5)
print(x)

# 当我们执行上面的代码，它产生以下结果：

[1] 0.610116

qbinom()

该函数采用概率值，并给出累积值与概率值匹配的数字。

# How many heads will have a probability of 0.25 will come out when a coin is tossed 51 times.
x <- qbinom(0.25, 51, 1/2)
print(x)

# 当我们执行上面的代码，它产生以下结果：

[1] 23

rbinom()

该函数从给定样本产生给定概率的所需数量的随机值。

# Find 8 random values from a sample of 150 with probability of 0.4.
x <- rbinom(8, 150, .4)
print(x)

# 当我们执行上面的代码，它产生以下结果：

[1] 63 62 64 54 54 47 57 59

泊松回归

泊松回归（Poisson regression）是用来为计数资料和列联表建模的一种回归分析，其中响应变量是计数而不是分数的形式。

例如，在一个足球系列赛中出线或获胜的次数。此外，响应变量的值遵循泊松分布。

泊松回归的一般数学方程为：

log(y) = a + b1x1 + b2x2 + bnxn ...

以下是所使用的参数的说明：

• x​ 是预测变量。
• y ​是响应变量。
• ​a ​和 ​b ​是数字系数。

用于创建泊松回归模型的函数是​ glm()​ 函数。

语法

在泊松回归中 ​glm()​ 函数的基本语法是：

glm(formula, data, family)

以下是所使用的参数的说明：

• formula ​是表示变量之间的关系的符号。
• ​data​ 是给出这些变量的值的数据集。
• ​family​ 是 R 语言对象来指定模型的细节。它的值是“泊松”的逻辑回归。

实例

我们有内置的数据集 ​warpbreaks，其描述了羊毛类型 wool（​A​ 或 ​B​）和张力 tension（低，中或高）对每个织机的经纱断裂数量 breaks 的影响。让我们考虑“断裂”作为响应变量，它是断裂次数的计数。羊毛“类型”和“张力”作为预测变量。

输入数据

input <- warpbreaks
print(head(input))

# 当我们执行上面的代码，它产生以下结果：

  breaks wool tension
1     26    A       L
2     30    A       L
3     54    A       L
4     25    A       L
5     70    A       L
6     52    A       L

创建回归模型

output <-glm(formula = breaks ~ wool + tension, data = warpbreaks, family = poisson)
print(summary(output))

# 当我们执行上面的代码，它产生以下结果：

Call:
glm(formula = breaks ~ wool + tension, family = poisson, data = warpbreaks)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.6871  -1.6503  -0.4269   1.1902   4.2616  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 297.37  on 53  degrees of freedom
Residual deviance: 210.39  on 50  degrees of freedom
AIC: 493.06

Number of Fisher Scoring iterations: 4

在 summary 中，我们查找最后一列中的 ​p ​值小于​ 0.05​，以考虑预测变量对响应变量的影响。如图所示，具有张力类型 ​M ​和​ H​ 的羊毛类型 ​B ​对断裂计数有影响。

协方差分析

我们使用回归分析创建模型，描述变量在预测变量对响应变量的影响。有时，如果我们有一个类别变量，如 Yes/No 或 Male/Female 等。简单的回归分析为分类变量的每个值提供多个结果。在这种情况下，我们可以通过将分类变量与预测变量一起使用并比较分类变量的每个级别的回归线来研究分类变量的效果。这样的分析被称为协方差分析，也称为 ANCOVA。

实例

考虑在 R 语言中内置的数据集 mtcars 。在其中我们观察到字段 am 表示传输的类型（自动或手动）。它是值为 0/1 的分类变量。汽车的每加仑英里数 mpg 也可以取决于马力 hp 的值。

我们研究 am 的值对 mpg 和 hp 之间回归的影响。它是通过使用 aov() 函数，然后使用 anova() 函数来比较多个回归来完成的。

输入数据

从数据集 mtcars 创建一个包含字段 mpg，hp 和 am 的数据框。这里我们取 mpg 作为响应变量，hp 作为预测变量，am 作为分类变量。

input <- mtcars[, c("am", "mpg", "hp")]
print(head(input))

# 当我们执行上面的代码，它产生以下结果：

                  am  mpg  hp
Mazda RX4          1 21.0 110
Mazda RX4 Wag      1 21.0 110
Datsun 710         1 22.8  93
Hornet 4 Drive     0 21.4 110
Hornet Sportabout  0 18.7 175
Valiant            0 18.1 105

协方差分析

我们创建一个回归模型，以 hp 作为预测变量，mpg 作为响应变量，考虑 am 和 hp 之间的相互作用。

模型与分类变量和预测变量之间的相互作用

# Get the dataset.
input <- mtcars

# Create the regression model.
result <- aov(mpg ~ hp * am, data = input)
print(summary(result))

# 当我们执行上面的代码，它产生以下结果：

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4  77.391 1.50e-09 ***
am           1  202.2   202.2  23.072 4.75e-05 ***
hp:am        1    0.0     0.0   0.001    0.981    
Residuals   28  245.4     8.8                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这个结果表明，马力 hp 和传输类型 am 对每加仑的英里 mpg 有显着的影响，因为两种情况下的 P 值都小于 0.05。但是这两个变量之间的相互作用不显着，因为 hp:am 的 P 值大于 0.05。

没有分类变量和预测变量之间相互作用的模型

# Get the dataset.
input <- mtcars

# Create the regression model.
result <- aov(mpg ~ hp + am, data = input)
print(summary(result))

# 当我们执行上面的代码，它产生以下结果：

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4   80.15 7.63e-10 ***
am           1  202.2   202.2   23.89 3.46e-05 ***
Residuals   29  245.4     8.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这个结果表明，马力 hp 和传输类型 am 对每加仑的英里 mpg 有显着的影响，因为两种情况下的 P 值都小于 0.05。

比较两个模型

现在我们可以比较两个模型来得出结论，变量的相互作用是否真正重要。为此，我们使用 anova() 函数。

# Get the dataset.
input <- mtcars

# Create the regression models.
result1 <- aov(mpg ~ hp * am, data = input)
result2 <- aov(mpg ~ hp + am, data = input)

# Compare the two models.
print(anova(result1, result2))

# 当我们执行上面的代码，它产生以下结果：

Analysis of Variance Table

Model 1: mpg ~ hp * am
Model 2: mpg ~ hp + am
  Res.Df    RSS Df  Sum of Sq     F Pr(>F)
1     28 245.43                           
2     29 245.44 -1 -0.0052515 6e-04 0.9806

由于 P 值大于 0.05，我们得出结论，马力 hp 和传播类型 am 之间的相互作用不显着。因此，在汽车和手动变速器模式下，每加仑的里程将以类似的方式取决于汽车的马力。

时间序列分析

时间序列是将统一统计值按照时间发生的先后顺序来进行排列，时间序列分析的主要目的是根据已有数据对未来进行预测。

一个稳定的时间序列中常常包含两个部分，那么就是：有规律的时间序列 + 噪声。所以，在以下的方法中，主要的目的就是去过滤噪声值，让我们的时间序列更加的有分析意义。

语法

时间序列分析中 ts() 函数的基本语法是：

timeseries.object.name <- ts(data, start, end, frequency)

以下是所使用的参数的说明：

• data 是包含在时间序列中使用的值的向量或矩阵。
• start 以时间序列指定第一次观察的开始时间。
• end 指定时间序列中最后一次观测的结束时间。
• frequency 指定每单位时间的观测数。

除了参数 data，所有其他参数是可选的。

时间序列的预处理

1. 平稳性检验

   拿到一个时间序列之后，我们首先要对其稳定性进行判断，只有非白噪声的稳定性时间序列才有分析的意义以及预测未来数据的价值。

   所谓平稳，是指统计值在一个常数上下波动并且波动范围是有界限的。如果有明显的趋势或者周期性，那么就是不稳定的。一般判断有三种方法：

   • 画出时间序列的趋势图，看趋势判断
   • 画自相关图和偏相关图，平稳时间序列的自相关图和偏相关图，要么拖尾，要么截尾。

   • 检验序列中是否存在单位根，如果存在单位根，就是非平稳时间序列。

     在 R 语言中，DF 检测是一种检测稳定性的方法，如果得出的 P 值小于临界值，则认为是数列是稳定的。

2. 白噪声检验

   白噪声序列，又称为纯随机性序列，序列的各个值之间没有任何的相关关系，序列在进行无序的随机波动，可以终止对该序列的分析，因为从白噪声序列中是提取不到任何有价值的信息的。

3. 平稳时间序列的参数特点

   均值和方差为常数，并且具有与时间无关的自协方差。

时间序列建模步骤

1. 拿到被分析的时间序列数据集。

2. 对数据绘图，观测其平稳性。若为非平稳时间序列要先进行 d 阶差分运算后化为平稳时间序列，此处的 d 即为 ARIMA(p,d,q) 模型中的 d；若为平稳序列，则用 ARMA(p,q) 模型。所以 ARIMA(p,d,q) 模型区别于 ARMA(p,q) 之处就在于前者的自回归部分的特征多项式含有 d 个单位根。

3. 对得到的平稳时间序列分别求得其自相关系数 ACF 和偏自相关系数 PACF，通过对自相关图和偏自相关图的分析，得到最佳的阶层 p 和阶数 q。由以上得到的 d、q、p，得到 ARIMA 模型。

4. 模型诊断。进行诊断分析，以证实所得模型确实与所观察到的数据特征相符。若不相符，重新回到第 3 步。

实例

考虑从 2012 年 1 月开始的一个地方的年降雨量细节。我们创建一个 R 时间序列对象为期 12 个月并绘制它。

# Get the data points in form of a R vector.
rainfall <- c(799, 1174.8, 865.1, 1334.6, 635.4, 918.5, 685.5, 998.6, 784.2, 985, 882.8, 1071)

# Convert it to a time series object.
rainfall.timeseries <- ts(rainfall, start = c(2012, 1), frequency = 12)

# Print the timeseries data.
print(rainfall.timeseries)

# Give the chart file a name.
png(file = "rainfall.png")

# Plot a graph of the time series.
plot(rainfall.timeseries)

# Save the file.
dev.off()

# 当我们执行上面的代码，它产生以下结果及图表：

        Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct    Nov    Dec
2012  799.0 1174.8  865.1 1334.6  635.4  918.5  685.5  998.6  784.2  985.0  882.8 1071.0

不同的时间间隔

ts() 函数中的频率参数值决定了测量数据点的时间间隔。值为 12 表示时间序列为 12 个月。其他值及其含义如下：

• frequency = 12：指定一年中每个月的数据点。
• frequency = 4：每年的每个季度的数据点。
• frequency = 6：每小时的 10 分钟的数据点。
• frequency = 24 * 6：将一天的每 10 分钟的数据点固定。

多时间序列

我们可以通过将两个系列组合成一个矩阵，在一个图表中绘制多个时间序列。

# Get the data points in form of a R vector.
rainfall1 <- c(799, 1174.8, 865.1, 1334.6, 635.4, 918.5, 685.5, 998.6, 784.2, 985, 882.8, 1071)
rainfall2 <- c(655, 1306.9, 1323.4, 1172.2, 562.2, 824, 822.4, 1265.5, 799.6, 1105.6, 1106.7, 1337.8)

# Convert them to a matrix.
combined.rainfall <- matrix(c(rainfall1, rainfall2), nrow = 12)

# Convert it to a time series object.
rainfall.timeseries <- ts(combined.rainfall, start = c(2012, 1), frequency = 12)

# Print the timeseries data.
print(rainfall.timeseries)

# Give the chart file a name.
png(file = "rainfall_combined.png")

# Plot a graph of the time series.
plot(rainfall.timeseries, main = "Multiple Time Series")

# Save the file.
dev.off()

# 当我们执行上面的代码，它产生以下结果及图表：

         Series 1 Series 2
Jan 2012    799.0    655.0
Feb 2012   1174.8   1306.9
Mar 2012    865.1   1323.4
Apr 2012   1334.6   1172.2
May 2012    635.4    562.2
Jun 2012    918.5    824.0
Jul 2012    685.5    822.4
Aug 2012    998.6   1265.5
Sep 2012    784.2    799.6
Oct 2012    985.0   1105.6
Nov 2012    882.8   1106.7
Dec 2012   1071.0   1337.8

非线性最小二乘

当模拟真实世界数据用于回归分析时，我们观察到，很少情况下，模型的方程是给出线性图的线性方程。大多数时候，真实世界数据模型的方程涉及更高程度的数学函数，如 3 的指数或 sin 函数。在这种情况下，模型的图给出了曲线而不是线。线性和非线性回归的目的是调整模型参数的值，以找到最接近您的数据的线或曲线。在找到这些值时，我们将能够以良好的精确度估计响应变量。

在最小二乘回归中，我们建立了一个回归模型，其中来自回归曲线的不同点的垂直距离的平方和被最小化。我们通常从定义的模型开始，并假设系数的一些值。然后我们应用 R 语言的 nls() 函数获得更准确的值以及置信区间。

语法

在 R 语言中创建非线性最小二乘测试的基本语法是：

nls(formula, data, start)

以下是所使用的参数的说明：

• formula 是包括变量和参数的非线性模型公式。
• data 是用于计算公式中变量的数据框。
• start 是起始估计的命名列表或命名数字向量。

实例

我们将考虑一个假设其系数的初始值的非线性模型。接下来，我们将看到这些假设值的置信区间是什么，以便我们可以判断这些值在模型中有多好。

所以让我们考虑以下的方程 a = b1*x^2 + b2。

让我们假设初始系数为 b1 = 1 和 b2 = 3，并将这些值拟合到 nls() 函数中。

x <- c(1.6, 2.1, 2, 2.23, 3.71, 3.25, 3.4, 3.86, 1.19, 2.21)
y <- c(5.19, 7.43, 6.94, 8.11, 18.75, 14.88, 16.06, 19.12, 3.21, 7.58)

# Give the chart file a name.
png(file = "nls.png")

# Plot these values.
plot(x, y)


# Take the assumed values and fit into the model.
model <- nls(y ~ b1*x^2 + b2, start = list(b1 = 1, b2 = 3))

# Plot the chart with new data by fitting it to a prediction from 100 data points.
new.data <- data.frame(x = seq(min(x), max(x), len = 100))
lines(new.data$x, predict(model, newdata = new.data))

# Save the file.
dev.off()

# Get the sum of the squared residuals.
print(sum(resid(model)^2))

# Get the confidence intervals on the chosen values of the coefficients.
print(confint(model))

# 当我们执行上面的代码，它产生以下结果及图表：

[1] 1.081935

Waiting for profiling to be done...
       2.5%    97.5%
b1 1.137708 1.253135
b2 1.497364 2.496484

我们可以得出结论，b1 的值更接近 1，而 b2 的值更接近 2 而不是 3。

决策树

决策树是以树的形式表示选择及其结果的图。图中的节点表示事件或选择，并且图的边缘表示决策规则或条件。它主要用于使用 R 的机器学习和数据挖掘应用程序。

决策树的使用的例子是：预测电子邮件是垃圾邮件或非垃圾邮件，预测肿瘤癌变，或者基于这些因素预测贷款的信用风险。通常，使用观测数据（也称为训练数据）来创建模型。然后使用一组验证数据来验证和改进模型。R 具有用于创建和可视化决策树的包。对于新的预测变量集合，我们使用此模型来确定 R 包 party 用于创建决策树。

安装 R 语言包

在 R 语言控制台中使用以下命令安装软件包。您还必须安装相关软件包（如果有）。

install.packages("party")

party 包具有用于创建和分析决策树的函数 ctree()。

语法

在 R 中创建决策树的基本语法是：

ctree(formula, data)

以下是所使用的参数的说明：

• formula 是描述预测变量和响应变量的公式。
• data 是所使用的数据集的名称。

输入数据

我们将使用名为 readingSkills 的 R 内置数据集来创建决策树。它描述了某人的 readingSkills 的分数，如果我们知道变量 年龄， shoesize，分数，以及该人是否为母语者。

# Load the party package. It will automatically load other dependent packages.
library(party)

# Print some records from data set readingSkills.
print(head(readingSkills))

# 当我们执行上面的代码，它产生以下结果：

  nativeSpeaker age shoeSize    score
1           yes   5 24.83189 32.29385
2           yes   6 25.95238 36.63105
3            no  11 30.42170 49.60593
4           yes   7 28.66450 40.28456
5           yes  11 31.88207 55.46085
6           yes  10 30.07843 52.83124

实例

我们将使用 ctree() 函数创建决策树并查看其图形。

# Load the party package. It will automatically load other dependent packages.
library(party)

# Create the input data frame.
input.dat <- readingSkills[c(1:105), ]

# Give the chart file a name.
png(file = "decision_tree.png")

# Create the tree.
output.tree <- ctree(nativeSpeaker ~ age + shoeSize + score, data = input.dat)

# Plot the tree.
plot(output.tree)

# Save the file.
dev.off()

当我们执行上面的代码，它产生以下图表：

结论

从上面显示的决策树，我们可以得出结论，其 readingSkills 分数低于 38.3 和年龄超过 6 的人不是一个母语者。

随机森林

在随机森林方法中，创建大量的决策树。每个观察被馈入每个决策树。每个观察的最常见的结果被用作最终输出。新的观察结果被馈入所有的树并且对每个分类模型取多数投票。

对构建树时未使用的情况进行错误估计。这称为 OOB（袋外） 误差估计，其被提及为百分比。

R 语言包 randomForest 用于创建随机森林。

安装 R 包

在 R 语言控制台中使用以下命令安装软件包。您还必须安装相关软件包（如果有）。

install.packages("randomForest")

包 randomForest 具有函数 randomForest()，用于创建和分析随机森林。

语法

在 R 语言中创建随机森林的基本语法是：

randomForest(formula, data)

以下是所使用的参数的说明：

• formula 是描述预测变量和响应变量的公式。
• data 是所使用的数据集的名称。

输入数据

我们将使用名为 readingSkills 的 R 语言内置数据集来创建决策树。它描述了某人的 readingSkills 的分数，如果我们知道变量 age，shoesize，score，以及该人是否是母语者。

# Load the party package. It will automatically load other required packages.
library(party)

# Print some records from data set readingSkills.
print(head(readingSkills))

# 当我们执行上面的代码，它产生以下结果：

  nativeSpeaker age shoeSize    score
1           yes   5 24.83189 32.29385
2           yes   6 25.95238 36.63105
3            no  11 30.42170 49.60593
4           yes   7 28.66450 40.28456
5           yes  11 31.88207 55.46085
6           yes  10 30.07843 52.83124

实例

我们将使用 randomForest() 函数来创建决策树并查看它的图。

# Load the party package. It will automatically load other required packages.
library(party)
library(randomForest)

# Create the forest.
output.forest <- randomForest(nativeSpeaker ~ age + shoeSize + score, data = readingSkills)

# View the forest results.
print(output.forest)

# Importance of each predictor.
print(importance(fit, type = 2))

# 当我们执行上面的代码，它产生以下结果：

Call:
 randomForest(formula = nativeSpeaker ~ age + shoeSize + score, data = readingSkills) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 1

        OOB estimate of  error rate: 1.5%
Confusion matrix:
    no yes class.error
no  99   1        0.01
yes  2  98        0.02

         MeanDecreaseGini
age              13.95406
shoeSize         18.91006
score            56.73051

结论 从上面显示的随机森林，我们可以得出结论，鞋码 shoeSize 和成绩 score 是决定如果某人是母语者或不是母语的重要因素。此外，该模型只有 1% 的误差，这意味着我们可以预测精度为 99%。

生存分析

生存分析处理预测特定事件将要发生的时间。它也被称为故障时间分析或分析死亡时间。例如，预测患有癌症的人将存活的天数或预测机械系统将失败的时间。

命名为 survival 的 R 语言包用于进行生存分析。此包包含函数 Surv()，它将输入数据作为 R 语言公式，并在选择的变量中创建一个生存对象用于分析。然后我们使用函数 survfit() 创建一个分析图。

安装软件包

install.packages("survival")

语法

在 R 语言中创建生存分析的基本语法是：

Surv(time, event)
survfit(formula)

以下是所使用的参数的说明：

• time 是直到事件发生的跟踪时间。
• event 指示预期事件的发生的状态。
• formula 是预测变量之间的关系。

实例

我们将考虑在上面安装的 survival 包中存在的名为 pbc 的数据集。它描述了关于受肝原发性胆汁性肝硬化 PBC 影响的人的生存数据点。在数据集中存在的许多列中，我们主要关注字段 time 和 status。时间表示在接受肝移植或患者死亡的患者的登记和事件的较早之间的天数。

# Load the library.
library("survival")

# Print first few rows.
print(head(pbc))

# 当我们执行上面的代码，它产生以下结果：

  id time status trt      age sex ascites hepato spiders edema bili chol albumin copper alk.phos    ast trig platelet protime stage
1  1  400      2   1 58.76523   f       1      1       1   1.0 14.5  261    2.60    156   1718.0 137.95  172      190    12.2     4
2  2 4500      0   1 56.44627   f       0      1       1   0.0  1.1  302    4.14     54   7394.8 113.52   88      221    10.6     3
3  3 1012      2   1 70.07255   m       0      0       0   0.5  1.4  176    3.48    210    516.0  96.10   55      151    12.0     4
4  4 1925      2   1 54.74059   f       0      1       1   0.5  1.8  244    2.54     64   6121.8  60.63   92      183    10.3     4
5  5 1504      1   2 38.10541   f       0      1       1   0.0  3.4  279    3.53    143    671.0 113.15   72      136    10.9     3
6  6 2503      2   2 66.25873   f       0      1       0   0.0  0.8  248    3.98     50    944.0  93.00   63       NA    11.0     3

从上述数据，我们正在考虑分析的时间和状态。

应用 Surv() 和 survfit() 函数

现在我们继续应用 Surv() 函数到上面的数据集，并创建一个将显示趋势图。

# Load the library.
library("survival")

# Create the survival object. 
survfit(Surv(pbc$time, pbc$status == 2) ~ 1)

# Give the chart file a name.
png(file = "survival.png")

# Plot the graph. 
plot(survfit(Surv(pbc$time, pbc$status == 2) ~ 1))

# Save the file.
dev.off()

# 当我们执行上面的代码，它产生以下结果及图表：

Call: survfit(formula = Surv(pbc$time, pbc$status == 2) ~ 1)

      n  events  median 0.95LCL 0.95UCL 
    418     161    3395    3090    3853 

上图中的趋势有助于我们预测在特定天数结束时的生存概率。

卡方检验

卡方检验是一种确定两个分类变量之间是否存在显着相关性的统计方法。这两个变量应该来自相同的人口，他们应该是类似：是/否，男/女，红/绿 等。

例如，我们可以建立一个观察人们的冰淇淋购买模式的数据集，并尝试将一个人的性别与他们喜欢的冰淇淋的味道相关联。如果发现相关性，我们可以通过了解访问的人的性别的数量来计划适当的味道库存。

语法

用于执行卡方检验的函数是 chisq.test()。

在 R 语言中创建卡方检验的基本语法是：

chisq.test(data)

以下是所使用的参数的说明：

• data 是以包含观察中变量的计数值的表的形式的数据。

实例

我们将在 MASS 包中获取 Cars93 数据，代表 1993 年不同型号汽车的销售额。

library("MASS")
print(str(Cars93))

# 当我们执行上面的代码，它产生以下结果：

'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...

上述结果表明数据集有很多因素变量，可以被认为是分类变量。对于我们的模型，我们将考虑变量 AirBags 和 Type。在这里，我们的目标是找出所售的汽车类型和安全气囊类型之间的任何显着的相关性。如果观察到相关性，我们可以估计哪种类型的汽车可以更好地卖什么类型的气囊。

# Load the library.
library("MASS")

# Create a data frame from the main data set.
car.data <- data.frame(Cars93$AirBags, Cars93$Type)

# Create a table with the needed variables.
car.data = table(Cars93$AirBags, Cars93$Type) 
print(car.data)

# Perform the Chi-Square test.
print(chisq.test(car.data))

# 当我们执行上面的代码，它产生以下结果：

                     Compact Large Midsize Small Sporty Van
  Driver & Passenger       2     4       7     0      3   0
  Driver only              9     7      11     5      8   3
  None                     5     0       4    16      3   6

	Pearson`s Chi-squared test

data:  car.data
X-squared = 33.001, df = 10, p-value = 0.0002723

Warning message:
In chisq.test(car.data) : Chi-squared approximation may be incorrect

结论 结果显示 P 值小于 0.05，这表明汽车类型和安全气囊类型之间具有相关性。