R语言学习笔记(15)——统计示例

平均值、中位数和模式

mean 平均值

通过求出数据集的和再除以求和数的总量得到平均值
函数mean()用于在R语言中计算平均值

mean()函数语法

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

x:输入向量
trim:用于设置计算均值前去掉两端数据的百分比,即计算结尾均值,取值在0~0.5之间,默认值为0
na.rm:指示是否允许有缺失值(NA)的情况,默认值为FALSE,即不允许
示例:
输入:

# 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()函数来计算此值。

median()函数语法

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模式

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

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"

附:
① unique()函数
用于从向量、 DataFrame 或数组中删除重复的元素/行。
② match()函数
匹配两个向量,返回第二个向量在第一个向量匹配位置的下标值。
③ tabulate()函数
用于计算向量中整数值的出现次数。
④ which.max()函数
查看向量中最大值对应的索引位置、如果向量中存在多个相等的最大、最小值返回的是第一个最大、最小值的索引

线性回归

回归分析是一种非常广泛使用的统计工具,用于建立两个变量之间的关系模型。 这些变量之一称为预测变量,其值通过实验收集。 另一个变量称为响应变量,其值从预测变量派生。
在线性回归中,这两个变量通过方程相关,其中这两个变量的指数(幂)为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(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(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()

输出:
在这里插入图片描述

多重回归

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

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

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

示例

输入数据

考虑在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

逻辑回归

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

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

y:响应变量
x1,x2,…:预测变量
a,b1,b2,…:数字常数的系数
用于创建回归模型的函数是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

结论

对于变量“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()函数

该函数给出给定平均值和标准偏差在每个点的概率分布的高度。即概率密度函数。d表示density
输入:

# 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()函数

该函数给出正态分布随机数的概率小于给定数的值。 它也被称为“累积分布函数”。即分布函数。p表示Probability
输入:

# 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()函数

该函数采用概率值,并给出累积值与概率值匹配的数字。qnorm()函数是pnorm()函数的反函数。q指的是quantile,即分位数
输入:

# 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

rnorm()函数

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

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

print(x)

输出:

[1] 57 60 64 53 68 66 56 70

泊松回归

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

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

y:响应变量。
a,b:数字系数。
x:预测变量。
用于创建泊松回归模型的函数是glm()函数
(该函数语法见逻辑回归处)

示例

我们有内置的数据集“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对断裂计数有影响。

协方差分析

我们使用回归分析创建模型,描述变量在预测变量对响应变量的影响。 有时,如果我们有一个类别变量,如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))

输出:

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

时间序列分析

时间序列是一系列数据点,其中每个数据点与时间相关联。 一个简单的例子是股票在某一天的不同时间点的股票价格。 另一个例子是一个地区在一年中不同月份的降雨量。 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

在这里插入图片描述

非线性最小二乘

当模拟真实世界数据用于回归分析时,我们观察到,很少情况下,模型的方程是给出线性图的线性方程。大多数时候,真实世界数据模型的方程涉及更高程度的数学函数,如3的指数或sin函数。在这种情况下,模型的图给出了曲线而不是直线。线性和非线性回归的目的是调整模型参数的值,以找到最接近您的数据的线或曲线。在找到这些值时,我们将能够以良好的精确度估计响应变量。
在最小二乘回归中,我们建立了一个回归模型,其中来自回归曲线的不同点的垂直距离的平方和被最小化。我们通常从定义的模型开始,并假设系数的一些值。然后我们应用R语言的nls()函数获得更准确的值以及置信区间。

nls()函数语法

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。

决策树

决策树是以树的形式表示选择及其结果的图。图中的节点表示事件或选择,并且图的边缘表示决策规则或条件。它主要用于使用R的机器学习和数据挖掘应用程序。
决策树的使用的例子是:预测电子邮件是垃圾邮件或非垃圾邮件,预测肿瘤癌变,或者基于这些因素预测贷款的信用风险。通常,使用观测数据(也称为训练数据)来创建模型。然后使用一组验证数据来验证和改进模型。 R具有用于创建和可视化决策树的包。对于新的预测变量集合,我们使用此模型来确定R包“party”用于创建决策树。

安装R包

install.packages("party")

“party”包具有用于创建和分析决策树的函数ctree()。
若安装报错,一般换个镜像即可。

ctree()函数语法

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包

install.packages("randomForest")

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

randomForest()函数语法

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
Loading required package: methods
Loading required package: grid
...............................
...............................

示例

我们将使用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%
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%。

生存分析

生存分析处理预测特定事件将要发生的时间。 它也被称为故障时间分析或分析死亡时间。 例如,预测患有癌症的人将存活的天数或预测机械系统将失败的时间。
命名为survival的R语言包用于进行生存分析。 此包包含函数Surv(),它将输入数据作为R语言公式,并在选择的变量中创建一个生存对象用于分析。 然后我们使用函数survfit()创建一个分析图。

安装R包

install.packages("survival")

Surv()和survfit()函数语法

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 

在这里插入图片描述
上图中的趋势有助于我们预测在特定天数结束时的生存概率。

卡方检验

卡方检验是一种确定两个分类变量之间是否存在显着相关性的统计方法。 这两个变量应该来自相同的人口,他们应该是类似 ——是/否,男/女,红/绿等。
例如,我们可以建立一个观察人们的冰淇淋购买模式的数据集,并尝试将一个人的性别与他们喜欢的冰淇淋的味道相关联。 如果发现相关性,我们可以通过了解访问的人的性别的数量来计划适当的味道库存。

chisq.test()函数语法

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 ...
NULL

上述结果表明数据集有很多因素变量,可以被认为是分类变量。 对于我们的模型,我们将考虑变量“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

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

	Pearson's Chi-squared test

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

结果显示p值小于0.05,这表明字符串相关。

  • 2
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

HanWLang

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值