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中數據挖掘方面的包已經足夠多了，參考包中的幫助文檔就足夠了。