悬赏R语言作业答案

# 一、R基本 *** 作

# 1、将数据文件mydata1.txt按照以下要求整理成标准形式。

#(1)读入数据文件mydata.txt命名为insurance。

insurance<-read.table("mydata1.txt")

head(insurance)

dim(insurance)#192个数据

#(2)将insurance转换为3列的矩阵。

insurance<-matrix(insurance$V1,nrow = 64,ncol = 3)#nrow =192/3=64

insurance

#(3)将insurance转换为数据框。

insurance<-as.data.frame(insurance)

class(insurance)

#(4)将列名命名为"District", "Holders"和"Claims"。

names(insurance)<-c("District", "Holders","Claims")

insurance

#(5)随机无放回抽取50行数据。

sub<-insurance[sample(1:nrow(insurance),50),]#无放回不用设置replace

sub

#(6)将抽样数据写入result1.txt。

write.table(sub,"result1.txt",row.names = FALSE)

######################################################################

# 2、将数据文件mydata2.txt按照以下要求整理成标准形式。

#(1)读入数据文件mydata2.txt命名为iris。

iris<-read.table("mydata2.txt")

head(iris)

dim(iris)#600个数据

#(2)将iris转换为4列的矩阵。

iris<-matrix(iris$V1,nrow = 150,ncol = 4)#nrow =600/3=150

iris

#(3)将iris转换为数据框。

iris<-as.data.frame(iris)

class(iris)

#(4)将列名命名为"Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"。

names(iris)<-c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")

iris

#(5)随机无放回抽取90行数据。

sub<-iris[sample(1:nrow(iris),90),]#无放回不用设置replace

sub

#(6)将抽样数据写入result2.txt。

write.table(sub,"result2.txt",row.names = FALSE)

######################################################################

# 3.将数据文件data.csv按照以下要求进行数据预处理。

#(1)读入数据文件data.csv命名为nhanes2。

nhanes2<-read.csv("data.csv")

#(2) 载入缺失值处理所需要的包。

install.packages("lattice")

install.packages("MASS")

install.packages("nnet")

library(lattice)

library(MASS)

library(nnet)

#(3) 判断nhanes2是否存在缺失值。

sum(is.na(nhanes2))

#(4) 利用插补法处理chl变量的缺失值。

sub=which(is.na(nhanes2[,4])==TRUE)#在数据集中chl变量是第4列，返回nhanes2数据集中第4列为NA的行

dataTR<-nhanes2[-sub,]#将第4列不为NA的数存入数据集dataTR

dataTE<-nhanes2[sub,]#将第4列为NA的数存入数据集dataTE中

dataTE[,4]<-sample(dataTR[,4],length(dataTE[,4]),replace = T)#在非缺失值中简单抽样

dataTE

#(5) 将插补法处理后的数据写入result3.txt。

write.table(dataTE,"result3.txt",row.names = FALSE)

######################################################################

######################################################################

#二、函数调用

#1、测得某班学术X（身高(cm)）与Y（体重(kg)）的数据如下,试画出散点图,建立线性回归方程，并作进一步分析。

# (1) 建立数据集，并画出散点图，考察数据点的分布趋势，看是否呈直线条状分布。

x1<-c(171,175,159,155,152,158,154,164,168,166,159,164)#身高

y1<-c(57,64,41,38,35,44,41,51,57,49,47,46)#体重

#构建数据集

model <- data.frame(x1,y1)

#探索性分析-做散点图查看数据的分布情况：

plot(x1,y1)

# (2)进行回归分析,列出回归方程，画拟合线，并对结果进行解读。

# 用lm()函数构建模型

lm.reg<-lm(y1~ x1)

# 添加回归曲线查看拟合效果

abline(lm.reg)

#模型解读

summary(lm.reg)

# (3)对回归系数进行假设检验。

anova(lm.reg) # 回归模型的方差分析

summary(lm.reg) #回归系数t检验：提取模型计算结果，其中有t检验的结果

# (4)对回归模型进行诊断。

#模型检验对方程进行进一步检验，以检查回归方程是否满足模型的先验条件及模型的稳健性。

par(mfrow=c(2,2))#画布分面

plot(lm.reg)

#结果解读：

#1.左上图：残差与拟合图，理论上散点应该散乱的分布在横线两侧；

#2.右上图：正太Q-Q图，用于检验因变量的正太分布性，若服从正太分布，则散点应分布在一条直线线

#3.左下图：齐方差检验，若满足其方差，则散点在水平线周围随机分布

#4.右下图：独立性检验，即一个样本是否会影响另一个样本

#################################################################

#2、研究某抗心律失常药对电刺激狗右心室致颤阙的影响,实验测得狗静脉注射不同剂量的抗心律失常药与右心室致颤阙的数据如下,试画出散点图,建立线性回归方程，并作进一步分析。

# (1) 建立数据集，并画出散点图，考察数据点的分布趋势，看是否呈直线条状分布。

x <- c(1,3,5,7,9)

y <- c(8.03, 14.97, 19.23, 27.83, 36.23)

#构建数据集

model <- data.frame(x,y)

#探索性分析-做散点图查看数据的分布情况：

plot(model)#画散点图

# (2)进行回归分析,列出回归方程，画拟合线，并对结果进行解读。

# 用lm()函数构建模型

fm <- lm(y ~ x)#建立回归模型

fm

# 添加回归曲线查看拟合效果

abline(fm)# 添加回归曲线至散点图

#模型解读

summary(fm)

# (3)对回归系数进行假设检验。

anova(fm) # 回归模型的方差分析

summary(fm) # 提取模型计算结果，其中有t检验的结果

# (4)对回归模型进行诊断。

#模型检验对方程进行进一步检验，以检查回归方程是否满足模型的先验条件及模型的稳健性。

par(mfrow=c(2,2))#画布分面

plot(fm)

#结果解读：

#1.左上图：残差与拟合图，理论上散点应该散乱的分布在横线两侧；

#2.右上图：正太Q-Q图，用于检验因变量的正太分布性，若服从正太分布，则散点应分布在一条直线线

#3.左下图：齐方差检验，若满足其方差，则散点在水平线周围随机分布

#4.右下图：独立性检验，即一个样本是否会影响另一个样本

##################################################################

# 3、countries数据集含有69个国家和地区的出生率与死亡率。

# (1) 请使用K-均值聚类将样本点聚为3个类别。

countries=read.csv("countries.csv")

head(countries)#查看前6行

names(countries)=c("country","birth","death")#修改变量名称

var=as.character(countries$country)#将变量country转为字符型并赋值给var

for(i in 1:69) row.names(countries)[i]=var[i]#将数据集的行名命名为国家名称

km1=kmeans(countries[,-1],center=3)#用kmeans算法对countries数据集进行聚类

# (2) 输出聚类结果及各类别的中心点坐标。

km1$cluster#获取类别

km1$centers#获取中心点坐标

# (3) 绘制聚类结果将中心点以星号标识。

#画出聚为四类的类别图，标注中心点。

plot(countries[,-1],pch=c(1,2,3))

#将中心点用星号标示出来

points(km1$centers,pch=8,col="red")

#对中心点添加标注

legend(km1$centers[1,1],km1$centers[1,2],"Center_1",bty="n",xjust=0.5,cex=0.8)

legend(km1$centers[2,1],km1$centers[2,2],"Center_2",bty="n",xjust=0.5,cex=0.8)

legend(km1$centers[3,1],km1$centers[3,2],"Center_3",bty="n",xjust=0.5,cex=0.8)

# (4) 判断与中国大陆同属于一个类别的国家和地区有哪些。

cluster_CHINA=km1$cluster[which(countries$country=="CHINA")]

which(km1$cluster==cluster_CHINA)

###############################################################

###############################################################

#三、数据分析

# 1、使用arules软件包中的Groceries数据集，该数据集是某一食品杂货店一个月的真实交易数据，使用R完成以下要求：(软件包：arules；数据集：Groceries； 函数：apriori())

# (1)利用apriori()函数进行关联分析，支持度为0.01，置信度为0.5。

install.packages("arules")

library(arules)

data("Groceries")

rules0<-apriori(Groceries,parameter=list(support=0.01,confidence=0.5))

inspect(rules0[1:10])

# (2)利用sort()函数按照支持度排序。

rules.sorted_sup<-sort(rules0,by="support")

inspect(rules.sorted_sup[1:5])

# (3)捆绑销售：寻找蛋黄酱（mayonnaise）的捆绑商品。（supp=0.001,conf=0.1,minlen=2, maxlen=6）

rules1=apriori(Groceries,parameter=list(minlen=2,maxlen=6,supp=0.001,conf=0.1),appearance=list(rhs="mayonnaise",default="lhs"))

inspect(rules1)

# (4)查看销量最高的商品。

itemsets_apr=apriori(Groceries,parameter=list(supp=0.001,target="frequent itemsets"),control=list(sort=-1))

inspect(itemsets_apr[1:5])

# (5)适合捆绑销售的商品。（supp=0.001,minlen=2, maxlen=3）

itemsets_apr1=eclat(Groceries,parameter=list(supp=0.001,minlen=2,maxlen=3,target="frequent itemsets"),control=list(sort=-1))

inspect(itemsets_apr1[1:5])

# (6)关联规则的可视化（support=0.001,con=0.5）

install.packages("arulesViz")

library(arulesViz)

rules5=apriori(Groceries,parameter=list(support=0.002,con=0.5))

rules5

plot(rules5)

#######################################################################

# 2、根据breast-cancer-wisconsin.csv威斯康星州乳腺癌数据集，通过对数据的分析，提取出关键特征来判断乳腺癌患病情况。(软件包：rpart；函数：rpart()。)

# (1)属性名依次设置为"编号","肿块厚度","肿块大小","肿块形状","边缘黏附","单个表皮细胞大小","细胞核大小","染色质","细胞核常规","有丝分裂","类别"),并将类别为2的设为"良性",为4的设为"恶性"。

install.packages("rpart")

library(rpart)

install.packages("rpart.plot")

library(rpart.plot)

#############加载数据

breast.cancer<-read.csv('breast-cancer-wisconsin.csv',header=F)

head(breast.cancer)

#数据整理

names(breast.cancer)=c("编号","肿块厚度","肿块大小","肿块形状","边缘黏附","单个表皮细胞大小","细胞核大小","染色质","细胞核常规","有丝分裂","类别")

breast.cancer$类别[breast.cancer$类别==2]="良性"

breast.cancer$类别[breast.cancer$类别==4]="恶性"

head(breast.cancer)

# (2)抽取训练数据集为原数据的70%，测试数据集取30%。

#数据预处理(分层抽样,划分训练集和测试集)

#分别计算良性和恶性组中应抽取测试集样本数，记为a，b

a=round(0.3*sum(breast.cancer$类别=="良性"))

b=round(0.3*sum(breast.cancer$类别=="恶性"))

ab #输出a,b值

install.packages("sampling")

library(sampling)

#使用strata函数对数据集中的“分组油耗”变量进行分层抽样

sub=strata(breast.cancer,stratanames="类别",size=c(b,a),method="srswor")

sub #所抽出的所有测试集样本信息

#生成训练集train1和测试集test1

train1=breast.cancer[-sub$ID_unit,]

test1=breast.cancer[sub$ID_unit,]

nrow(train1)nrow(test1) #显示训练集和测试集的行数，检查两者比例是否为7：3

# (3) minsplit=5，建立决策树。

#CART建立分类树

formula_cla=类别~肿块厚度+肿块大小+肿块形状+边缘黏附+单个表皮细胞大小+细胞核大小+染色质+细胞核常规+有丝分裂

cla1=rpart(formula_cla,train1,method="class",minsplit=5)#

cla1

# (4)选择cp=0.05来剪枝。

######修改cp的值

cla2=rpart(formula_cla,train1,method="class",minsplit=5,cp=0.05)

cla2

# (5)画出type为2和4的树图。

rpart.plot(cla1,type=2)#修改type

rpart.plot(cla1,type=4)

# (6)测试数据进行预测，并输出混淆矩阵，给出模型准确率为。

#预测

pre1=predict(cla1,test1,type="class")

pre1

table(test1$类别,pre1)#获取混淆矩阵

#计算样本错误率

error1<-sum(as.numeric(pre1!=test1$类别))/nrow(test1)

error1

###################################################################

# 3、美国科罗拉多州某加油站连续 57 天的OVERSHORTS序列“OVERSHORTS.csv”

# (1) 判断该序列的平稳性与纯随机性。

# （时序图检验、白噪声检验）

install.packages("fUnitRoots")

install.packages("TSA")

install.packages("forecast")

install.packages("zoo")

library(fUnitRoots)

library(TSA)

library(forecast)

library(zoo)

#读取数据

c<-read.csv("OVERSHORTS.csv")

#转换为时间序列

overshort<-ts(c$overshort,start = 1)

#平稳性，纯随机（白噪声检验）

## 绘制序列的时间序列图

plot.ts(overshort, xlab = "time", ylab = "prop")

##对序列做单位根检验

unitrootTest(overshort)

##对序列做白噪声检验

Box.test(overshort, lag = 1, type = "Ljung-Box")

# (2) 如果序列平稳且非白噪声，选择适当模型拟合该序列的发展。(10分)

# （模型的识别、参数估计（模型显著性、模型参数的显著性））

#模型识别

##观察自相关，偏自相关图，模型定阶

par(mfrow=c(1,2))

acf(overshort)###衰减到零是突然的，所以自相关系数1阶截尾

pacf(overshort)### 衰减到零不是突然的，所以偏相关系数托尾

# 推荐模型为 MA(1)

##或者对序列进行模型识别，自动定阶

auto.arima(overshort)# 推荐模型为 MA(1)

#参数估计

###模型检验

x.fit<-arima(overshort,order=c(0,0,1),method="ML")

x.fit

##对残差x.fit$residual进行白噪声检验

for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

##P>0.05,接受原假设，即残差为白噪声，所以拟合模型显著有效

####参数检验

###模型参数的显著性检验

t1<--0.8477/0.1206

pt(t1,df=56,lower.tail=T) ###p<0.05参数显著非零

t0<--4.7942/1.0253

pt(t0,df=56,lower.tail=T) ###p<0.05参数显著非零

# (3) 利用拟合模型，预测该加油站未来5天的OVERSHORTS。(10分)

# （模型预测、绘制预测图）

####模型预测

c<-read.csv("OVERSHORTS.csv")

x<-ts(c$overshort,start=1)

x.fit<-arima(x,order=c(0,0,1))

x.fit

x.fore<-forecast(x.fit,h=5)#预测

x.fore

plot(x.fore)

##############################################################

#4、使用是survival软件包中的“pbc”数据集，该数据集记录的是肝硬化数据， 使用R完成一下要求：(软件包：survival；数据集：pbc； 函数：Surv()、survfit()、survdiff()、coxph()、cox.zph(), 将答案保存在“姓名.doc”文件中。)

# (1)生成生存分析对象，拟合生存曲线模型。

install.packages("survival") #安装survival包

library(survival) #加载survival包

#使用survival包自带的“pbc”数据集为例（418*20）

data("pbc")

str(pbc)

head(pbc)

#生成生存分析对象

Sur_Obj<-Surv(pbc$time,pbc$status)

Sur_Obj

#拟合曲线模型

model<-survfit(Sur_Obj~1)

summary(model)

# (2)两种方法绘制生存曲线。

plot(model,ylab = "生存率",xlab="天")

#用survminer进行漂亮的展示

install.packages("survminer")

library(survminer)

ggsurvplot(model, data = pbc)

# (3)进行单因素比较分析，并进行结果解释。

#survdiff(formula)函数进行log-rank检验。

survdiff(Sur_Obj~pbc$trt) #trt是分组条件

# (4)考虑年龄，性别以及trt是否会影响肝硬化的生存时间,进行多因素分析Cox模型的建立，并进行结果解释。

coxmodel<-coxph(Sur_Obj~pbc$age+pbc$sex+pbc$bili)

coxmodel

# (5)模型诊断——PH检验。

zphmodel<-cox.zph(coxmodel)

zphmodel

##############################################################

# 5、life.csv为50位急性淋巴细胞白血病病人的数据，包括：入院治疗时取得外辕血中细胞数X1，淋巴结浸润等级X2，出院后有无巩固治疗X3(1表示有巩固治疗，0表示无巩固治疗)；随访后，变量Y=0表示生存期在1年以内，Y=1表示生存时间在1年以上，使用R完成一下要求：(函数：glm()，predict()。)

# (1)建立全变量logistic回归，对模型结果进行解释。

life<-read.csv("life.csv")

#建立全变量logistic回归

glm.sol<-glm(Y~X1+X2+X3, family=binomial, data=life)

#回归模型解读

summary(glm.sol)

# (2)预测当X1=5,X2=2,X3=0时，y的概率是多少？

pre<-predict(glm.sol, data.frame(X1=5,X2=2,X3=0))

p<-exp(pre)/(1+exp(pre))

p

# (3)预测当X1=5,X2=2,X3=1时，y的概率是多少？(6分)

pre<-predict(glm.sol, data.frame(X1=5,X2=2,X3=1))

p<-exp(pre)/(1+exp(pre))

p

# (4)对回归模型参数进行检验，用step()函数做变量筛选。

step(glm.sol)

glm.new<-glm(Y~X2+X3, family=binomial, data=life)

summary(glm.new)

# (5)对筛选后的变量进行建模，预测。

pre<-predict(glm.new, data.frame(X2=2,X3=0))

p<-exp(pre)/(1+exp(pre))

p

pre<-predict(glm.new, data.frame(X2=2,X3=1))

p<-exp(pre)/(1+exp(pre))

p

在日常学习或工作中经常会使用线性回归模型对某一事物进行预测，例如预测房价、身高、GDP、学生成绩等，发现这些被预测的变量都属于连续型变量。然而有些情况下，被预测变量可能是二元变量，即成功或失败、流失或不流失、涨或跌等，对于这类问题，线性回归将束手无策。这个时候就需要另一种回归方法进行预测，即Logistic回归。

在实际应用中，Logistic模型主要有三大用途：

1）寻找危险因素，找到某些影响因变量的"坏因素"，一般可以通过优势比发现危险因素；

2）用于预测，可以预测某种情况发生的概率或可能性大小；

3）用于判别，判断某个新样本所属的类别。

Logistic模型实际上是一种回归模型，但这种模型又与普通的线性回归模型又有一定的区别：

1）Logistic回归模型的因变量为二分类变量；

2）该模型的因变量和自变量之间不存在线性关系；

3）一般线性回归模型中需要假设独立同分布、方差齐性等，而Logistic回归模型不需要；

4）Logistic回归没有关于自变量分布的假设条件，可以是连续变量、离散变量和虚拟变量；

5）由于因变量和自变量之间不存在线性关系，所以参数(偏回归系数)使用最大似然估计法计算。

logistic回归模型概述

广义线性回归是探索“响应变量的期望”与“自变量”的关系，以实现对非线性关系的某种拟合。这里面涉及到一个“连接函数”和一个“误差函数”，“响应变量的期望”经过连接函数作用后，与“自变量”存在线性关系。选取不同的“连接函数”与“误差函数”可以构造不同的广义回归模型。当误差函数取“二项分布”而连接函数取“logit函数”时，就是常见的“logistic回归模型”，在0-1响应的问题中得到了大量的应用。

Logistic回归主要通过构造一个重要的指标：发生比来判定因变量的类别。在这里我们引入概率的概念，把事件发生定义为Y=1，事件未发生定义为Y=0，那么事件发生的概率为p,事件未发生的概率为1-p，把p看成x的线性函数；

回归中，最常用的估计是最小二乘估计，因为使得p在[0,1]之间变换，最小二乘估计不太合适，有木有一种估计法能让p在趋近与0和1的时候变换缓慢一些（不敏感），这种变换是我们想要的，于是引入Logit变换,对p/(1-p)也就是发生与不发生的比值取对数，也称对数差异比。经过变换后，p对x就不是线性关系了。

使用R做回归分析整体上是比较常规的一类数据分析内容，下面我们具体的了解用R语言做回归分析的过程。

首先，我们先构造一个分析的数据集

x<-data.frame(y=c(102,115,124,135,148,156,162,176,183,195),

var1=runif(10,min=1,max=50),

var2=runif(10,min=100,max=200),

var3=c(235,321,412,511,654,745,821,932,1020,1123))

接下来，我们进行简单的一元回归分析，选择y作为因变量，var1作为自变量。

---