需要解释gc()函数，可以查看内存使用情况。同样，在清除了大的对象之后，使用gc()以释放内存使用空间。
李航在”参考网址2“中提到，对于大矩阵的操作，尽量避免使用cbind和rbind之类，因为这会让内存不停地分配空间。“对于长度增加的矩阵，尽量先定义一个大矩阵，然后逐步增加”和“注意清除中间对象”。
使用bigmemory家族：bigmemory, biganalytics, synchronicity, bigtabulate and bigalgebra， 同时还有
biglm。
bigmemory package的使用：
1. 建立big.memory对象
bigmemory采用C++的数据格式来“模仿”R中的matrix。
编写大数据格式文件时候，可以先建立filebacked.big.matrix
big.matrix(nrow, ncol, type = options()$bigmemory.default.type, init = NULL, dimnames = NULL, separated = FALSE, backingfile = NULL, backingpath = NULL, descriptorfile = NULL, shared = TRUE)
filebacked.big.matrix(nrow, ncol, type = options()$bigmemory.default.type, init = NULL, dimnames = NULL, separated = FALSE, backingfile = NULL, backingpath = NULL, descriptorfile = NULL)
as.big.matrix(x, type = NULL, separated = FALSE, backingfile = NULL, backingpath = NULL, descriptorfile = NULL, shared=TRUE)
使用注意：
big.matrix采用两种方式储存数据：一种是big.matrix默认的方式，如果内存空间比较大，可以尝试使用；另外一种是filebacked.big.matrix，这种储存方法可能会备份文件（file-backings），而且需要descriptor file；
“init”指矩阵的初始化数值，如果设定，会事先将设定的数值填充到矩阵中；如果不设置，将处理为NA
"type"是指在big.matrix中atomic element的储存格式，默认是“double”(8 byte)，可以改为“integer”(4 byte), "short"(2 byte) or "char"(1 byte)。注意：这个包不支持字符串的储存，type = "char"是指ASCII码字母。
在big.matrix非常大的时候，避免使用rownames和colnames(并且bigmemory禁止用名称访问元素)，因为这种做法非常占用内存。如果一定要改变，使用options(bigmemory.allow.dimnames=TRUE)，之后colnames, rownames设置。
直接在命令提示符后输入x（x是一个big matrix），将返回x的描述，不会出现所有x中所有内容。因此，注意x[ , ](打印出矩阵全部内容)；
如果big.matrix有很多列，那么应该将其转置后储存；（不推荐）或者将参数“separated”设置为TRUE，这样就将每一列分开储存。否则，将用R的传统方式（column major的方式）储存数据。
如果建立一个filebacked.big.matrix，那么需要指定backingfile的名称和路径+descriptorfile。可能多个big.matrix对象对应唯一一个descriptorfile，即如果descriptorfile改变，所以对应的big.matrix随之改变；同样，decriptorfile随着big.matrix的改变而改变；如果想维持一种改变，需要重新建立一个filebacked.big.matrix。attach.big.matrix(descriptorfile or describe(big.matrix))函数用于将一个descriptorfile赋值给一个big.matrix。这个函数很好用，因为每次在创建一个filebacked.big.matrix后，保存R并退出后，先前创建的矩阵会消失，需要再attach.big.matrix以下
2. 对big.matrix的列的特定元素进行条件筛选
对内存没有限制；而且比传统的which更加灵活（赞！）
mwhich(x, cols, vals, comps, op = 'AND')
x既可以是big.matrix，也可以是传统的R对象；
cols：行数
vals：cutoff，可以设定两个比如c(1, 2)
comps：'eq'(==), 'neq'(!=), 'le'(<), 'lt'(<=), 'ge'(>) and 'gt'(>=)
op：“AND”或者是“OR”
可以直接比较NA，Inf和-Inf
3.bigmemory中其他函数
nrow, ncol, dim, dimnames, tail, head, typeof继承base包
big.matrix, is.big.matrix, as.big.matrix, attach.big.matrix, describe, read.big.matrix, write.big.matrix, sub.big.matrix, is.sub.big.matrix为特有的big.matrix文件操作；filebacked.big.matrix, is.filebacked（判断big.matrix是否硬盘备份） , flush(将filebacked的文件刷新到硬盘备份上)是filebacked的big.matrix的操作。
mwhich增强base包中的which， morder增强order，mpermute（对matrix中的一列按照特定序列操作，但是会改变原来对象，这是为了避免内存溢出）
big.matrix对象的copy使用deepcopy(x, cols = NULL, rows = NULL, y = NULL, type = NULL, separated = NULL, backingfile = NULL, backingpath = NULL, descriptorfile = NULL, shared=TRUE)
biganalytics package的使用
biganalytics主要是一些base基本函数的扩展，主要有max, min, prod, sum, range, colmin, colmax, colsum, colprod, colmean, colsd, colvar, summary, apply（只能用于行或者列，不能用行列同时用）等
比较有特色的是bigkmeans的聚类
剩下的biglm.big.matrix和bigglm.big.matrix可以参考Lumley's biglm package。
bigtabulate package的使用
并行计算限制的突破：
使用doMC家族：doMC, doSNOW, doMPI, doRedis, doSMP和foreach packages.
foreach package的使用
foreach(..., .combine, .init, .final=NULL, .inorder=TRUE, .multicombine=FALSE, .maxcombine=if (.multicombine) 100 else 2, .errorhandling=c('stop', 'remove', 'pass'), .packages=NULL, .export=NULL, .noexport=NULL, .verbose=FALSE)
foreach的特点是可以进行并行运算，如在NetWorkSpace和snow？
%do%严格按照顺序执行任务（所以，也就非并行计算），%dopar%并行执行任务
...：指定循环的次数；
.combine：运算之后结果的显示方式，default是list，“c”返回vector， cbind和rbind返回矩阵，"+"和"*"可以返回rbind之后的“+”或者“*”
.init：.combine函数的第一个变量
.final：返回最后结果
.inorder：TRUE则返回和原始输入相同顺序的结果（对结果的顺序要求严格的时候），FALSE返回没有顺序的结果（可以提高运算效率）。这个参数适合于设定对结果顺序没有需求的情况。
.muticombine：设定.combine函数的传递参数，default是FALSE表示其参数是2，TRUE可以设定多个参数
.maxcombine：设定.combine的最大参数
.errorhandling：如果循环中出现错误，对错误的处理方法
.packages：指定在%dopar%运算过程中依赖的package（%do%会忽略这个选项）。
getDoParWorkers( ) ：查看注册了多少个核，配合doMC package中的registerDoMC( )使用
getDoParRegistered( ) ：查看doPar是否注册；如果没有注册返回FALSE
getDoParName( ) ：查看已经注册的doPar的名字
getDoParVersion( )：查看已经注册的doPar的version
===================================================
# foreach的循环次数可以指定多个变量，但是只用其中最少？的
>foreach(a = 1:10, b = rep(10, 3)) %do% (a*b)
[[1]]
[1] 10
[[2]]
[1] 20
[[3]]
[1] 30
# foreach中.combine的“+”或者“*”是cbind之后的操作；这也就是说"expression"返回一个向量，会对向量+或者*
>foreach(i = 1:4, .combine = "+") %do% 2
[1] 8
>foreach(i = 1:4, .combine = "rbind") %do% rep(2, 5)
[,1] [,2] [,3] [,4] [,5]
result.122222
result.222222
result.322222
result.422222
>foreach(i = 1:4, .combine = "+") %do% rep(2, 5)
[1] 8 8 8 8 8
>foreach(i = 1:4, .combine = "*") %do% rep(2, 5)
[1] 16 16 16 16 16
=============================================
iterators package的使用
iterators是为了给foreach提供循环变量，每次定义一个iterator，它都内定了“循环次数”和“每次循环返回的值”，因此非常适合结合foreach的使用。
iter(obj, ...)：可以接受iter, vector, matrix, data.frame, function。
nextElem(obj, ...)：接受iter对象，显示对象数值。
以matrix为例，
iter(obj, by=c('column', 'cell', 'row'), chunksize=1L, checkFunc=function(...) TRUE, recycle=FALSE, ...)
by：按照什么顺序循环；matrix和data.frame都默认是“row”，“cell”是按列依次输出（所以对于“cell”，chunksize只能指定为默认值，即1）
chunksize：每次执行函数nextElem后，按照by的设定返回结果的长度。如果返回结构不够，将取剩余的全部。
checkFunc=function(...) TRUE：执行函数checkFun，如果返回TRUE，则返回；否则，跳过。
recycle：设定在nextElem循环到底（“错误: StopIteration”）是否要循环处理，即从头再来一遍。
以function为例
iter(function()rnorm(1))，使用nextElem可以无限重复；但是iter(rnorm(1))，只能来一下。
更有意思的是对象如果是iter，即test1 <- iter(obj)test2 <- iter(test1)，那么这两个对象是连在一起的，同时变化。
==============================================
>a
[,1] [,2] [,3] [,4] [,5]
[1,]159 13 17
[2,]26 10 14 18
[3,]37 11 15 19
[4,]48 12 16 20
>i2 <- iter(a, by = "row", chunksize=3)
>nextElem(i2)
[,1] [,2] [,3] [,4] [,5]
[1,]159 13 17
[2,]26 10 14 18
[3,]37 11 15 19
>nextElem(i2) #第二次iterate之后，只剩下1行，全部返回
[,1] [,2] [,3] [,4] [,5]
[1,]48 12 16 20
>i2 <- iter(a, by = "column", checkFunc=function(x) sum(x) >50)
>nextElem(i2)
[,1]
[1,] 13
[2,] 14
[3,] 15
[4,] 16
>nextElem(i2)
[,1]
[1,] 17
[2,] 18
[3,] 19
[4,] 20
>nextElem(i2)
错误: StopIteration
>colSums(a)
[1] 10 26 42 58 74
>testFun <- function(x){return(x+2)}
>i2 <- iter(function()testFun(1))
>nextElem(i2)
[1] 3
>nextElem(i2)
[1] 3
>nextElem(i2)
[1] 3
>i2 <- iter(testFun(1))
>nextElem(i2)
[1] 3
>nextElem(i2)
错误: StopIteration
>i2 <- iter(testFun(1))
>i3 <- iter(i2)
>nextElem(i3)
[1] 3
>nextElem(i2)
错误: StopIteration
============================================
iterators package中包括
irnorm(..., count)；irunif(..., count)；irbinom(..., count)；irnbinom(..., count)；irpois(..., count)中内部生成iterator的工具，分别表示从normal，uniform，binomial，negativity binomial和Poisson分布中随机选取N个元素，进行count次。其中，negative binomial分布：其概率积累函数(probability mass function)为掷骰子，每次骰子为3点的概率为p，在第r+k次恰好出现r次的概率。
icount(count)可以生成1:conunt的iterator；如果count不指定，将从无休止生成1:Inf
icountn(vn)比较好玩，vn是指一个数值向量（如果是小数，则向后一个数取整，比如2.3 -->3）。循环次数为prod(vn)，每次返回的向量中每个元素都从1开始，不超过设定 vn，变化速率从左向右依次递增。
idiv(n, ..., chunks, chunkSize)返回截取从1:n的片段长度，“chunks”和“chunkSize”不能同时指定，“chunks”为分多少片段（长度从大到小），“chunkSize”为分段的最大长度（长度由大到小）
iapply(X, MARGIN)：与apply很像，MARGIN中1是row，2是column
isplit(x, f, drop=FALSE, ...)：按照指定的f划分矩阵
=============================================
>i2 <- icountn(c(3.4, 1.2))
>nextElem(i2)
[1] 1 1
>nextElem(i2)
[1] 2 1
>nextElem(i2)
[1] 3 1
>nextElem(i2)
[1] 4 1
>nextElem(i2)
[1] 1 2
>nextElem(i2)
[1] 2 2
>nextElem(i2)
[1] 3 2
>nextElem(i2)
[1] 4 2
>nextElem(i2)
错误: StopIteration
# 一、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=64insurance#(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),]#无放回不用设置replacesub#(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=150iris#(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),]#无放回不用设置replacesub#(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的数存入数据集dataTRdataTE<-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转为字符型并赋值给varfor(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))rules5plot(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，ba=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和测试集test1train1=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)#修改typerpart.plot(cla1,type=4) # (6)测试数据进行预测，并输出混淆矩阵，给出模型准确率为。#预测pre1=predict(cla1,test1,type="class")pre1table(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.1206pt(t1,df=56,lower.tail=T) ###p<0.05参数显著非零t0<--4.7942/1.0253pt(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.fitx.fore<-forecast(x.fit,h=5)#预测x.foreplot(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))ppre<-predict(glm.new, data.frame(X2=2,X3=1))p<-exp(pre)/(1+exp(pre))p
