r语言lattice脚本编写

⑴ 加权后的数据怎么用r转换出来

地理加权回归（GWR）在R里面怎么实现？
121 人关注0 条评论

写回答
查看全部 5 个回答

写回答

叶山Shan Ye
GIS/地质/人文地理/可持续发展
A2A 谢邀，

我和我认识的一些人，刚开始用哪中则R做空间分析的时候，也遇到过这个问题。R这种开源的东西，优点是各种包很丰富，缺点是有些包的说明写得很乱，地理加李棚权回归（GWR）的R包其实功能很培高强大，但大部分说明都不大靠谱。

GWR在R里面可以用好几个不同的包来实现，其中步骤最简单的是spgwr。思路就两步：建立窗口、用窗口扫全局。这其实就是GWR本质上的两步。比如我要在全美国范围内统计某两个（或多个）变量之间的回归关系，我可以做一个全局回归（global regression），但因为这些变量在空间分布上或许会有异质性（heterogeneity），表现在统计结果上就是空间不稳定性（nonstationarity），因此只看全局的统计，可能看不出什么结果来。举个不完全恰当但是很容易领会精神的例子，你比如说，我要分析亚洲范围内，经济发展程度与牛肉销量之间的关系，经济越发达的地方，人们就越吃得起牛肉。可是等我统计到印度的时候，坏了，印度大部分人不吃牛肉，这不是经济状况导致的，这一下就影响了全局统计的参考价值，那怎么办呢？我们可以建立一个窗口（正规说法是带宽窗口，bandwidth window），每次只统计窗口范围内的经济与牛肉销量的关系，然后用这个窗口去扫过全局的范围。等统计到印度的时候，印度内部的各地和印度自己比，吃牛肉的人的比例就不会突然减少，这样就能减少这种空间不稳定性对全局统计的影响。

所以，第一步就是要建立这样一个‘窗口’。当然了，首先要安装包，我们要用到的R包有：

library(spgwr)
library(rgdal)
library(sf)
library(spData)
library(sp)
library(lattice)
library(ggplot2)
library(ggthemes)
其中，spgwr是做GWR的包，rgdal是用来读取矢量要素的，sf，sp和spData都是用来处理矢量数据的，别的基本都是画图用。

以下默认你会R和GWR的基本操作。并且，以下只展现方法，不要纠结我的数据和结果，我随便找的数据，这个数据本身没有什么意义，所以做出的统计看起来很‘壮观’。

我们先导入数据。这里我用的是美国本土48州各个县（county，也有翻译成郡的）的人口普查数据和农业数据，来源是ESRI Online数据库。为啥用这个数据呢？因为...我电脑里面就存了这么个可以用来做GWR的数据...

我们用rgdal读取数据，然后把它画出来看看

require(rgdal)
usa_agri <- readOGR(dsn = "~/Documents/Spatial", layer = "usa_counties")
plot(usa_agri)
会得到这个东西：

readOGR里面，dsn后面加储存shp的路径（加到文件夹为止），layer后面写shp的文件名（不加.shp）。不喜欢rgdal的同学可以不用，用maptools或者spData等别的处理shp的R包代替。不过如果用maptools，要注意处理一下参考系。

我们看一下这个shp里面的列联表都有什么：

可见，shp里面有3108个县的数据，数据有61种。然后再看data下面有什么：

总之就是各种人口普查的数据，后面截不完图，还有经济、房地产和农业之类的数据。那我们就随便选两个来当变量。我就随便挑了，因变量选AVESIZE12，即2012年各个县农场的平均占地面积。自变量选POP_SQMI，也就是人口密度（每平方英里的人口）。

现在正式建立窗口，调用的是spgwr里面的gwr.sel函数：

bw <- gwr.sel( AVE_SIZE12 ~ POP_SQMI, data = usa_agri, gweight = gwr.Gauss,
verbose = FALSE, method = "cv")
其中~前后分别是因变量和自变量。GWR里因变量只能有1个，但自变量可以选多个，如果需要多个自变量的话，就在代码POP_SQMI之后用+号连接就行。gweight是你的空间加权的函数（随空间距离增大而不断衰减的函数，衰减率由下面要提到的带宽控制），这里用的是比较常用的高斯函数，其余的还有gwr.bisquare等函数可以调用。verbose决定是否汇报制定窗口的过程。method是决定构建带宽窗口模型的方法，这里用的cv指的是cross validation，即交叉验证法，也是最常用的方法，简单说就是把数据分成不同的组，分别用不同的方法来做回归计算，计算完了之后记录下结果，然后打乱重新分组，再回归计算，再看结果，周而复始，最后看哪种计算方法的结果最靠谱，这种方法就是最优解。还有一种很常见的选择最佳拟合模型的方法是AIC optimisation法，把method后面的cv改成aic就可以用。具体AIC optimisation是什么：AIC（赤池信息准则）_网络。总之，空间加权函数和带宽窗口构建方法的选择是GWR里面十分重要的步骤。

以上便是固定带宽窗口的示意图。比如我在对佐治亚做GWR，这一轮的regression target是红色的这个县，根据做出来的窗口，圆圈以内的县都要被算为红色县的邻县，其权重根据高斯函数等空间权重函数来赋值，而圆圈以外的县，空间权重都赋为0。

不喜欢固定带宽窗口的同学也可以不用它，而是用符合Tobler地理学第一定律的非固定带宽邻域统计，操作方法是在gwr.sel里面加一个命令adapt = TRUE，这样的情况下，根据你设置的k邻居数，每一轮统计的时候，和本轮对象在k以内相邻的多边形的权重参数会被赋值为0到1之间的一个数，比如下图:

我在对佐治亚做GWR，这一轮的regression target是红色的这个县，那么图上标为1的县就是红色县的1阶邻县，标为2的是2阶（邻县的邻县），标为3的是3阶（邻县的邻县的邻县）。如果用非固定带宽邻域统计，k为3，那么1、2、3都被定义为红色县的邻县，它们的权重从3到1依次增加，会按比例被赋上0和1之间的值，而其它没有标注的县，权重为0。

下一步就是用前一步做出的窗口去扫过全局区域：

gwr_result <- gwr(AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
这一步如果数据量大，可能会要跑一阵，跑完之后我们看看结果里面有什么：

Call:
gwr(formula = AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 205880.3
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 7.3883e+01 2.1081e+02 3.2802e+02 6.6691e+02 8.5705e+03 625.5656
POP_SQMI -8.0085e+01 -4.5983e-01 -1.4704e-01 -7.3703e-02 -2.1859e-03 -0.0426
Number of data points: 3108
Effective number of parameters (resial: 2traceS - traceS'S): 119.6193
Effective degrees of freedom (resial: 2traceS - traceS'S): 2988.381
Sigma (resial: 2traceS - traceS'S): 1048.78
Effective number of parameters (model: traceS): 84.90185
Effective degrees of freedom (model: traceS): 3023.098
Sigma (model: traceS): 1042.741
Sigma (ML): 1028.4
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 52109.55
AIC (GWR p. 96, eq. 4.22): 52017.7
Resial sum of squares: 3287040139
Quasi-global R2: 0.4829366
基本上你做GWR该需要的结果这里都有了。比如窗口大小（Fixed bandwidth）是205880.3，意思是前一步构建的带宽窗口是半径205.88千米的圆。Effective number of parameters显示的是你带宽窗口的大小合不合适。Sigma是残差的标准差，这个值要尽量小。Resial sum of squares（RSS）也是对拟合程度的一个评估值。最重要的是最后那个R2，越靠近1说明统计的拟合度越好。我这里面Sigma很大，R2也不是很大，因为我这里只是呈现方法，用的数据本来就是互不相干、没什么太大意义的，所以不用太纠结。如果你是真正的统计数据要来做GWR，就需要注意这些值了。

然后，我们就可以把每个县的R2画在地图上。首先，前面报告里的这些数据，比如R2，要先自己去生成的GWR结果里面去找，然后自己再算一下每个县的local R2，并把它们赋值到shp里面去：

⑵ 悬赏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$类别=="恶性"))
a;b #输出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$resial进行白噪声检验
for(i in 1:2) print(Box.test(x.fit$resial,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

⑶ 怎么学习用 R 语言进行数据挖掘

什么是R语言？应该如何开始学习/使用R语言呢？

学习R有几个月了，总算是摸着了一点门道。
写一些自己的心得和经验，方便自己进一步鼓捣R。如果有人看到我写的东西而得到了帮助，那就更好了。
什么是R？R的优点何在？
R是一个数据分析软件。简单点说，R可以看做MATLAB的“替代品”，而且具有免费开源的优势。R可以像MATLAB一样解决有关数值计算的问题，而且具有强大的数据处理，绘图功能。
R拥有大量的统计分析工具包，我的感觉是——只有我们没听说过的工具，绝对没有R没有的工具包。配合着各种各样的工具包，你可以毁灭任何关于数据和统计的问题。因为数据包的数量庞大，所以查找自己需要的数据包，可能很烦恼。
如果有以下技能，学R会很方便：
1.已经了解些高级程序语言（非常重要）
2.英语不坏
3.概率统计理论基础
4.看数据不头疼
5.看cmd or terminal 也不头疼
你需要一本适合你的R语言教材
我开始学习R的时候，找到了这个帖子

非常强大的关于R语言教材综述。我非常感谢原帖作者。你可以参考这个帖子选一本适合你的教材。
我这里在说一下我主要使用的几本教材的心得：
1. 统计建模与R软件（薛毅着）：非常优秀的R语言入门教材，涵盖了所有R的基础应用&方法，示例代码也很优秀。作为一本中文的程序语言教材，绝对是最优秀的之一。但是要看懂这本书，还是需要“已经了解些高级程序语言”。PS：我亲爱的吉林大学图书馆，有两本该教材流通，我常年霸占一本。
2. R in Nutshell：从讲解内容上看，与上一本差别不大，在R语言的应用上都是比较初级的入门，但是有些R软件&语言上的特性，写得比薛毅老师的教材深刻。这本书最大的优点就是工具书，方便开始入门时候，对有些“模棱两可”的东西的查询。PS：我将这本书打印了出来，简单的从头到尾翻过，最大的用途就是像一本字典一样查询。察掘
3. ggplot2 Elegant Graphics for Data：这是一本介绍如何使用ggplot2包，进行绘图的书。ggplot2包，非常强大的绘图工具，几乎可以操作任何图中的元素，而且是提供添加图层的方式让我们可以一步步的作图。提到ggplot2包，应该提到悄没孙一个词——“潜力无穷”，每一个介绍
ggplot2的人，都会用这个形容词。这本书最大的作用也是当做一本绘图相关的工具书，书中讲解详细，细致，每个小参数的变动都会配图帮你理解。PS：这本书我也打印出来了，非常适合查询。
几个可以逐步提高R能力的网站
1.R-bloggers: 这里有关于R和数据的一切讨论，前沿的问题，基础的问题，应有尽有。可以说这些家伙们让R变得越来越强大。我RSS了这个网站，每天都看一下有什么我感兴趣的方法和话题，慢慢的积累一些知识，是一个很有意思的过程。
2.统计之都： 这是一个有大量R使用者交流的论坛，你可以上去提问题，总有好心人来帮助你的。
3.R客： 是关于R的一个博客，更新不快，偏重国内R的一些发展。
R的使用环境
如果你看见terminal or cmd就打怵的话，一定要使用Rstudio。Rstudio的优点是，集成了Rconsole、脚本编辑器、可视化的数据查询、历史命令、帮助查询等，还有的完美的脚本和console的互动。毕竟是可视化的界面，有许多按钮可以用。R 的脚本编辑器很蛋疼，就比记事本多了个颜色高亮吧，不适合编写脚本，但适合调试脚本。
最后，说一下，刚开始学习R或者其他什么语言，都有一个通病，就是一些小细节的不知道，或者是记得不清楚，往往一个蛋疼的bug就可以耗掉大量的时间，这是一个让人想砸电脑的过程。我往后，会在博客里记录一些让我蛋很疼的小细节。本文分为6个部分，分别介绍初级入门，高级入门，绘图与可视化，计量经济学，时间序列分析，金融等。
1.初级入门
《An Introction to R》，这是官方的入门小册子。其有中文版，由丁国徽翻译，译名为《R导论》。《R4Beginners》，这本小册子有中文版应该叫《R入门》。除此之外，还可以去读刘思喆的《153分钟学会R》。这本书收集了R初学者提问频率最高的153个问题。为什么叫153分钟呢？因为最初作者写了153个问题，阅读一个问题花费1分钟时间，全局下来也就启链是153分钟了。有了这些基础之后，要去读一些经典书籍比较全面的入门书籍，比如《统计建模与R软件》，国外还有《R Cookbook》和《R in action》，本人没有看过，因此不便评论。
最后推荐，《R in a Nutshell》。对，“果壳里面的R”！当然，是开玩笑的，in a Nutshell是俚语，意思大致是“简单的说”。目前，我们正在翻译这本书的中文版，大概明年三月份交稿！这本书很不错，大家可以从现在开始期待，并广而告知一下！
2.高级入门
读了上述书籍之后，你就可以去高级入门阶段了。这时候要读的书有两本很经典的。《Statistics with R》和《The R book》。之所以说这两本书高级，是因为这两本书已经不再限于R基础了，而是结合了数据分析的各种常见方法来写就的，比较系统的介绍了R在线性回归、方差分析、多元统计、R绘图、时间序列分析、数据挖掘等各方面的内容，看完之后你会发现，哇，原来R能做的事情这么多，而且做起来是那么简洁。读到这里已经差不多了，剩下的估计就是你要专门攻读的某个方面内容了。下面大致说一说。
3.绘图与可视化
亚里斯多德说，“较其他感觉而言，人类更喜欢观看”。因此，绘图和可视化得到很多人的关注和重视。那么，如何学习R画图和数据可视化呢？再简单些，如何画直方图？如何往直方图上添加密度曲线呢？我想读完下面这几本书你就大致会明白了。
首先，画图入门可以读《R Graphics》，个人认为这本是比较经典的，全面介绍了R中绘图系统。该书对应的有一个网站，google之就可以了。更深入的可以读《Lattice：Multivariate Data Visualization with R》。上面这些都是比较普通的。当然，有比较文艺和优雅的——ggplot2系统，看《ggplot2：Elegant Graphics for Data Analysis》。还有数据挖掘方面的书：《Data Mining with Rattle and R》，主要是用Rattle软件，个人比较喜欢Rattle!当然，Rattle不是最好的，Rweka也很棒！再有就是交互图形的书了，着名的交互系统是ggobi，这个我已经喜欢两年多了，关于ggobi的书有《Interactive and Dynamic Graphics for Data Analysis With R and GGobi》，不过，也只是适宜入门，更多更全面的还是去ggobi的主页吧，上面有各种资料以及包的更新信息！
特别推荐一下，中文版绘图书籍有《现代统计图形》。
4.计量经济学
关于计量经济学，首先推荐一本很薄的小册子:《Econometrics In R》，做入门用。然后，是《Applied Econometrics with R》，该书对应的R包是AER，可以安装之后配合使用，效果甚佳。计量经济学中很大一部分是关于时间序列分析的，这一块内容在下面的地方说。
5.时间序列分析
时间序列书籍的书籍分两类，一种是比较普适的书籍，典型的代表是：《Time Series Analysis and Its Applications ：with R examples》。该书介绍了各种时间序列分析的经典方法及实现各种经典方法的R代码，该书有中文版。如果不想买的话，建议去作者主页直接下载，英文版其实读起来很简单。时间序列分析中有一大块儿是关于金融时间序列分析的。这方面比较流行的书有两本《Analysis of financial time series》，这本书的最初是用的S-plus代码，不过新版已经以R代码为主了。这本书适合有时间序列分析基础和金融基础的人来看，因为书中关于时间序列分析的理论以及各种金融知识讲解的不是特别清楚，将极值理论计算VaR的部分就比较难看懂。另外一个比较有意思的是Rmetrics推出的《TimeSeriesFAQ》，这本书是金融时间序列入门的东西，讲的很基础，但是很难懂。对应的中文版有《金融时间序列分析常见问题集》，当然，目前还没有发出来。经济领域的时间序列有一种特殊的情况叫协整，很多人很关注这方面的理论，关心这个的可以看《Analysis of Integrated and Cointegrated Time Series with R》。最后，比较高级的一本书是关于小波分析的，看《Wavelet Methods in Statistics with R》。附加一点，关于时间序列聚类的书籍目前比较少见，是一个处女地，有志之士可以开垦之！
6.金融
金融的领域很广泛，如果是大金融的话，保险也要被纳入此间。用R做金融更多地需要掌握的是金融知识，只会数据分析技术意义寥寥。我觉得这些书对于懂金融、不同数据分析技术的人比较有用，只懂数据分析技术而不动金融知识的人看起来肯定如雾里看花，甚至有人会觉得金融分析比较低级。这方面比较经典的书籍有：《Advanced Topics in Analysis of Economic and Financial Data Using R》以及《Modelling Financial Time Series With S-plus》。金融产品定价之类的常常要用到随机微分方程，有一本叫《Simulation Inference Stochastic Differential Equations：with R examples》的书是关于这方面的内容的，有实例，内容还算详实!此外，是风险度量与管理类。比较经典的有《Simulation Techniques in Financial Risk Management》、《Modern Actuarial Risk Theory Using R》和《Quantitative Risk Management：Concepts, Techniques and Tools》。投资组合分析类和期权定价类可以分别看《Portfolio Optimization with R》和《Option Pricing and Estimation of Financial Models with R》。
7.数据挖掘
这方面的书不多，只有《Data Mining with R:learing with case studies》。不过，R中数据挖掘方面的包已经足够多了，参考包中的帮助文档就足够了。