多元统计分析——基于R语言的单车使用情况可视化分析

注：基于R语言的单车使用情况可视化分析为实验记录，存在不足，自行改进。

一、提出问题（要解决或分析的问题）

1 、用户对共享单车的使用习惯，环境对共享单车运营带来的影响？

2 、共享单车的租赁量主要与哪些环境有关？

二、数据来源及选取方法（数据完整可靠）

数据说明：所使用的数据是共享单车运营数据，记录了共享单车租赁的时间、地点、环境（包括季节，温度，湿度）等数据。

数据来源： https://www.kaggle.com/c/bike-sharing-demand/data

三、数据可视化分析

3.1 数据基本信息描述（图或表）

datetime ( 时间 ): hourly date + timestamp

Season( 季节 ) : 1 = spring, 2 = summer, 3 = fall, 4 = winter

holiday ( 是否是周末 ): whether the day is considered a holiday

workingday ( 是否是工作日 ):whether the day is neither a weekend nor holiday

weather ( 天气 ): 1: Clear, Few clouds, Partly cloudy, Partly cloudy

2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist

3: Light Snow, Light Rain + Thunderstorm +Scattered clouds,Light Rain+Scattered clouds

4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog

temp ( 温度 ) : temperature in Celsius

atemp ( 最高温度 ) :"feels like" temperature in Celsius

humidity ( 湿度 ): relative humidity

windspeed ( 风速 ):wind speed

casual ( 损坏车辆 ) : number of non-registered user rentals initiated

registered ( 登记使用车辆 ) :number of registered user rentals initiated

count ( 总租赁车辆数 ):number of total rentals

3.2.采用的多元分析方法（一种或多种，必须理论阐述该研究方法的基本原理）

3.2.1采用主成分分析

主成分分析（principal components analysis ）也称主分量分析，是由霍特林于 1933 年首次提出的。主成分分析是利用降维的思想，在损失很少信息的前提下，把多个指标转化为几个综合指标的多元统计方法。

简单地说，主成分分析就是在研究问题选取的指标比较多、研究的问题比较复杂时，可以把原来研究的指标做几个线性组合，即主成分，来解释原来变量绝大对数信息的一种多元统计法。

3.2.2主成分分析的意义

通过主成分分析，可以从事物间错综复杂的关系中找出一些主要成分，从而能有效利用大量统计数据进行定量分析，揭示变量之间的内在关系，得到对事物特征发展规律的一些深层次的启发，把研究工作引向深入。

3.3数据分析结果及解释（图或表，必须含分析解释）

（1）读取数据

data <- read.csv("C:\\Users\\leglon\\Desktop\\假期r\\train(1).csv", header = TRUE,stringsAsFactors = F)

（2）缺失值处理

#缺失值处理
na.omit(data,cols,invert)

（3）从data中读取温度、最高温度、湿度、风速、损坏使用数据、登记使用数据，租赁总数据

# 从data中读取温度、最高温度、湿度、风速、损坏使用数据、登记使用数据，租赁总数据
suppressWarnings(temperature <- as.numeric(data[,6]))
suppressWarnings(atemp <- as.numeric(data[,7]))
suppressWarnings(humidity <- as.numeric(data[,8]))
suppressWarnings(windspeed <- as.numeric(data[,9]))
suppressWarnings(casual <- as.numeric(data[,10]))
suppressWarnings(register <- as.numeric(data[,11]))
suppressWarnings(count <- as.numeric(data[,12]))

（4）求温度、湿度和风速的均值与标准差

# 温度的均值和标准差
mean(temperature, na.rm = T)
sd(temperature, na.rm = T)

# 湿度的均值和标准差
mean(humidity, na.rm = T)
sd(humidity, na.rm = T)

# 风速的均值和标准差
mean(windspeed, na.rm = T)
sd(windspeed, na.rm = T)

（5）可视化湿度与租赁量关系图

plot(x = humidity,y=count,xlab="湿度",ylab="租赁量",main="湿度与租赁量关系图",cex=1,pch=1)

由湿度-租赁量散点图可知，湿度在 30-60 区间内，租赁量较多，湿度低于 20 ，且大于 90 ，

租赁量较少。

（6）可视化温度与租赁量关系图

plot(x=temperature,y=count,xlab="温度",ylab="租赁量",main="温度与租赁量关系图",cex=1,pch=1)

由温度-租赁量散点图可知，当温度在 15-25 之间时，共享单车租赁量较高；当温度在 15 度以下时，温度越低，租赁量越低；当温度较高时，仍能保持一定的租赁量。

（7）查看拟合效果

shapiro.test(temperature[0:5000])
#W的值越接近1就越表明数据和正态分布拟合得越好,可以认为样本数据服从正态分布.

w值越接近1就越表明数据和正态分布拟合得越好,可以认为样本数据服从正态分布，从图中可以看出w的值为0.97539，服从正态分布。

（8）可视化工作日与休息日单车使用情况

hours = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,22, 23, 24)

count_workday = c(36.73225806451613, 16.003236245954692, 8.436065573770492, 
                  4.892733564013841, 5.363636363636363, 24.529032258064515,
                  102.57741935483871, 290.69032258064516, 479.9451612903226,
                  242.29354838709676, 133.59677419354838, 157.0193548387097, 
                  199.34726688102893, 197.16077170418006, 180.36655948553056,
                  198.62700964630224, 292.4662379421222, 529.2090032154341,
                  495.4855305466238, 349.2829581993569, 249.36334405144694, 
                  184.85530546623795, 138.34405144694534, 88.9967845659164)
count_weekday = c(94.48965517241379, 71.9103448275862, 53.74825174825175,
                  25.53472222222222, 8.544827586206896, 9.373239436619718,
                  19.99310344827586, 47.26896551724138, 112.2551724137931, 
                  177.9241379310345, 263.80689655172415, 325.3862068965517, 
                  379.11034482758623, 387.82068965517243, 378.7310344827586, 
                  373.70344827586206, 367.64827586206894, 339.1241379310345, 
                  292.24827586206897, 242.3448275862069, 183.80689655172415, 
                  148.73793103448276, 123.35172413793103, 90.60689655172413)

# 绘制租赁量和时刻的柱形图
barplot(height = count_workday, xlab = "时刻", 
        ylab = "租赁量", main = "工作日使用情况", names.arg = hours)
#由工作日-时刻租赁柱状图/饼状图可知，工作日的共享单车使用高峰期集中在上午7-10点
#和下午18-20点，与上下班交通高峰期正好重合；
barplot(height = count_weekday, xlab = "时刻",
        ylab = "租赁量", main = "休息日使用情况", names.arg = hours)
# 绘制租赁量和时刻的饼状图
pie(x = count_workday, labels = hours, main = "工作日使用情况", radius = 1)
pie(x = count_weekday, labels = hours, main = "休息日使用情况", radius = 1)
#由休息日-时刻租赁柱状图/饼状图可知，休息日的共享单车使用集中在11-19点；

由工作日-时刻租赁柱状图可知，工作日的共享单车使用高峰期集中在上午 7-10 点和下午 18-20 点，与上下班交通高峰期正好重合。饼状图为工作日不同时间段共享单车使用占比情况。

由休息日-时刻租赁柱状图可知，休息日的共享单车使用集中在 11-19 点。饼状图为休息日不同时间段共享单车使用占比情况

（9）做主成分可视化分析

X = data.frame(
  x1=c(temperature),
  x2=c(atemp),
  x3=c(humidity),
  x4=c(windspeed),
  x5=c(casual),
  x6=c(register)
)
str(X)
cor(X) #计算相关矩阵
#### 作主成分分析
PCA=princomp(X,cor = T)
PCA
PCA$loadings  #主成分负荷
eigen(cor(X)) #计算相关矩阵的特征值
PCA$scores   #主成分得分
summary(PCA,loadings = TRUE)
screeplot(PCA,type='lines')

由于前三个主成分的累积贡献了已达到 82.04% ，所以另外三个主成分可以舍去，达到降维 目的。将温度（ temp ）、最高温度（ atemp ）和湿度（ humidity ）作为三个主要成分，可以 认为共享单车的租赁量主要与这三个环境有关。

3.4作回归诊断图

#### 将主成分作线性回归
lm.sol<-lm(count~temperature+atemp+humidity,data=X)
lm.sol$coefficients  #输出估计的回归系数
summary(lm.sol)
#回归方程通过了回归参数的检验与回归方程的检验，得到回归方程为：
#Y = 178.32869+3.08058*x1+5.26462*x2-2.80534*x3
#回归诊断
influence.measures(lm.sol)
op<-par(mfrow=c(2,2),mar=0.4+c(4,4,1,1),oma=c(0,0,2,0))
plot(lm.sol,1:4)
par(op)

四、讨论和小结：联系实际，分析问题（对应第一部分提出的问题）

1 、用户对共享单车的使用习惯，环境对共享单车运营带来的影响？

由工作日- 时刻租赁柱状图可知，工作日的共享单车使用高峰期集中在上午 7-10 点和下午 18-20 点，与上下班交通高峰期正好重合；由休息日- 时刻租赁柱状图可知，休息日的共享单车使用集中在 11-19 点；温度和湿度都会对共享单车的使用带来影响，由温度- 租赁量散点图可知，当温度在 15-25 之间时，共享单车租赁量较高；当温度在 15 度以下时，温度越低，租赁量越低；当温度较高时，仍能保持一定的租赁量；由湿度- 租赁量散点图可知，湿度在 30-60 区间内，租赁量较多，湿度低于 20 ，且大于 90 ，租赁量较少。

2 、共享单车的租赁量主要与哪些环境有关？

由碎石图可知，由于前三个主成分的累积贡献了已达到 82.04% ，所以另外三个主成分可以舍去，达到降维目的。将温度（temp ）、最高温度（ atemp ）和湿度（ humidity ）作为三个主要成分，可以认为共享单车的租赁量主要与这三个环境有关。

五、源代码

# 读取数据
data <- read.csv("C:\\Users\\leglon\\Desktop\\假期r\\train(1).csv", header = TRUE,stringsAsFactors = F)
#缺失值处理
na.omit(data,cols,invert)
# 从data中读取温度、最高温度、湿度、风速、损坏使用数据、登记使用数据，租赁总数据
suppressWarnings(temperature <- as.numeric(data[,6]))
suppressWarnings(atemp <- as.numeric(data[,7]))
suppressWarnings(humidity <- as.numeric(data[,8]))
suppressWarnings(windspeed <- as.numeric(data[,9]))
suppressWarnings(casual <- as.numeric(data[,10]))
suppressWarnings(register <- as.numeric(data[,11]))
suppressWarnings(count <- as.numeric(data[,12]))

# 温度的均值和标准差
mean(temperature, na.rm = T)
sd(temperature, na.rm = T)

# 湿度的均值和标准差
mean(humidity, na.rm = T)
sd(humidity, na.rm = T)

# 风速的均值和标准差
mean(windspeed, na.rm = T)
sd(windspeed, na.rm = T)

plot(x = humidity,y=count,xlab="湿度",ylab="租赁量",main="湿度与租赁量关系图",cex=1,pch=1)
plot(x=temperature,y=count,xlab="温度",ylab="租赁量",main="温度与租赁量关系图",cex=1,pch=1)
#由温度-租赁量散点图可知，当温度在15-25之间时，共享单车租赁量较高；
#当温度在15度以下时，温度越低，租赁量越低；当温度较高时，仍能保持一定的租赁量；
shapiro.test(temperature[0:5000])
#W的值越接近1就越表明数据和正态分布拟合得越好,可以认为样本数据服从正态分布.
hours = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,22, 23, 24)

count_workday = c(36.73225806451613, 16.003236245954692, 8.436065573770492, 
                  4.892733564013841, 5.363636363636363, 24.529032258064515,
                  102.57741935483871, 290.69032258064516, 479.9451612903226,
                  242.29354838709676, 133.59677419354838, 157.0193548387097, 
                  199.34726688102893, 197.16077170418006, 180.36655948553056,
                  198.62700964630224, 292.4662379421222, 529.2090032154341,
                  495.4855305466238, 349.2829581993569, 249.36334405144694, 
                  184.85530546623795, 138.34405144694534, 88.9967845659164)
count_weekday = c(94.48965517241379, 71.9103448275862, 53.74825174825175,
                  25.53472222222222, 8.544827586206896, 9.373239436619718,
                  19.99310344827586, 47.26896551724138, 112.2551724137931, 
                  177.9241379310345, 263.80689655172415, 325.3862068965517, 
                  379.11034482758623, 387.82068965517243, 378.7310344827586, 
                  373.70344827586206, 367.64827586206894, 339.1241379310345, 
                  292.24827586206897, 242.3448275862069, 183.80689655172415, 
                  148.73793103448276, 123.35172413793103, 90.60689655172413)

# 绘制租赁量和时刻的柱形图
barplot(height = count_workday, xlab = "时刻", 
        ylab = "租赁量", main = "工作日使用情况", names.arg = hours)
#由工作日-时刻租赁柱状图/饼状图可知，工作日的共享单车使用高峰期集中在上午7-10点
#和下午18-20点，与上下班交通高峰期正好重合；
barplot(height = count_weekday, xlab = "时刻",
        ylab = "租赁量", main = "休息日使用情况", names.arg = hours)
# 绘制租赁量和时刻的饼状图
pie(x = count_workday, labels = hours, main = "工作日使用情况", radius = 1)
pie(x = count_weekday, labels = hours, main = "休息日使用情况", radius = 1)
#由休息日-时刻租赁柱状图/饼状图可知，休息日的共享单车使用集中在11-19点；

X = data.frame(
  x1=c(temperature),
  x2=c(atemp),
  x3=c(humidity),
  x4=c(windspeed),
  x5=c(casual),
  x6=c(register)
)
str(X)
cor(X) #计算相关矩阵
#### 作主成分分析
PCA=princomp(X,cor = T)
PCA
PCA$loadings  #主成分负荷
eigen(cor(X)) #计算相关矩阵的特征值
PCA$scores   #主成分得分
summary(PCA,loadings = TRUE)
screeplot(PCA,type='lines')
#z1 = 0.540*X1 + 0.535*X2 - 0.211*X3 - 0.471*X5 + 0.389*X6
#z2 = 0.308*X1 + 0.334*X2 + 0.612*X3 - 0.574*X4 - 0.186*X5 - 0.233*X6
#z3 = 0.301*X1 + 0.273*X2 + 0.684*X4 - 0.286*X5 - 0.531*X6
#由于前三个主成分的累积贡献了已达到82.04%，所以另外三个主成分可以舍去，
#达到降维目的。

#### 将主成分作线性回归
lm.sol<-lm(count~temperature+atemp+humidity,data=X)
lm.sol$coefficients  #输出估计的回归系数
summary(lm.sol)
#回归方程通过了回归参数的检验与回归方程的检验，得到回归方程为：
#Y = 178.32869+3.08058*x1+5.26462*x2-2.80534*x3
#回归诊断
influence.measures(lm.sol)
op<-par(mfrow=c(2,2),mar=0.4+c(4,4,1,1),oma=c(0,0,2,0))
plot(lm.sol,1:4)
par(op)