上一篇:[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是预测变量。
a和b被称为系数常数。
建立回归的步骤
回归的简单例子是当人的身高已知时预测人的体重。 为了做到这一点,我们需要有一个人的身高和体重之间的关系。
创建关系的步骤是 -
进行收集高度和相应重量的观测值的样本的实验。
使用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)
以下是所使用的参数的说明 -
公式是表示x和y之间的关系的符号。
数据是应用公式的向量。
创建关系模型并获取系数
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()
当我们执行上面的代码,它产生以下结果 -
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()
当我们执行上面的代码,它产生以下结果 -
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()
当我们执行上面的代码,它产生以下结果 -
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()
当我们执行上面的代码,它产生以下结果 -
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
时间序列图 -
不同的时间间隔
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
多时间序列图 -
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
我们可以得出结论,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
结论
从上面显示的决策树,我们可以得出结论,其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
上图中的趋势有助于我们预测在特定天数结束时的生存概率。
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,这表明字符串相关。