数据建模与预测期末复习

万瑾琳

信息管理与电子商务系

办公室：7＃305

R语言基础

安装包

install.packages("name")

加载包

Library(name)

变量赋值

a=1; a<-1

向量:

赋值

a=c(1,2); a<-c("xxxx","yyyy")

min(x), max(x), range(x)分别表示向量x的最小分量、最大分量和向量x的范围

与min(x), max(x)有关的函数是which. min(), which. max(),表示在第几个分量求到其值

等差数列

x<-1:30

等差数列运算大于加减乘除运算

等间隔数列

seq(from,to,by = 间隔)

检测是否为空值并替换

s<-is.na(data)

data[is.na(data)]<-num

字符型向量

分隔用的字符可用sep参数指定，如：

paste("result.", 1:4, sep="")

[1] “result.1” “result.2” “result.3” “result.4”

image-20211213192203407

多维数组和矩阵

生成数组矩阵

dim(a)<-c(3,4)
column.names <- c("col1","col2","col3") #定义列名
row.names <- c("row1","row2","row3") #定义行名
dimnames = list(column.names,row.names)
#矩阵行数
nrow(X)
#矩阵列数
ncol(X)

生成数组或矩阵－用matrix()函数构造矩阵

A<- matrix(1:15,nrow=3,ncol=5,byrow=TRUE)

矩阵转置

T()函数

如果矩阵A和B具有相同的维数，则A＊B表示矩阵中对应的元素的乘积，A％＊％B表示通常意义下的两个矩阵的成积（要求矩阵A的列数等于矩阵B的行数

生成对角阵和矩阵取对角阵运算

函数diag()依赖于它的变量：

当v是一个向量时，diag(v)表示以v的元素为对角线元素的对角阵.当M是一个矩阵时，则diag(M)表示的是取M对角线上的元素的向量.

解线性方程组和求矩阵的逆矩阵

#解方程组Ax＝b的解x
A<-t(array(c(1:8,10),dim=c(3,3)))
b<-c(1,1,1)
x<-solve(A,b)

#求逆矩阵
b<-solve(A)

求矩阵的特征值与特征向量

#求矩阵特征根与特征向量
ev<-eigen(A)
ev

ev存放着对称矩阵Sm特征值和特征向量，是由列表形式给出的，其中ev $values是Sm的特征值构成的向量，ev$ vectors是Sm的特征向量构成的矩阵.

image-20211213194050021

矩阵合并

#求矩阵的合并 函数cbind()把其向量纵向合并
x1<-c(1,2,3,4,5,6,7)
x2<-c(2,3,4,6,7,9,10)
cbind(x1,x2)

#函数rbind()把其向量横向合并
x1<-c(1,2,3,4,5,6,7)
x2<-c(2,3,4,6,7,9,10)
rbind(x1,x2)

image-20211213194445173

数据框

第二章 P63

构建数据框

data.frame()

X<-data.frame(
name=c("amy","james","jeffrey","john"),
sex=c("F","M","F","M"), 
age=c(13,12,13,12),
height=c(56.5,57.3,62.5,59.0),
weight=c(84.0,83.0,84.8,99.5))

image-20211213194945928

数据的调用⭐

从剪切板调用

#从剪贴板读取数据
#打开excel数据，用ctrl+C复制数据
#然后调用read.delim函数，从剪贴板读取数据

b <-read.delim("clipboard")

从文本调用

read.table(file, header = FALSE, sep = "")

常用参数说明：

file 文件路径

header 第一行是否为列名称

sep 字段之间的分隔符

quote字符串的标记，默认为"

dec小数点的符号，默认为.

读取excel/csv

X<-read.xlsx("F:\\OneDrive - chenxiangxi\\桌面\\2.xlsx",sheetIndex = 1) #安装xlsx
X<-read.csv("textdata.csv")

储存文件

#存储到csv
write.csv(mydata,"G:/xxx/Inventory.csv")

#存储数据到txt文件,制表符分隔
write.table(mydata,""G:/xxx/Inventory.txt", sep="\t")

#存储数据到xlsx文件
write.xlsx(mydata, "G:/xxx/Inventor..xlsx", sheetName="Sheet1")

多元数据直观表示

均值

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

x是对象，trim是在计算均值前去掉x两端观察值的比例，默认值为0，即包括全部数据。当na. rm＝TRUE时，允许数据中有缺失数据。

求矩阵平均值

用apply()函数

apply(x, 1, mean)  #计算矩阵各行的均值
apply(x, 2, mean)  #计算矩阵各列的均值

方差

var(x)计算样本方差

sd(x)计算样本标准差

变异系数CV

当CV>15%时，说明数据离散程度过高，需要注意

cv<-100*sd(x)/mean(x)

偏度

偏度系数是刻画数据的对称性指标。关于均值对称的数据其偏度系数为0，越趋近于0，数据越服从正态分布

右侧更分散的数据偏度系数为正，左侧更分散的数据偏度系数为负

n<-length(x)
m<-mean(x)
s<-sd(x)
g1<-n/((n-1)*(n-2))*sum((x-m)^3)/s^4

峰度

当数据的总体分布为正态分布时，峰度系数近似为0

当分布较正态分布的尾部更分散时，峰度系数为正；否则为负

当峰度系数为正时，两侧数据较多，当峰度系数为负，两端数据较少

g2<-((n*(n+1))/(n-1)*(n-2)*(n-3))*sum((x-m)^4/s^4-(3*(n-1)^2)/((n-2)*(n-3)))

描述统计函数

summary()
library(Hmisc)  #加载Hmisc包
describe(mydata)  #调用函数

该函数统计样本量、缺失样本量、样本唯一值数、均值、百分之五位数、百分之十位数、四分之一位数、中位数、四分之三位数、百分之九十位数、百分之九十五位数和最小、最大五个数值，分别为：n, nmiss, unique, mean, 5,10,25,50,75,90,95th percentiles, 5 lowest and 5 highest scores

library(psych)      #加载psych安装包
describe(mydata)    #调用describe函数

该函数统计变量名、变量序号、有效样本量、均值、样本标准差、中位数、众数、最小值、最大值、偏度、峰度、标准误差，分别为：item name ,item number, nvalid, mean, sd, median, mad, min, max, skew, kurtosis, se

library(pastecs) 
stat.desc（x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95）

数据框或者时间序列的所有值、空值、缺失值的数量，以及最小值、最大值、值域、总和、均值、中位数、标准误、平均数值信度为95%的置信区间、方差、标准差以及变异系数#nbr.val, nbr.null, nbr.na, min max, range, sum

分组统计函数

library(psych) #加载安装包
describeBy(data, newdata$Gender) #分组统计

频数统计

分析各年龄层所占人数

#绑定数据
attach(data)
#一维列联表
table(年龄)

分析各年龄层性别比例

#二维列联表
table(年龄,性别)
#以年龄、性别排列的结果频数三维列联表
ftable(年龄,性别,结果)

交叉表

频数交叉表

mytable <- table(A, B, C) 
ftable(mytable)

调包实现

library(gmodels)
CrossTable(A,B)

绘制多元分析图

多元分析（下）P13

#安装ggplot包
library(ggplot2)
ggplot(data)+geom_type()

绘制条形图

ggplot(BOD) + geom_bar(aes(x=Time,y=demand),stat="identity")

将x轴由连续值变成不连续值

#factor()函数
BOD$Time<-factor(BOD$Time)
#将BOD中time列变成分类变量，就不再是连续值了

簇状条形图和堆积面积图

簇状条形图和堆积面积图包含两个定性变量（分类变量）以及一个定量变量的情况，用fill函数将不同类别填充为不同的颜色

ggplot(cabbage_exp)+geom_bar(aes(x=Cultivar,y=Weight,fill=date),stat=“identity”)

ggplot(cabbage_exp)+geom_bar(aes(x=Cultivar,y=Weight,fill=date),stat=“identity”, position=“dodge”) #dodge表示并排显示

饼图

ggplot(mydata,aes(x=factor(1),fill=Loan.Purpose))+geom_bar()+coord_polar(theta = "y")

factor(1) 是虚拟x轴的变量，原本的y轴被映射到极坐标的theta，因此y没有实际的意义

散点图

ggplot(data, aes(x=wt, y=mpg)) + geom_point(color="blue")

折线图

多元分析（下）P25

ggplot(data,aes(x=wt,y=mpg)) + geom_line()

分组折线图

ggplot(mtcars,aes(x=wt,y=mpg,group=factor(gear)))+geom_line(aes(color=factor(gear)))
#将gear变成因子变量

频数直方图

'''
binwidth = binsize 用来设置直方图条形的宽度，组距越宽，直方图越粗糙。
#fill="pink" 填充的颜色
colour="blue" 边缘的颜色
'''

mydata <- read.csv("F:/R/data/credit risk.csv")
binsize <- diff(range((mydata$Checking)))/20 #将x取值切分为20组
ggplot(mydata,aes(x=Checking)) + geom_histogram(binwidth = binsize, fill="pink", colour = "blue")

概率密度估计图

ggplot(mtcars,aes(x=wt)) + geom_density()

多元线性回归

读取数据

rdata <- read.csv("F:/Rdata/crime.csv")
library(xlsx)
rdata <- read.xlsx("F:/Rdata/crime.xlsx")

基本描述统计

#汇总分析
summary(rdata)
library(psych)
a=describe(rdata[,c("crime","poverty","single")])
# 散点图 
pairs(~crime+poverty+single,data=rdata, main="Simple Scatterplot Matrix")
#or
ggplot(data, aes(x=wt, y=mpg)) + geom_point(color="blue")

建立回归模型并估计参数

ols<- lm(crime ~ poverty + single, data = rdata))
summary(ols)

image-20211225233930765

对模型拟合效果进行分析

opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(ols, las = 1)

查看偏离较大的数据

P 27

library(MASS)
d1 <- cooks.distance(ols)          #计算偏离程度           
r <- stdres(ols)                   #标准化回归误差
a <- cbind(rdata, d1, r)           #绑定在一起，进行分析
a[d1 > 4/51, ]                     #前4个数据

库克距离：库克距离用来判断强影响点是否为Y的异常值点。一般认为，当D<0.5时认为不是异常值点；当D>0.5时认为是异常值点

在判断Cook距离大小的时候，通常采用的经验分界点是Cook距离序列的4/n处，其中n是观测值的个数。

讨论模型是否可以通过robust回归修正

Huber方法

残差较小的观测值被赋予的权重为1，残差较大的观测值的权重随着残差的增大而递减

library(MASS)
summary(rr.huber <- rlm(crime ~ poverty + single, data = rdata))

也可以查看加权情况

hweights <- data.frame(state = rdata$state, resid = rr.huber$resid, weight = rr.huber$w) 
hweights2 <- hweights[order(rr.huber$w), ] 
hweights2[1:15, ]

Bisquare方法

所有的非0残差所对应观测值的权重都是递减的

rr.bisquare <- rlm(crime ~ poverty + single, data=rdata, psi = psi.bisquare) 
summary(rr.bisquare)

也可以查看加权情况

biweights <- data.frame(state = cdata$state, resid = rr.bisquare$resid, weight = rr.bisquare$w) 
biweights2 <- biweights[order(rr.bisquare$w), ]
biweights2[1:15, ]

Huber方法的软肋在于无法很好的而处理极端离群点，而bisquare方法的软肋在于回归结果不易收敛，以至于经常有多个最优解

逐步回归

#先建立线性回归模型
steplm<-lm(crime~pctmetro+pctwhite+pcths+poverty+single,data=rdata)
#逐步回归
step(steplm, direction = c("both", "backward", "forward"))

the mode of stepwise search, can be one of "both", "backward", or "forward", with a default of "both". If the scope argument is missing the default for direction is "backward".

logitstic回归

函数形式

$logit(y)=ln\frac{P}{1-P}=\beta_0+\beta_1x_1+\beta_2x_2+...+\beta_px_p=X\beta$

logistic回归系数的意义

p 20 例二全面一点p 31

读取数据

d1=read.table(“clipboard”,header=T)  #读取例1数据

描述统计

同上

建立全变量logistic回归模型

logit.glm<-glm(y~x1+x2+x3,family=binomial,data=d1)  #Logistic回归模型
summary(logit.glm)  #Logistic回归模型结果

AIC为模型拟合效果

得到初步的logistic回归模型：
$p=\frac{exp(0.5976-1.4961x_1-0.0016x_2+0.3159x_3)}{1+exp(0.5976-1.4961x_1-0.0016x_2+0.3159x_3)}$

$logit(p)=0.5976-1.4961x_1-0.0016x_2+0.3159x_3$

逐步筛选变量logistic回归模型：

logit.step<-step(logit.glm,direction="both")  #逐步筛选法变量选择
summary(logit.step) #逐步筛选法变量选择结果

得到最佳逻辑回归模型

预测

pre1<-predict(logit.step, data.frame(x1=1))  #预测视力正常司机Logistic回归结果
p1<-exp(pre1)/(1+exp(pre1))  #预测视力正常司机发生事故概率
pre2<-predict(logit.step,data.frame(x1=0))  #预测视力有问题的司机Logistic回归结果
p2<-exp(pre2)/(1+exp(pre2))  #预测视力有问题的司机发生事故概率
c(p1,p2)  #结果显示

多类别logitstic回归

需要r包

library(foreign) 
library(nnet) 
library(ggplot2)
library(reshape2)

读取数据

同上

描述统计

同上

设定参照类，做多分类logistic回归分析

Library(nnet)
#设定参照类,新生成一列
ml$prog2 <- relevel(ml$prog, ref = "academic")
#做多分类logistic回归
test <- multinom(prog2 ~ ses + write, data = ml)

获取多分类logistic回归分析结果

summary(test)

image-20211227223556328

ln\left(\frac{P(prog=general)}{P(prog=academic)}\right) = b_{10} + b_{11}(ses=2) + b_{12}(ses=3) + b_{13}write

$ln\left(\frac{P(prog=vocation)}{P(prog=academic)}\right) = b_{20} + b_{21}(ses=2) + b_{22}(ses=3) + b_{23}write$

b13 变量 write 增加一个单位与参加普通课程与学术课程的对数几率降低 0.058 相关。
b23 变量 write 增加一个单位与参加职业计划与学术计划的对数几率降低有关。 0.1136 的数量。
b12 如果从 ses="low" 变为 ses="high"，则在普通课程与学术课程的对数几率将降低 1.163。
b11 如果从 ses="low" 移动到 ses="middle"，则普通课程与学术课程的对数几率将降低 0.533，尽管该系数不显着。
b22 如果从 ses="low" 变为 ses="high"，则参加职业计划与学术计划的对数几率将降低 0.983。
b21 如果从 ses="low" 移动到 ses="middle"，职业计划与学术计划的对数几率将增加 0.291，尽管该系数不显着。

计算标准化得分 z-score及p-value

z <- summary(test)$coefficients/summary(test)$standard.errors 
z
p <- (1 - pnorm(abs(z), 0, 1)) * 2 
p

基本原理如下：

原假设 H0：βj=0

备择假设 H1：βj≠0

构造统计量：
$z_{score}=\frac{\bar{\beta}_j-0}{std.errors}$
计算p值，即z_score发生的概率，双侧检验

提取参数估计值并转换为数值

exp(coef(test))

image-20211227225400405

ses社会地位变量从1变为3时，普通项目vs学术项目的相对优势比是0.3126

Write变量增加一个单位，普通项目vs学术项目的相对优势比是0.9437

查看各样本所属分类的概率

head(pp <- fitted(test))

预测

dses <- data.frame(ses = c("low", "middle","high"),write=mean(ml$write))
predict(test, newdata = dses, "probs")
#预测：Write不变，检测社会地位的改变对预测值的影响

泊松回归

研究因变量为计数变量，自变量为数值变量或分类变量之间关系的回归模型

require(ggplot2) 
require(sandwich)

读取数据

p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")

描述统计

同上

p <- within(p, { prog <- factor(prog, levels=1:3, labels=c("General", "Academic", "Vocational"))id <- factor(id)  })
summary(p)
with(p, tapply(num_awards, prog, function(x) {sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))}))

画图

ggplot(p, aes(num_awards, fill = prog))+geom_histogram(binwidth=.5, position="dodge")

建立回归模型

summary(m1 <- glm(num_awards ~ prog + math, family="poisson", data=p))

计算p值及置信区间

library(sandwich)
cov.m1 <- vcovHC(m1, type="HC0")
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(Estimate= coef(m1), "Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)
r.est

模型的拟合优度，通过卡方统计量来判断模型的拟合优度

with(m1, cbind(res.deviance = deviance, df = df.residual,p = pchisq(deviance, df.residual, lower.tail=FALSE)))

预测

p$phat <- predict(m1, type="response")
# order by program and then by math
p <- p[with(p, order(prog, math)), ]

主成分分析

导入数据计算相关矩阵

X=read.table("clipboard",header=T) 
cor(X)

求相关矩阵的特征值和主成分负荷

PCA = princomp(X,cor=T)  #主成分分析
PCA                  #特征值开根号结果 
summary(PCA)
PCA$loadings  #主成分载荷

最后一行为累计贡献率，一般取累计贡献率达到85%~95%的特征值。

image-20211228201428843

从载荷矩阵可以看出：

Comp1=-0.353*x1-0.249*x2-0.374*x3-0.302*x4-0.376*x5-0.404*x6-0.371*x7-0.374*x8

依次可以写出每个主成分的构成形式

确定主成分

screeplot(PCA,type="lines") #绘制碎石图

按照累积方差贡献率大于80%原则，选入了两个主成分，其累积方差贡献率为80.7%，从碎石图上也可以看出m取2比较合适。

主成分得分

PCA$scores[,1:2]  #主成分得分

image-20211228201658597

主成分分析与多元线性回归结合

P 32

导入数据

conomy<-data

同上

建立多元回归模型

lm.sol<-lm(y~x1+x2+x3, data=conomy)
summary(lm.sol)

主成分分析

b<-data.frame(conomy[,1:3])  #对三个自变量作主成分分析
conomy.pr<-princomp(b, cor=T)
summary(conomy.pr, loadings=TRUE)

image-20211228205605313

将原始数据准换为用主成分变量表达

上一步可知主成分为comp 1和comp 2，将第1主成分的预测值和第2主成分的预测值存放在数据框里

pre<-predict(conomy.pr)
conomy$z1<-pre[,1]
conomy$z2<-pre[,2]

作主成分回归

lm.sol<-lm(y~z1+z2, data=conomy)
summary(lm.sol)

$Y=21.8909+2.994z_1-0.7616z_2$

这是因变量与主成分的关系，应用起来不方便，需变换成因变量与自变量的关系

变换公式代码：

beta<-coef(lm.sol)       #提取回归系数
A<-loadings(conomy.pr)     #提取主成分对应的特征向量
x.m<-conomy.pr$center        #数据的中心,即各类数据均值
x.sd<-conomy.pr$scale       #数据的标准差         
coef<-(beta[2]*A[,1]+beta[3]*A[,2])/x.sd #
beta0<-beta[1]-sum(x.m*coef)
c(beta0,coef)

因子分析

P 26

读入数据

X=read.table("clipboard",header=T)

数据标准化

data<-data.frame(X[,3:9])
X<-scale(data)

计算相关系数矩阵

cor(X)

计算特征值、因子载荷及共同度

(FA0=factanal(X,3,rot="none")) #极大似然法因子分析

library(mvstats) 
(Fac=factpc(X,3)) #主成份法因子分析

因子旋转

如果求出主因子后，各个主因子的典型代表变量不很突出，还需要进行因子旋转，通过适当的旋转得到比较满意的主因子。

进行因子旋转，就是要使因子载荷矩阵中因子载荷的绝对值向0和1两个方向分化，使大的载荷更大，小的载荷更小

(Fa1=factanal(X,3,rot="varimax")) #varimax法旋转因子分析

image-20211228213917398

因子得分

#使用回归估计法的极大似然法因子分析 
Fa1=factanal(X,3,scores="regression") 
Fa1$scores

或

#使用回归估计法的主成份法因子分析 
Fac1=factpc(X,3,scores="regression") 
Fac1$scores

数据建模与预测_R语言基础

数据建模与预测期末复习

R语言基础

多维数组和矩阵

数据框

数据的调用⭐

多元数据直观表示

均值

方差

变异系数CV

偏度

峰度

描述统计函数

分组统计函数

频数统计

交叉表

绘制多元分析图

绘制条形图

簇状条形图和堆积面积图

饼图

散点图

折线图

频数直方图

概率密度估计图

多元线性回归

读取数据

基本描述统计

建立回归模型并估计参数

对模型拟合效果进行分析

查看偏离较大的数据

讨论模型是否可以通过robust回归修正

Huber方法

Bisquare方法

逐步回归

logitstic回归

函数形式

logistic回归系数的意义

读取数据

描述统计

建立全变量logistic回归模型

逐步筛选变量logistic回归模型：

预测

多类别logitstic回归

读取数据

描述统计

设定参照类，做多分类logistic回归分析

获取多分类logistic回归分析结果

计算标准化得分 z-score及p-value

提取参数估计值并转换为数值

查看各样本所属分类的概率

预测

泊松回归

读取数据

描述统计

建立回归模型

计算p值及置信区间

模型的拟合优度，通过卡方统计量来判断模型的拟合优度

预测

主成分分析

导入数据 计算相关矩阵

求相关矩阵的特征值和主成分负荷

确定主成分

主成分得分

主成分分析与多元线性回归结合

导入数据

建立多元回归模型

主成分分析

将原始数据准换为用主成分变量表达

作主成分回归

因子分析

读入数据

数据标准化

计算相关系数矩阵

计算特征值、因子载荷及共同度

因子旋转

因子得分

导入数据计算相关矩阵