最终出图如下
这里自动做统计检验的函数是 stat_compare_means()
读入数据
作图
这个函数来自于ggpubr这个包，只需要指定根据那一列来分组就可以了
默认的是Wilcoxon Rank Sum and Signed Rank Tests，如果要用t检验指定method参数
如果想把P值改成星号，直接加label=“p.signif”参数
这里如果不显著会在图上显示ns，如果不想要ns，可以加 hide.ns = TRUE 参数
星号的位置可以手动指定，用 label.y = c(26,31) 参数
使用到的是 ggsignif 这个包
小明的数据分析笔记本
在r中看函数源代码：在R中，代码可以分为如下几个级别：
首先，是你输入了函数对象名称，你可以直接看到代码的，如要获得函数对象fivenum的代码，就只需要在Console中键入函数对象名称fivenum就可以得到如下结果：
function (x, na.rm = TRUE)
{
xna <- is.na(x)
if (na.rm)
x <- x[!xna]
else if (any(xna))
return(rep.int(NA, 5))
x <- sort(x)
n <- length(x)
if (n == 0)
rep.int(NA, 5)
else {
n4 <- floor((n + 3)/2)/2
d <- c(1, n4, (n + 1)/2, n + 1 - n4, n)
0.5 * (x[floor(d)] + x[ceiling(d)])
}
}
<environment: namespace:stats>
从上面的例子可以看出，这类函数对象的代码是最容易看到的，也是我们学习的最好的材料了，而R中最大多数的函数对象是以这种方式出现的。
其次，我们在输入mean这类函数名次的时候，会出现如下结果：
function (x, ...)
UseMethod("mean")
<environment: namespace:base>
这表示函数作者把函数“封”起来了。这个时候我们可以先试一试methods(mean)，利用methods函数看看mean这个函数都有哪些类型的，我们得到的结果如下：
[1] mean.data.frame mean.Date mean.defaultmean.difftime mean.POSIXct mean.POSIXlt
其实对此可以有一个简单的理解，虽然不够精确。因为在R中，mean函数可以求得属于不同类型对象的平均值，而不同类型对象平均值的求法还是有一些小小差 异的，比如说求一个向量的平均值和求一个数据框的平均值就有所差异，就要编写多个mean函数，然后“封”起来，以一个统一的mean出现，方便我们使 用。这正好也反映了R有一种类似泛型编程语言的性质。
既然我们已经知道mean中还有这么多种类，我们可以输入mean.default试一试就可以得到：
function (x, trim = 0, na.rm = FALSE, ...)
{
if (!is.numeric(x) &&!is.complex(x) &&!is.logical(x)) {
warning("argument is not numeric or logical: returning NA")
return(as.numeric(NA))
}
if (na.rm)
x <- x[!is.na(x)]
trim <- trim[1]
n <- length(x)
if (trim >0 &&n >0) {
if (is.complex(x))
stop("trimmed means are not defined for complex data")
if (trim >= 0.5)
return(stats::median(x, na.rm = FALSE))
lo <- floor(n * trim) + 1
hi <- n + 1 - lo
x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]
n <- hi - lo + 1
}
.Internal(mean(x))
}
<environment: namespace:base>
同样就可以得到mean.data.frame、mean.Date、mean.difftime、mean.POSIXct、mean.POSIXlt 的具体内容了。值得注意的是，在R中，出现有多个同样近似功能的函数封装为一个函数的时候（这时候在函数中多半会出现类似UseMethod函数使用的情 况），我们不妨先输入mean.default试一试。这种形式的函数在R中一般作为默认的函数表示。
第三，这是一种特殊的情况，有人认为应该和第二种是一类，但是我还是要提出来单独归类。在这种情况也和第二种的原因有些类似，但并不是完全一致。
也许我们大家都很熟悉plot函数了吧，输入函数名plot的时候，我们会得到如下结果：
function (x, y, ...)
{
if (is.null(attr(x, "class")) &&is.function(x)) {
nms <- names(list(...))
if (missing(y))
y <- {
if (!"from" %in% nms)
0
else if (!"to" %in% nms)
1
else if (!"xlim" %in% nms)
NULL
}
if ("ylab" %in% nms)
plot.function(x, y, ...)
else plot.function(x, y, ylab = paste(deparse(substitute(x)),
"(x)"), ...)
}
else UseMethod("plot")
}
<environment: namespace:graphics>
请注意plot函数中也出现了UseMethod这个函数，但是和mean不同的是，前面有相当多的语句用于处理其他一些事情。这个时候，我们也使用methods(plot)来看看，得到如下结果：
plot.acf* plot.data.frame*plot.Date* plot.decomposed.ts* plot.default
plot.dendrogram*plot.densityplot.ecdf plot.factor*plot.formula*
plot.hclust*plot.histogram* plot.HoltWinters* plot.isoreg*plot.lm
plot.medpolish* plot.mlmplot.POSIXct* plot.POSIXlt* plot.ppr*
plot.prcomp*plot.princomp* plot.profile.nls* plot.spec plot.spec.coherency
plot.spec.phase plot.stepfunplot.stl* plot.table* plot.ts
plot.tskernel* plot.TukeyHSD
不看不知道，一看吓一跳，还以为我们输入plot的输出就是函数本身，结果也许不是如此。可能有人已经理解了，其实最后的UseMethod函数实在默认的调用plot.default函数，赶快再看看plot.default函数吧，发现它再调用plot.xy函数，再看看plot.xy函数，再plot.xy函数中调用了一个.Internal(plot.xy(xy, type, pch, lty, col, bg, cex, lwd, ...))函数，也许这就是真正起作用的函数了吧。思路基本上就是如此了，是否这个时候您可以获得一些阅读查找R函数内容的乐趣。
除了直接输入FUN.default形式外，还可以使用getS3method(FUN,"default")来获得代码。这样就解决了绝大多数函数代码查看的工作了。
在第二种情况种，我们说了一般可以通过FUN.default获得想要的结果。但是只有称为generic的函数才有这种“特权”。而lm等则没有，不过我们也可以尝试使用methods(lm)来看看结果如何，发现：
[1] lm.fit lm.fit.null lm.influence lm.wfit lm.wfit.null
Warning message:
function 'lm' appears not to be generic in: methods(lm)
出现了警告信息，表示说lm不是泛型函数，但是还是给出了结果lm.fit等，大致上可以看成是和lm相关的系列函数吧。这样子就出现了有趣的局面，比如说既有plot.ts，也有ts.plot。
依照第三种情况，我们发现竟然有的函数用星号标识了的，比如plot.stl*等，当我们输入plot.stl，甚至是plot.stl*的时候都会给出 要么找不到这个对象，要么干脆是代码错误的信息。原来凡是用了*标识的函数，都是隐藏起来的函数，估计是怕被人看见（其实这是玩笑话）！我们要看这些函数 的代码，我们该怎么办呢？其实也很容易，使用功能强大的getAnywhere(FUN)，看看这个函数的名称，就可以猜想到它的功能估计是很强大的， Anywhere的内容都可以找到！getAnywhere(plot.stl)的结果如下：
A single object matching 'plot.stl' was found
It was found in the following places
registered S3 method for plot from namespace stats
namespace:stats
with value
function (x, labels = colnames(X), set.pars = list(mar = c(0,
6, 0, 6), oma = c(6, 0, 4, 0), tck = -0.01, mfrow = c(nplot,
1)), main = NULL, range.bars = TRUE, ..., col.range = "light gray")
{
sers <- x$time.series
ncomp <- ncol(sers)
data <- drop(sers %*% rep(1, ncomp))
X <- cbind(data, sers)
colnames(X) <- c("data", colnames(sers))
nplot <- ncomp + 1
if (range.bars)
mx <- min(apply(rx <- apply(X, 2, range), 2, diff))
if (length(set.pars)) {
oldpar <- do.call("par", as.list(names(set.pars)))
on.exit(par(oldpar))
do.call("par", set.pars)
}
for (i in 1:nplot) {
plot(X[, i], type = if (i <nplot)
"l"
else "h", xlab = "", ylab = "", axes = FALSE, ...)
if (range.bars) {
dx <- 1/64 * diff(ux <- par("usr")[1:2])
y <- mean(rx[, i])
rect(ux[2] - dx, y + mx/2, ux[2] - 0.4 * dx, y -
mx/2, col = col.range, xpd = TRUE)
}
if (i == 1 &&!is.null(main))
title(main, line = 2, outer = par("oma")[3] >0)
if (i == nplot)
abline(h = 0)
box()
right <- i%%2 == 0
axis(2, labels = !right)
axis(4, labels = right)
axis(1, labels = i == nplot)
mtext(labels[i], side = 2, 3)
}
mtext("time", side = 1, line = 3)
invisible()
}
<environment: namespace:stats>
注意到前面有一段解释型的语言，描述了我们要找的这个函数放在了什么地方等等。其实对任意我们可以在R中使用的函数，都可以先试一试getAnywhere，看看都有些什么内容。算是一个比较“霸道”的函数。
在上面plot.xy函数中，我们还可以看到.Internal这个函数，类似的也许还可以看到.Primitive、.External、.Call等函数这就和R系统内部工作方式和与外部接口的定义有关了，如果对这些函数有兴趣的话，就要学习组成R系统的源代码了。
最后，如果真的想阅读组成R系统本身的源代码，在各个CRAN中均有下载。你可以得到组成R系统所需要的材料。其中很多C语言（还有就是F）的源代码，均 是精心挑选过的算法，哪怕就是想学从头到尾编写具体的算法，也是学习的好材料。同时，你可以看到R系统内部是如何构成的，理解了这些对于高效使用R有至关 重要的作用。这个范畴的材料就要着重看一看R-Lang和R-inits了。
至此，R中阅读代码的内容就依照我的理解介绍了一下。随后将有一些R代码示例的分析注解、语言本身、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
