R语言:统计应用篇

->点击访问个人博客地址,相互交流学习<-

上一篇:[R语言——进阶篇](https://www.jianshu.com/p/c0908bf810fe

1.平均值,中位数和模式

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

我们在本章中讨论的功能是平均值,中位数和模式。

Mean平均值

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

函数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 = 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

Median中位数

数据系列中的最中间值称为中值。 在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

Mode模式

模式是一组数据中出现次数最多的值。 Unike平均值和中位数,模式可以同时包含数字和字符数据。

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"

2.线性回归

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

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

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

y = ax + b

以下是所使用的参数的描述 -

  • y是响应变量。

  • x是预测变量。

  • ab被称为系数常数。

建立回归的步骤

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

创建关系的步骤是 -

  • 进行收集高度和相应重量的观测值的样本的实验。

  • 使用R语言中的lm()函数创建关系模型。

  • 从创建的模型中找到系数,并使用这些创建数学方程

  • 获得关系模型的摘要以了解预测中的平均误差。 也称为残差。

  • 为了预测新人的体重,使用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)

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

  • 公式是表示xy之间的关系的符号。

  • 数据是应用公式的向量。

创建关系模型并获取系数

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 = "linearregression.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()

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

R中线性回归

3.多重回归

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

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

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

以下是所使用的参数的描述 -

  • y是响应变量。

  • a,b1,b2 ... bn是系数。

  • x1,x2,... xn是预测变量。

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

lm()函数

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

语法

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

lm(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.15+(-0.000937)*x1+(-0.0311)*x2+(-3.8008)*x3

应用方程预测新值

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

Y = 37.15+(-0.000937)*221+(-0.0311)*102+(-3.8008)*2.91 = 22.7104

4.逻辑回归

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

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

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

以下是所使用的参数的描述 -

  • Ÿ是响应变量。

  • X是预测变量。

  • 一和b是作为数字常数的系数。

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

语法

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

glm(formula,data,family)

以下是所使用的参数的描述 -

  • 式是表示变量之间的关系的符号。

  • 数据是给出这些变量的值的数据集。

  • family是R语言对象来指定模型的细节。它的值是二项逻辑回归。

内置数据集“mtcars”描述具有各种发动机规格的汽车的不同型号。在“mtcars”数据集中,传输模式(自动或手动)由am列描述,它是一个二进制值(0或1)。我们可以在列“是”和其他3列(马力,重量和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

结论

在总结中,对于变量“cyl”和“hp”,最后一列中的p值大于0.05,我们认为它们对变量“am”的值有贡献是无关紧要的。只有重量(wt)影响该回归模型中的“是”值。

5.标准分布

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

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

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

以下是在上述功能中使用的参数的描述 -

  • X是数字的向量。

  • p是概率的向量。

  • Ñ是观察的数量(样本大小)。

  • mean是样本数据的平均值。它的默认值为零。

  • 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()

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

dnorm()图

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()

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

pnorm()图

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()

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

qnorm()图

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()

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

RNORM()图

6.二项分布

二项分布模型处理在一系列实验中仅发现两个可能结果的事件的成功概率。 例如,掷硬币总是给出头或尾。 在二项分布期间估计在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()

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

dbinom()图

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] 58 61 59 66 55 60 61 67

7.泊松回归

泊松回归包括回归模型,其中响应变量是计数而不是分数的形式。 例如,足球比赛系列中的出生次数或胜利次数。 此外,响应变量的值遵循泊松分布。

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

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

以下是所使用的参数的描述 -

  • y是响应变量。

  • a和b是数字系数。

  • x是预测变量。

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

语法

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

glm(formula,data,family)

以下是在上述功能中使用的参数的描述 -

  • formula是表示变量之间的关系的符号。

  • data是给出这些变量的值的数据集。

  • family是R语言对象来指定模型的细节。 它的值是“泊松”的逻辑回归。

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

输入数据

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

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

8.协方差分析

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

考虑在数据集mtcars中内置的R语言。 在其中我们观察到字段“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

这个结果表明,马力和传输类型对每加仑的英里有显着的影响,因为两种情况下的p值都小于0.05。 但是这两个变量之间的相互作用不显着,因为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

这个结果表明,马力和传输类型对每加仑的英里有显着的影响,因为两种情况下的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))

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

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,我们得出结论,马力和传播类型之间的相互作用不显着。 因此,在汽车和手动变速器模式下,每加仑的里程将以类似的方式取决于汽车的马力。

9.时间序列分析

时间序列是一系列数据点,其中每个数据点与时间戳相关联。 一个简单的例子是股票在某一天的不同时间点的股票价格。 另一个例子是一个地区在一年中不同月份的降雨量。 R语言使用许多函数来创建,操作和绘制时间序列数据。 时间序列的数据存储在称为时间序列对象的R对象中。 它也是一个R语言数据对象,如矢量或数据帧。

使用ts()函数创建时间序列对象。

语法

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

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

以下是所使用的参数的描述 -

  • data是包含在时间序列中使用的值的向量或矩阵。

  • start以时间序列指定第一次观察的开始时间。

  • end指定时间序列中最后一次观测的结束时间。

  • frequency指定每单位时间的观测数。

除了参数“data”,所有其他参数是可选的。

考虑从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
2012  799.0  1174.8  865.1  1334.6  635.4  918.5  685.5  998.6  784.2
        Oct    Nov    Dec
2012  985.0  882.8 1071.0

时间序列图 -

时间序列,使用R

不同的时间间隔

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

  • 频率= 12指定一年中每个月的数据点。

  • 频率= 4每年的每个季度的数据点。

  • 频率= 6每小时的10分钟的数据点。

  • 频率= 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

多时间序列图 -

结合时间序列,使用R

10.非线性最小二乘

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

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

语法

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

nls(formula, data, start)

以下是所使用的参数的描述 -

  • formula是包括变量和参数的非线性模型公式。

  • data是用于计算公式中变量的数据框。

  • start是起始估计的命名列表或命名数字向量。

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

所以让我们考虑下面的方程为这个目的 -

a = b1*x^2+b2

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

xvalues <- c(1.6,2.1,2,2.23,3.71,3.25,3.4,3.86,1.19,2.21)
yvalues <- 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(xvalues,yvalues)

# Take the assumed values and fit into the model.
model <- nls(yvalues ~ b1*xvalues^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(xvalues = seq(min(xvalues),max(xvalues),len = 100))
lines(new.data$xvalues,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
非线性至少方块R

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

11.决策树

决策树是以树的形式表示选择及其结果的图。图中的节点表示事件或选择,并且图的边缘表示决策规则或条件。它主要用于使用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
Loading required package: methods
Loading required package: grid
...............................
...............................

我们将使用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()

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

null device 
          1 
Loading required package: methods
Loading required package: grid
Loading required package: mvtnorm
Loading required package: modeltools
Loading required package: stats4
Loading required package: strucchange
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading required package: sandwich
决策树,使用R

结论

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

12.随机森林算法

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

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

[R语言包“随机森林”用于创建随机森林。

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

install.packages("randomForest)

包“随机森林”具有函数随机森林(),用于创建和分析随机森林。

语法

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

randomForest(formula, data)

以下是所使用的参数的描述 -

  • 式是描述预测变量和响应变量的公式。

  • 数据是所使用的数据集的名称。

输入数据

我们将使用名为readingSkills的R语言内置数据集来创建决策树。它描述了某人的阅读技能的分数,如果我们知道变量“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
Loading required package: methods
Loading required package: grid
...............................
...............................

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

# 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%
Confusion matrix:
    no yes class.error
no  99   1        0.01
yes  1  99        0.01
         MeanDecreaseGini
age              13.95406
shoeSize         18.91006
score            56.73051

结论

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

13.生存分析

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

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

安装软件包

install.packages("survival")

语法

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

Surv(time,event)
survfit(formula)

以下是所使用的参数的描述 -

  • time是直到事件发生的跟踪时间。

  • event指示预期事件的发生的状态。

  • formula是预测变量之间的关系。

我们将考虑在上面安装的生存包中存在的名为“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
1  1  400      2   1 58.76523   f       1      1       1   1.0 14.5  261
2  2 4500      0   1 56.44627   f       0      1       1   0.0  1.1  302
3  3 1012      2   1 70.07255   m       0      0       0   0.5  1.4  176
4  4 1925      2   1 54.74059   f       0      1       1   0.5  1.8  244
5  5 1504      1   2 38.10541   f       0      1       1   0.0  3.4  279
6  6 2503      2   2 66.25873   f       0      1       0   0.0  0.8  248
  albumin copper alk.phos    ast trig platelet protime stage
1    2.60    156   1718.0 137.95  172      190    12.2     4
2    4.14     54   7394.8 113.52   88      221    10.6     3
3    3.48    210    516.0  96.10   55      151    12.0     4
4    2.54     64   6121.8  60.63   92      183    10.3     4
5    3.53    143    671.0 113.15   72      136    10.9     3
6    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 
,使用R生存分析

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

14.卡方检验

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

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

语法

用于执行卡方检验的函数是chisq.test()。
在R语言中创建卡方检验的基本语法是 -

chisq.test(data)

以下是所使用的参数的描述 -

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

我们将在“大众”图书馆中获取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,这表明字符串相关。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 219,366评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,521评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,689评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,925评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,942评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,727评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,447评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,349评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,820评论 1 317
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,990评论 3 337
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,127评论 1 351
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,812评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,471评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,017评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,142评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,388评论 3 373
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,066评论 2 355

推荐阅读更多精彩内容