- 利用统计学方法进行：描述统计（利用图表形式对数据进行汇总和展示，包括计算一些统计量如均值、方差、比率等；描述统计包括：图表描述，统计量描述）和推断统计（主要利用样本信息推断总体特征，主要包含两类一类是参数估计（利用样本信息推断所关心的总体特征），一类是假设检验（利用样本信息推断总体的某个假设是否成立））。
-
变量：描述所观察事物的某种特征；
类别变量（定性变量）：无序类别变量，有序类别变量
数值变量（定量变量）：离散变量，连续变量
-
数据：变量的观测值。
数据来源：二手数据，亲自调查，概率抽样，简单随机抽样，分层抽样，系统抽样（等
距抽样），整群抽样。
- R的基本数据类型：
数值型：a<-c(2,5,8,3,9)
字符型：b<-c("甲","乙","丙","丁")
逻辑型：c<-c("TRUE","FALSE","FALSE","TRUE")
因子型
- R的基本数据结构：
向量：一些创建向量的方法：
v1<-1:6 # 产生1~6的等差数列
v2<-seq(from=2,to=4,by=0.5) # 在2~4之间产生步长为0.5的等差数列
v3<-rep(1:3,times=3) # 将1~3的向量重复3次
v4<-rep(1:3,each=3) # 将1~3的向量中每个元素重复3次
矩阵：
a<-1:6 # 生成1到6的数值向量
mat<-matrix(a, nrow=2, ncol=3, byrow=TRUE)
rownames(mat)=c("甲","乙") # 添加行名
colnames(mat)=c("A","B","C") # 添加行名
# 用t函数对矩阵转置
t(mat) # 矩阵转置
数组
数据框：
# 写入姓名和分数向量
names<-c("刘文涛","王宇翔","田思雨","徐丽娜","丁文斌") # 写入姓名向量
stat<-c(68,85,74,88,63) #写入各门课程分数向量
math<-c(85,91,74,100,82)
econ<-c(84,63,61,49,89)
# 将向量组织成数据框形式
table1_1<-data.frame(学生姓名=names,统计学=stat,数学=math,经济学=econ)
因子 factor
列表
- R数据查看：head(t,3), tail(t,3), nrow(t), ncol(t) dim(t), class(t), str(t)#查看t的数据结构
- 选中数据框中特定列：mean(df$统计学), mean(df[,1])
- 读取R数据：load("t.RData")
- 保存为R数据：save(matrix1,file="./data/mat.RData")
- 读取含标题csv文件：t<-read_csv("1.csv")
- 读取不含标题的csv文件：t<-read_csv("1.csv",header=FALSE)
- 保存csv文件：write.csv(t,file="1.csv")
- 读取xlsx文件：t<-read_excel("1.xlsx",as.data.frame=TRUE)
- 保存.xlsx文件：write_xlsx(x=t,"./data/1.xlsx")
- 数据排序:table=Table[order(Table$数学,decreasing=TRUE),]#降序
- table=Table[order(-Table$数学),]#降序
- table=Table[order(Table$数学),]#升序
- 对矩阵行或者列求和
- colSums(t[,1:5])
- rowSums(t)#对t的每一行求和
- rbind(t,totals=colSums(t[,1:5]))#将1-5列求和的结果命名为totals并以行添加到t中
- cbind(t,totals=rowSums(t[1:10,]))#将1-10行的求和结果以列添加到t中
-
apply(matrix1_1, 1, sum) # 对矩阵所有行求和
-
apply(matrix1_1, 2, sum) # 对矩阵所有列求和
-
apply(matrix1_1, 1, mean)
-
apply(matrix1_1, 2, mean)
-
apply(matrix_1, 2, sd)
- 编辑数据框——变量重命名
-
rename(table1_1, c(姓名="name",统计学="stat"))
- 缺失值处理
- x=c(1,2,NA)
- sum(x,na.rm=TRUE)
-
- 数据合并
- rbind(df1,df2)
- cbind(df1,df2[2:5])#按列合并数据
-
短格式->长格式
-
library(reshape2) # 加载reshape2包
-
tab.long<-melt(table1_1,id.vars="姓名",variable.name="课程",value.name="分数")# 融合table1_1与id变量，并命名variable.name="课程",value.name="分数"
-
-
·生成随机数
-
rnorm(n, mean=0, sd=1)
-
#产生n个服从均值为mean，方差为sd正态分布的随机数
-
-
runif(n, min, max) #产生n个在min~max之间服从均匀分布随机数
rexp(n) #产生n个服从指数分布的随机数
rchisq(n, df) #产生n个服从自由度为df的卡方分布的随机数
mytable1<-table(example1_1$社区);mytable1
# 将频数表转化成百分比表
prop.table(mytable1)*100
# 生成社区和态度的二维列联表
mytable2<-table(example1_1$态度,example1_1$社区);mytable2
addmargins(mytable2) # 为列联表添加边际和
addmargins(prop.table(mytable2)*100) # 将列联表转换成百分比表
抽样与数据筛选
简单随机抽样
samples(table1,10) #默认无放回
samples(table1,10,replace=TRUE) #有放回抽样
samples(table1$姓名[table1$成绩>90])#筛选出成绩大于90分的姓名
分层抽样
result=Strata(table1,stratanames="姓名",size=c(round(lenght(table1$性别[table1$性别==“男”])*0.2)，round(lenght(table1$性别[table1$性别=="女"])*0.2)）,method="srswor") #按照性别分层抽取20%样本
系统抽样
result=sampleBy(~1,frac=0.2,relpace=FALSE,data=table1,systematic=TRUE)
#~1表示数据没有分组，frac表述数据抽样的间隔1/0.2=5个，replace=FALSE表示无放回抽样，systematic=TRUE表示做系统抽样。
•将列联表转换为原始数据
library(DescTools)
mytable<-ftable(example1_1,row.vars=c("社区"), col.vars=c("性别","态度"))
# 将列联表转化成原始数据框
df<-Untable(mytable)
head(df,3)
# 将列联表转换成带有交叉频数标签的数据框
dff<-as.data.frame(mytable)
数值数据如何生成生成频数分布表
确定要分组的组数：大致计算出组数。
确定各组的组距
统计出各组的频数
library(DescTools)
example1_2<-read.csv(".//data//example//chap01//example1_2.csv")
tab<-Freq(example1_2$销售额); tab
或者
# 指定组距，不含上限值
tab1<-Freq(example1_2$销售额,breaks=c(160,175,190,205,220,235,250,265,280),
right=FALSE)
tab2<-data.frame(分组=tab1$level,频数=tab1$freq,
频数百分比=tab1$perc*100,累积频数=tab1$sumfreq,
累积百分比=tab1$sumperc*100) #重新命名频数表中的变量
基本绘图函数
| plot() | 散点图，条形图，箱线图。。。。。。 |
| barplot() | 条形图 |
| boxplot() | 箱线图 |
| coplot() | 条件图 |
| dotplot() | 点图 |
| hist() | 直方图 |
| mosaicplot() | 马赛克图 |
| pairs() | 散点图矩阵 |
| abline() | 为图形添加解决为a，斜率为b的直线 |
| lines() | 在坐标点(x,y)之间添加直线 |
| layout() | 布局图形 |
| legend() | 添加图例 |
| mtext() | 在图形区域的边距或区域的外部边距添加文本 |
| text() | 在坐标点(x,y)处添加文本 |
| points() | 在坐标点(x,y)处添加点 |
| title() | 为图形添加标题 |
| 函数 | 数据类型 | 图形 |
| plot() | （一个）数值 | 散点图 |
| plot() | （一个）因子 | 条形图 |
| plot() | 一维列联表 | 条形图 |
| plot() | 数值 vs 数值 | 散点图 |
| plot() | 二维列联表 | 马赛克图 |
| plot() | 数值 vs 因子 | 箱线图 |
| plot() | 因子 vs 数值 | 脊形图 |
| plot() | 因子 vs 因子 | 马赛克 |
绘图
##绘制2*2的图，第二列为1个图
#设置布局环境，2*2，宽度比为2:1
dev.new()
layout(matrix(c(1,2,3,3),nrow=2,ncol=2),widths=c(2,1))
#设置图形边框，文字和绘图符号大小
par(mai=c(0.6,0.6,0.1,0.1),cex=0.7)
#mai表示设置图形边距大小c(底部，左侧，上方，右侧)，cex控制文字和绘图符号的大小
x=rnorm(5000) #标准正态分部
y=rchisq(5000,10) #卡方分布
hist(x,prob=TRUE,col="lightblue",xlab="x",ylab="Density",ylim=c(0,0.4),main="")
hist(y,freq=TRUE,col="pink",xlab="y",ylab="Density",main="")
boxplot(x,col="red",lwd=1)title(main=“My Title”，sub=“sub-title”,xlab=“…”,ylab=“…”
axis( side=, at=, labels=,pos=, lty=, col=, las=, tck=,……)
text(location, “text to place”, pos, ………) #图里面的点添加文字
mtext(“text to place”,side, line=n ,………) #图外面的边框添加文字
legend(location=, title=, legend,…………)
直方图：
用矩形的宽度和高度来表示频数分布，通过直方图可以观察出数据的大致分布情况如是否对称。
hist(table$销售额,xlab="x",ylab="y",main="销售额")
hist(table$销售额, breaks=20,xlab="x",ylab="y",main="销售额(分成20组)")
hist(table$销售额,prob=TRUE,breaks=20,xlab="x",ylab="y",main="销售额分成20组加上核密度线")
lines(density(table$销售额),col="red")
hist(table$销售额,prob=TRUE,breaks=20,xlab="x",ylab="y",main="销售额分成20组加上正态分布曲线")
curve(dnorm(x,mean(table$销售额),sd(table$销售额),add=T,col="red"))
茎叶图
反映数据的分布特征，且能保留原始数据的基本信息。，当数据较少时，用茎叶图更容易观察数据的分布。
箱线图
反映一组数据的分布特征，还可以用于多组数据分布特征进行比较。
小提琴图
（将分布的核密度曲线和与箱线图结合在一起，在箱线图上以镜像的形式叠加一条核密度估计曲线）从核密度估计曲线可以大致看出数据的形状。
三、数据的描述统计量
描述水平的统计量——平均数
能够反映数据的总体特征
· 消除了观测值的随机波动，容易受到极端值的影响。
由总体计算出来的均值称为总体均值(mu)，用样本数据计算出来的均值称为样本均值(x_bar)
计算普通均值：mean(table$成绩)
计算加权均值：weighted.mean(table$组中值,table$人数)，公式中mi表示组的平均数，fi表示频数。
描述水平的统计量——分位数
将数据按从小到大排序后，位于某一位置上的数值，分位数有：中位数，（上/下）四分位数，百分位数。分位数不受极端值的影响，在研究收入分配时很有用。
中位数计算： median median(table$销售额)
25%位置为下四分位数，75%称为上四分位数
四分位数位置的计算：
Q25%的位置=（n+1）/4 , Q75%的位置为=3（n+1）/4,当计算出的位置不是整数时，若小数位是0.5则取两边数值的平均数作为分位数；若是0.25则应在整数位置（m）的数值上加0.25*（值m+1-值m）。 quantile(table1$销售额,probs=c(0.25,0.75),type=6)
描述水平的统计量——百分位数
将数据用99点均分为100份，处于各个位置的数值就是对应的分位数；反映了各项数据在最大数据与最小数据之间的分布关系。
描述差异的统计量
差异量反映了数据的离散程度，描述差异的统计量值越小说明描述水平的统计量数据的代表性就越好，反之说明描述水平的统计量对该数据的代表性就越差。
描述差异的统计量——极差（全距）
最大值与最小值之差，容易受到极端值的影响，不能全面反映差异情况。
描述差异的统计量——四分位差
上四分位数-下四分位数，反映了中间50%数据的离散程度，能够说明中间数据的集中情况，不受极端值的影响。
描述差异的统计量——平均离差（数据与均值之差称为离差）
将离差取取绝对值后在平均
描述差异的统计量——方差(var)
将离差平方取平均后在平均（分母是n-1）
描述差异的统计量——标准差(sd)
将方差开方，标准差与原始数据的单位相同，其实际意义要比方差清楚，因此实际问题时常常使用标准差。
方差和标准差存在的问题: 当数据的观测值越大时，方差和标准差的值也会越大，采用不同的计量单位计算出的方差和标准差不一样，对于不同样本的数据，若数据的观测值相差较大或者计量单位不同时不能使用标准差直接比较离散程度。
描述差异的统计量——变异系数 CV
数据 标准差/均值，主要用于比较不同样本数据的离散程度，变异系数越小说明数据的相对离散程度越小，数值越大说明数据的相对立离散程度越大；消除了数值大小和计量单位对标准差的影响，反映了一组数据的相对离散程度。
描述差异的统计量——标准分数 z
也称z分数或者标准化值，用于度量数值在数据中的相对位置，并判断一组数据是否有离群点,可以把数据转化均值为0，标准差为1的新数据。计算方法为：zi=(xi-x_bar)/s （数据与平均值相差几倍方差）
描述分布形状的统计量——偏度系数 SK （skewness(table1$销售额)）
偏度是指数据分布的不对称性，测量数据的分布不对称性的统计量称为偏度系数SK，若取值大于1或者小于-1则称为严重偏斜分布，若-1<SK<-0.5或者0.5<SK<1则称为中等偏斜分布，若-0.5<SK<0.5则称为轻微偏斜。取负值表示左偏分布（左侧有长尾），取正值表示右偏分布（右侧有长尾）；越接近于零表示越接近对称分布。
描述分布形状的统计量——峰度系数 K （kurtosis(table1$销售额)）
峰度是指数据分布峰值的高低(通常与标准正态分布相比较而言)，测度一组数据分布峰值高低的统计量称为峰度系数。当K大于零时为尖峰分布，数据分布的峰值比标准正态分布高，数据相对集中；K<0时为扁平分布，数据分布的峰值比标准正态分布低，数据较为分散。
四、随机变量的概率分布
在n次伯努利实验中成功x次的概率
标准差是根据原始观测值计算出来的，反映了原始数据的离散程度，而标准误是根据样本统计量计算出来的，反映的是统计量的离散程度。标准误越小，表名样本统计量与总体参数的值越接近，样本对总体越具有代表性。
中心极限定理：如果总体不是正态分布，随着样本量n的增大，样本均值的概率分布仍趋于正态分布，其分布的期望值为总体均值mu，方差为总体方差的1/n。
五、参数估计
#总体均值的置信区间(大样本)
setwd("E:\\大三下学期\\R\\R_code\\data\\data\\")
t1=read.csv(".\\example\\chap05\\example5_1.csv")
install.packages("BSDA")
install.packages("BSDA")
library(BSDA)
z.test(t1$耗油量,mu=0,sigma.x=sd(t1$耗油量),conf.level=0.9)
#一个总体均值的置信区间（小样本,方差未知）
t1=read.csv(".\\example\\chap05\\example5_2.csv")
t.test(t1,conf.level=0.9)
#计算两个总均值之差的置信区间（大样本）
t1=read.csv(".\\example\\chap05\\example5_3.csv")
z.test(t1$男性工资,t1$女性工资,mu=0,sigma.x=sd(t1$男性工资),sigma.y=sd(t1$女性工资),conf.level=0.95)$conf.int
#上面只显示95置信水平下的置信区间
#计算两个总体均值之差的置信区间（小样本）
#假设方差相等
t1=read.csv(".\\example\\chap05\\example5_4.csv")
t.test(t1$方法一,t1$方法二,var.equal=TRUE,conf.level=0.95)
#假设方差不等
t.test(t1$方法一,t1$方法二,var.equal=FALSE,conf.level=0.95)
#配对样本估计（大样本用z检验,小样本用t检验）
t1=read.csv(".\\example\\chap05\\example5_5.csv")
t.test(t1$试卷A,t1$试卷B,paired=TRUE)
#z.test(t1$差值d,mu=mean(t1$差值d),sigma.x=sd(t1$差值d))
#z检验书上没有，应该是这样写，不用记这一个
#总体比例的区间估计(大样本)
#某城市想要进行一项交通措施改革，为征求市民对该项改革措施的意见，
#在成年人中随机调查了500个市民，其中325人赞成改革措施。用97.5%的
#置信水平估计该城市成年人口中赞成该项改革的人数比例的置信区间
n=500
x=325
p=x/n
q=qnorm(0.975)
#95%的置信水平，则alpha=0.05,那么在计算分位数时应该
#计算1-alpha/2位置的分位数，也就是0.975
LCI=p-q*sqrt(p*(1-p)/n)
UCI=p+q*sqrt(p*(1-p)/n)
data.frame(LCI,UCI)
#结果是：0.6081925 0.6918075
#任意大小样本的比例估计（最重要的是n,p的更换后的定义）
n=500+4
x=325+2
p=x/n
q=qnorm(0.975)
#95%的置信水平，则alpha=0.05,那么在计算分位数时应该
#计算1-alpha/2位置的分位数，也就是0.975
LCI=p-q*sqrt(p*(1-p)/n)
UCI=p+q*sqrt(p*(1-p)/n)
data.frame(LCI,UCI)
#结果是：0.6071358 0.6904833
#两个总体比例的置信区间估计
#在某个电视节目的收视率调查中，女性观众随机调查了500人，
#有225人收看了该节目；男性观众随机调查了400人，有128人
#收看了该节目。用95%的置信水平估计女性与男性收视率差值的置信区间
#利用大样本的估计方法
n1=500
n2=400
x1=225
x2=128
p1=x1/n1
p2=x2/n2
q=qnorm(0.975)
LCI=p1-p2-q*sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
UCI=p1-p2+q*sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
data.frame(LCI,UCI)
#结果是：0.06682346 0.1931765，
#需要注意在计算两个总体的上下置信区间时需要p1-p2
#利用任意样本的区间估计
n1=500+2
n2=400+2
x1=225+1
x2=128+1
p1=x1/n1
p2=x2/n2
q=qnorm(0.975)
LCI=p1-p2-q*sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
UCI=p1-p2+q*sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
data.frame(LCI,UCI)
#结果是：0.06624396 0.1923634#特别注意在比例区间估计中q的计算是1-alpha/2
六、假设检验
概率P小于显著性水平则拒绝H0
单样本t检验效应量的计算公式为：d=|样本均值-假设的总体均值|/样本标准差,0.2<d<0.5是小的效应量，0.5<d<0.8是中的效应量，d>0.8是大的效应量
library(lsr)
cohensD(t1$厚度,mu=5)#一个总体的效应量
cohensD(t1$甲企业,t1$乙企业)#两个独立样本t检验的效应量
cohensD(t1$旧款饮料,t1$新款饮料,paired=TRUE)#配对样本t检验的效应量
#一个总体均值检验（大样本）:当备择检验是mu小于81时
library(BSDA)
z.test(t1$数据,mu=81,sigma.x=sd(t1$数据),alternative="less",conf.level=0.95)
#因为备择假设是小于，所以用alternative="less"
#一个总体均值检验(小样本)
t.test(t1$数据,mu=81)
library(lsr)
cohensD(t1$数据1,mu=81)
#两个总体均值之差检验（大样本），当备择假设是不等于号时
z.test(t1$数据1,t1$数据2,sigma.x=sd(t1$数据1),sigma.y=sd(t1$数据2),alternative="two.sided")
#两个总体均值之差检验（小样本），方差相等
t.test(t1$数据1,t1$数据2,var.equal=TRUE)
#两个总体均值之差检验（小样本），方差不相等
t.test(t1$数据1,t1$数据2,var.equal=FALSE)
#效应量计算
library(lsr)
cohensD(t1$数据1,,t1$数据2)
#配对样本检验
t.test(t1$数据1,t1$数据2,paired=TRUE)
#计算效应量
library(lsr)
cohensD(t1$数据1,t1$数据2,method="paired")
#一个总体比例的检验
#备择假设是pai>0.25
n=2000
p=450/2000
pi0=0.25
z=(p-pi0)/sqrt(pi0*(1-pi0)/n)
p_value=1-pnorm(z) #因为备择假设是>0.25
data.frame(z,p_value)
#两个总体比例之差的检验（pi1-pi2=0）
n1=200
n2=200
p1=0.27
p2=0.35
p=(p1*n1+p2*n2)/(n1+n2)
z=(p1-p2)/sqrt(p*(1-p)*(1/n1+1/n2))
p_value=pnorm(z)
data.frame(z,p_value)
#两个总体比例之差的检验（pi1-pi2=d0）
n1=300
n2=300
p1=33/300
p2=84/300
z=(p1-p2-d0)/sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
p_value=pnorm(z)
data.frame(z,p_value)
#一个总体的方差检验
#假设备择假设是大于16
sigma.test(t1$数据,sigmasq=16,alternative="greater")
#两个总体的方差检验
#假设备择假设是不等于
var.test(t1$数据1,t1$数据2,alternative="two.sided")
八、方差分析
方差分析是分析类别自变量对数值因变量影响的一种统计方法，（自变量对因变量的影响简称效应），方差分析是通过对数据误差的分析来检验这种效应是否显著,分析数据处理误差是否显著存在。
#第八章代码
##实验8.2:单因子方差分析
#假设不同管理者对评分影响的效应量分别为alpha1,alpha2,alpha3
#H0:alpha1=alpha2=alpha3
#H1：alpha1!=alpha2!=alpha3
setwd("E:\\大三下学期\\R\\R_code\\data\\data\\")
T1<-read.csv(".\\exercise\\chap08\\exercise8_2.csv")
#将短格式数据改为长格式
library(reshape)
T2=melt(T1,variable_names="评分")
T2=rename(T2,c(variable="管理者水平",value="评分"))#这个语句要注意
#(1)方差分析
model1=aov(评分~管理者水平,data=T2)#注意这一项数值在前，类别变量在后
summary(model1)
#P=0.000633,因此有理由说明管理者水平对评分有影响
#计算效应量
library(DescTools)
EtaSq(model1,anova=T)
#结果表明评分总误差中被管理者水平解释的比例为62.54%
#(2)LSD,HSD方法比较管理者之间评分的差异,找出mu1,mu2,mu3之间到底哪两个均值不等，
#这种对均值的配对检验称之为方差分析中的多重比较。
library(agricolae)
LSD=LSD.test(model1,"管理者水平")
LSD
HSD=HSD.test(model1,"管理者水平")
HSD
#(3)正态性检验，QQ图，或者shapior-Wilk,K-S检验
qqnorm(T2$评分,xlab=c("观察值"),ylab=c("期望正态值"),main="QQ图")
qqline(T2$评分)
#shapiro检验正态性，P大于0.05则说明服从正态分布
shapiro.test(T2$评分)
#ks检验正态性
ks.test(T2$评分,"pnorm",mean(T2$评分),sd(T2$评分))
#方差齐性检验:画图，Bartlett方差齐性检验
plot(model1,which=c(1,2))
bartlett.test(评分~管理者水平,data=T2)
#也可以写成barlett(model1)
#P=0.3726>0.05表示不同管理者水平的评分满足方差齐性。
#8.5双因子方差分析
setwd("E:\\大三下学期\\R\\R_code\\data\\data\\")
T1<-read.csv(".\\exercise\\chap08\\exercise8_5.csv")
library(reshape)
T2=melt(T1,variable_name="路段")#将类别变量那一列进行命名
T2=rename(T2,c(value="车流量"))
#主效应方差分析
model2=aov(车流量~时段+路段,data=T2)
summary(model2)
#两个P值都小于0.05说明时段和路段对车流量有显著影响
#主效应方差分析的参数估计
model2$coefficients
#交互效应量分析
#设时段对车流量的影响效应为alpha1,alpha2，路段对车流量的影响效应为beta1,beta2,beta3,
#路段和时段对车流量影响的交互效应为gama_ij
#检验路段效应的假设为：H0：alpha1=alpha2=0,H1:alpha1,alpha2至少有一个不等于0
#检验时段效应的假设为：H0：bata1=beta2=beta3=0,H1:beta1，beta2，beta3至少有一个不等于0
#检验交互效应的假设为：H0:gama_ij=0（i=1,2,j=1,2,3）,H1:gama_ij至少有一个不等于0
model3=aov(车流量~时段+路段+时段:路段,data=T2)
summary(model3)
#时段：路段的P值>0.05，说明时段和路段的交互效应不显著
#计算交互效应方差的效应量
library(DescTools)
EtaSq(model3)
#看第二列数据，在排除时段和路段交互影响后，时段对车流量的误差的解释比例为74.5%
#排除时段和路段交互影响后，路段解释了车流量误差的比例为64.7%
#排除路段和时段的影响后，路段和时段的交互效应解释了车流量误差的0.02%,
#如此小的偏效应量也证明了前面交互效应对车流量没有显著影响的结论
九、一元线性回归
#第九章
#线性回归建模的基本步骤
#1.检验变量之间的关系
#2.建立线性关系模型，并进行估计和检验
#3.用建立好的线性模型进行预测
#4.对模型进行诊断（包括正态性检验，ks,shapiro,plot和方差齐性检验ncvTest,plot(位置尺图)）
setwd("E:\\大三下学期\\R\\R_code\\data\\data\\")
T1<-read.csv(".\\example\\chap09\\example9_1.csv")
#确定变量之间的关系
cor.test(T1$销售收入,T1$广告支出)
#p-value = 1.161e-09<0.05说明两个变量具有较强的线性相关性
model=lm(销售收入~广告支出,data=T1)
summary(model)
#看coefficients这一列的estimate参数，第一个为截距，第二个为斜率,
#同时第二行的t值和p_value反映了x对y的影响程度（回归系数检验）。
#这里p=1.16e-9表示广告支出是影响销售收入的重要因素
#Multiple R_Squared表示R2值（决定系数），越接近一则表示拟合效果越好，这里的决定系数R2=87.82%
#表示在销售收入取值的总误差中，有87.82%是可以由销售收入和广告支出之间的线性关系解释的
#输出方差分析表，用F检验，检验线性关系检验
anova(modle)
#结果里面有F值P值，P<0.05表示拒绝H0，即销售收入与广告支出之间的线性关系显著
#对新数据进行预测
x0=data.frame(广告支出=500)
predict(model,newdata=x0)
#求出均值置信区间：
predict(model,data.frame(广告支出=500),interval="confidence",level=0.95)
#求出预测区间：（预测区间大于置信区间）
predict(model,data.frame(广告支出=500),interval="prediction",level=0.95)
#模型的诊断：正态性检验，方差齐性检验
#正态性检验
ks.test(T1$广告支出,"pnorm",mean(T1$广告支出),sd(T1$广告支出))
shapiro.test(T1$广告支出)
#或者直接绘制模型的四个图看正态QQ图
plot(model1)
#方差齐性检验,在上面的四个图中的第三个（位置尺图）就可以看方差齐性，数据点在水平线周围随机分布
library(car)
ncvTest(model)#看P值是否大于0.05,大于则说明满足方差齐性
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
