電子發(fā)燒友App
R語言簡介
R是用于統(tǒng)計分析、繪圖的語言和操作環(huán)境。R是屬于GNU系統(tǒng)的一個自由、免費(fèi)、源代碼開放的軟件，它是一個用于統(tǒng)計計算和統(tǒng)計制圖的優(yōu)秀工具。
R語言是一個開源的數(shù)據(jù)分析環(huán)境，起初是由數(shù)位統(tǒng)計學(xué)家建立起來，以更好的進(jìn)行統(tǒng)計計算和繪圖，這篇wiki中包含了一些基本情況的介紹。由于R可以通過安裝擴(kuò)展包（Packages）而得到增強(qiáng)，所以其功能已經(jīng)遠(yuǎn)遠(yuǎn)不限于統(tǒng)計分析。
R語言的特點(diǎn)
R作為一種統(tǒng)計分析軟件，是集統(tǒng)計分析與圖形顯示于一體的。它可以運(yùn)行于UNIX，Windows和Macintosh的操作系統(tǒng)上，而且嵌入了一個非常方便實(shí)用的幫助系統(tǒng)，相比于其他統(tǒng)計分析軟件，R還有以下特點(diǎn)：1.R是自由軟件。這意味著它是完全免費(fèi)，開放源代碼的。可以在它的網(wǎng)站及其鏡像中下載任何有關(guān)的安裝程序、源代碼、程序包及其源代碼、文檔資料。標(biāo)準(zhǔn)的安裝文件身自身就帶有許多模塊和內(nèi)嵌統(tǒng)計函數(shù)，安裝好后可以直接實(shí)現(xiàn)許多常用的統(tǒng)計功能。
2.R是一種可編程的語言。作為一個開放的統(tǒng)計編程環(huán)境，語法通俗易懂，很容易學(xué)會和掌握語言的語法。而且學(xué)會之后，我們可以編制自己的函數(shù)來擴(kuò)展現(xiàn)有的語言。這也就是為什么它的更新速度比一般統(tǒng)計軟件，如，SPSS，SAS等快得多。大多數(shù)最新的統(tǒng)計方法和技術(shù)都可以在R中直接得到。
3. 所有R的函數(shù)和數(shù)據(jù)集是保存在程序包里面的。只有當(dāng)一個包被載入時，它的內(nèi)容才可以被訪問。一些常用、基本的程序包已經(jīng)被收入了標(biāo)準(zhǔn)安裝文件中，隨著新的統(tǒng)計分析方法的出現(xiàn)，標(biāo)準(zhǔn)安裝文件中所包含的程序包也隨著版本的更新而不斷變化。在另外版安裝文件中，已經(jīng)包含的程序包有：base一R的基礎(chǔ)模塊、mle一極大似然估計模塊、ts一時間序列分析模塊、mva一多元統(tǒng)計分析模塊、survival一生存分析模塊等等。
4.R具有很強(qiáng)的互動性。除了圖形輸出是在另外的窗口處，它的輸入輸出窗口都是在同一個窗口進(jìn)行的，輸入語法中如果出現(xiàn)錯誤會馬上在窗口口中得到提示，對以前輸入過的命令有記憶功能，可以隨時再現(xiàn)、編輯修改以滿足用戶的需要。輸出的圖形可以直接保存為JPG，BMP，PNG等圖片格式，還可以直接保存為PDF文件。另外，和其他編程語言和數(shù)據(jù)庫之間有很好的接口。［2］ 5.如果加入R的幫助郵件列表一，每天都可能會收到幾十份關(guān)于R的郵件資訊?？梢院腿蛞涣鞯慕y(tǒng)計計算方面的專家討論各種問題，可以說是全世界最大、最前沿的統(tǒng)計學(xué)家思維的聚集地。
R是基于S語言的一個GNU項(xiàng)目，所以也可以當(dāng)作S語言的一種實(shí)現(xiàn)，通常用S語言編寫的代碼都可以不作修改的在R環(huán)境下運(yùn)行。 R的語法是來自Scheme。R的使用與S-PLUS有很多類似之處，這兩種語言有一定的兼容性。S-PLUS的使用手冊，只要稍加修改就可作為R的使用手冊。所以有人說：R，是S-PLUS的一個“克隆”。但是請不要忘了：R是免費(fèi)的（R is free）。R語言源代碼托管在github，具體地址可以看參考資料。
R語言的下載可以通過CRAN的鏡像來查找。
R語言有域名為.cn的下載地址，有六個，其中兩個由Datagurn，由中國科學(xué)技術(shù)大學(xué)提供的。R語言Windows版，其中由兩個下載地點(diǎn)是Datagurn和USTC提供的。
R語言基礎(chǔ)入門教程一：
1、學(xué)習(xí)前提
在繼續(xù)學(xué)習(xí)本教程之前，您應(yīng)該基本了解計算機(jī)編程術(shù)語。 對任何編程語言的基本理解將幫助您理解R語言編程概念，并在學(xué)習(xí)軌道上快速移動
R語言適用人群
本教程是為期待使用R編程開發(fā)統(tǒng)計軟件的軟件程序員，統(tǒng)計學(xué)家和數(shù)據(jù)挖掘者設(shè)計的。 如果你試圖理解R編程語言作為一個初學(xué)者，本教程將給你足夠的了解語言的幾乎所有的概念，從那里你可以把自己的更高水平的專業(yè)知識。
2 為什么要學(xué)習(xí)R語言
可能你想說，“我已經(jīng)學(xué)會了spss/sas/stata.。。，為什么還要去學(xué)習(xí)R呢？”如下幾方面可能會吸引到你：
R是免費(fèi)開源軟件：現(xiàn)在很多學(xué)術(shù)期刊都對分析軟件有版權(quán)要求，而免費(fèi)的分析工具可以使你在這方面不會有什么擔(dān)心。另一方面，如果學(xué)術(shù)界出現(xiàn)一種新的數(shù)據(jù)分析方法，那么要過很長一段時間才會出現(xiàn)在商業(yè)軟件中。但開源軟件的好處就在于，很快就會有人將這種方法編寫成擴(kuò)展包，或者你自己就可以做這件工作。
命令行工作方式：許多人喜歡類似SPSS菜單式的操作，這對于初學(xué)者來說很方便入門，但對于數(shù)據(jù)分析來說，命令行操作會更加的靈活，更容易進(jìn)行編程和自動化處理。而且命令行操作會更容易?？?，不是嘛，一般人看到你在狂敲一推代碼后得到一個分析結(jié)果，對你投來的目光是會不一樣的。
小巧而精悍：R語言的安裝包更小，大約不到40M，相比其它幾個大家伙它算是非常小巧精悍了。目前R語言非常受到專業(yè)人士歡迎，根據(jù)對數(shù)據(jù)挖掘大賽勝出者的調(diào)查可以發(fā)現(xiàn)，他們用的工具基本上都是R語言。此外，從最近幾次R語言大會上可以了解到，咨詢業(yè)、金融業(yè)、醫(yī)藥業(yè)都在大量的使用R語言，包括google/facebook的大公司都在用它。因此，學(xué)習(xí)R語言對你的職業(yè)發(fā)展一定是有幫助的。
3 R語言的下載和GUI界面
R語言安裝包可以在官方網(wǎng)站下載，windows版可直接點(diǎn)擊這個連接
在ubuntu下面安裝R則更容易，在終端里頭運(yùn)行如下命令即可
sudo apt-get update
sudo apt-get install r-base
此外，學(xué)習(xí)R語言時強(qiáng)烈推薦安裝Rstudio做為R的圖形界面，關(guān)于Rstudio之前的博文有過簡單介紹，點(diǎn)這里可能轉(zhuǎn)到它的官方網(wǎng)站。
4 R語言的學(xué)習(xí)方法
學(xué)習(xí)R并不是一件非常輕松的事情，初學(xué)者需要記住的就是：
親手鍵入代碼并理解其意義
在筆記里記下一些重點(diǎn)或心得（個人推薦Evernote）
堅持練習(xí)，對手邊的數(shù)據(jù)進(jìn)行應(yīng)用分析
理解背景知識，細(xì)節(jié)很重要。
5 哪里可以得到參考資料
1.官方網(wǎng)站 http://cran.csdb.cn/index.html （官方文獻(xiàn)集中地）
2.統(tǒng)計之都論壇
3.人大經(jīng)濟(jì)論壇－R子論壇 （免費(fèi)資料也不少）
4.http://library.nu/ 這是網(wǎng)上電子書最多的地方，其中有一個R語言專門書柜（也就是一個shelves）
5.關(guān)于R語言的教材小結(jié)
6.筆者在verycd上發(fā)的一個書單
7.一個國外著名的R語言群博 http://www.r-bloggers.com/
8.展示R語言的各類繪圖 http://addictedtor.free.fr/graphiques/
本人博客里也有一些關(guān)于R語言的資料：xccds1977.blogspot.com （需FQ）
如果有一些簡單的入門問題，也可以在推特上follow me twitter： @xccds
6 本系列博文的目的
本系列入門的目的是為初學(xué)者提供最簡潔清晰的資料，以迅速入門。所針對的讀者人群是那些正在大學(xué)里學(xué)習(xí)初級統(tǒng)計學(xué)的同學(xué)。本系列計劃包括內(nèi)容有：基本命令，數(shù)據(jù)操作；描述統(tǒng)計和繪圖；重要的R語言函數(shù)計算；統(tǒng)計推斷和估計；非參數(shù)統(tǒng)計方法；方差分析；線性回歸和一般線性模型。
#e#
R語言入門教程二：數(shù)據(jù)導(dǎo)入和描述統(tǒng)計
1 數(shù)據(jù)導(dǎo)入
對初學(xué)者來講，面對一片空白的命令行窗口，第一道真正的難關(guān)也許就是數(shù)據(jù)的導(dǎo)入。數(shù)據(jù)導(dǎo)入有很多途徑，例如從網(wǎng)頁抓取、公共數(shù)據(jù)源獲得、文本文件導(dǎo)入。為了快速入門，建議初學(xué)者采取R語言協(xié)同Excel電子表格的方法。也就是先用較為熟悉的Excel讀取和整理你要處理的數(shù)據(jù)，然后“粘貼”到R中。
例如我們先從這個地址下載iris.csv演示數(shù)據(jù)，在Excel中打開，框選所有的樣本然后“復(fù)制”。在R語言中輸入如下命令：
data=read.table（‘clipboard’，T）
這的里read.table是R讀取外部數(shù)據(jù)的常用命令，T表示第一行是表頭信息，整個數(shù)據(jù)存在名為data的變量中。另一種更方便的導(dǎo)入方法是利用Rstudio的功能，在workspace菜單選擇“import dataset”也是一樣的。
2 Dataframe操作
在數(shù)據(jù)導(dǎo)入R語言后，會以數(shù)據(jù)框（dataframe）的形式儲存。dataframe是一種R的數(shù)據(jù)格式，可以將它想象成類似統(tǒng)計表格，每一行都代表一個樣本點(diǎn)，而每一列則代表了樣本的不同屬性或特征。初學(xué)者需要掌握的基本操作方法就是dataframe的編輯、抽取和運(yùn)算。
盡管建議初學(xué)者在Excel中就把數(shù)據(jù)處理好，但有時候還是需要在R中對數(shù)據(jù)進(jìn)行編輯，下面的命令可以讓你有機(jī)會修改數(shù)據(jù)并存入到新的變量newdata中：
newdata=edit（data）
另一種情況就是我們可能只關(guān)注數(shù)據(jù)的一部分，例如從原數(shù)據(jù)中抽取第20到30號樣本的Sepal.Width變量數(shù)據(jù)，因?yàn)镾epal.Width變量是第2個變量，所以此時鍵入下面的命令即可：
newdata=data［20:30，2］
如果需要抽取所有數(shù)據(jù)的Sepal.Width變量，那么下面兩個命令是等價的：
newdata=data［，2］
newdata=data$Sepal.Width
第三種情況是需要對數(shù)據(jù)進(jìn)行一些運(yùn)算，例如需要將所有樣本的Sepal.Width變量都放大10倍，我們先將原數(shù)據(jù)進(jìn)行一個復(fù)制，再用$符號來提取運(yùn)算對象即可：
newdata=data
newdata$Sepal.Width=newdata$Sepal.Width*10
3 描述統(tǒng)計
描述統(tǒng)計是一種從大量數(shù)據(jù)中壓縮提取信息的工具，最常用的就是summary命令，運(yùn)行summary（data）得到結(jié)果如下：對于數(shù)值變量計算了五個分位點(diǎn)和均值，對于分類變量則計算了頻數(shù)。
也可以單獨(dú)計算Sepal.Width變量的平均值和標(biāo)準(zhǔn)差
mean（data$Sepal.Width）
sd（data$Sepal.Width）
計算分類數(shù)據(jù)Species變量的頻數(shù)表和條形圖
table（data$Species）
barplot（table（data$Species））
對于一元數(shù)值數(shù)據(jù)，繪制直方圖和箱線圖觀察其分布是常用的方法：
hist（data$Sepal.Width）
boxplot（data$Sepal.Width）
對于二元數(shù)值數(shù)據(jù)，則可以通過散點(diǎn)圖來觀察規(guī)律
plot（data$Sepal.Width，Sepal.Length）
如果需要保存繪圖結(jié)果，建議使用Rstudio中的plot菜單命令，選擇save plot as image
R語言入門教程三：常用統(tǒng)計函數(shù)運(yùn)算
在R語言中經(jīng)常會用到函數(shù)，例如上節(jié)中講到的求樣本統(tǒng)計量就需要均值函數(shù)（mean）和標(biāo)準(zhǔn)差函數(shù)（sd）。對于二元數(shù)值數(shù)據(jù)還用到協(xié)方差（cov），對于二元分類數(shù)據(jù)則可以用交叉聯(lián)列表函數(shù)（table）。下文講述在初級統(tǒng)計學(xué)中最常用到的三類函數(shù)。
一、數(shù)據(jù)匯總函數(shù)
我們還是以R中自帶的iris數(shù)據(jù)為例，輸入head（iris）你可以獲得數(shù)據(jù)的前6個樣本及對應(yīng)的5個變量。取出最后兩列數(shù)據(jù)作為講解的對象：Species表示花的種類，Petal.Width表示花瓣寬度
data=iris［，c（4，5）］
下一步我們想計算不同種類花瓣的平均寬度，可以使用tapply函數(shù)，在計算前先用attach命令將data這個數(shù)據(jù)框解包以方便直接操作其變量，而不需再用$符號。
attach（data）
tapply（X=Petal.Width，INDEX=Species，F(xiàn)UN=mean）
結(jié)果如下
setosa versicolor virginica
0.246 1.326 2.026
和tapply類似的還有sapply函數(shù)，在進(jìn)一步講解前初學(xué)者還需搞清楚兩種數(shù)據(jù)表現(xiàn)方式，即stack（堆疊數(shù)據(jù)）和unstack（非堆疊數(shù)據(jù)），上面的data就是一個堆疊數(shù)據(jù)，每一行表示一個樣本。而非堆疊數(shù)據(jù)可以根據(jù)unstack函數(shù)轉(zhuǎn)換而來
data.unstack=unstack（data）
head（data.unstack）
你應(yīng)該明白這二者之間的區(qū)別了，如果要對非堆疊數(shù)據(jù)計算不同種類花瓣的平均寬度，可以利用如下函數(shù)。
sapply（data.unstack，F(xiàn)UN=mean）
結(jié)果是一樣的，也就是說tapply對應(yīng)于stack數(shù)據(jù)，而sapply對應(yīng)于unstack數(shù)據(jù)
二、概率計算函數(shù)
如果給定一種概率分布，通常會有四類計算問題：
計算其概率密度density （d）
計算其概率分布probability（p）
計算其百分位數(shù)quantile （q）
隨機(jī)數(shù)模擬random （r）
記住上面四類計算對應(yīng)的英文首字母，再對照下表就很容易計算各種概率問題了。
舉例來講，我們求標(biāo)準(zhǔn)正態(tài)分布曲線下小于1的面積p（x《1），正態(tài)分布是norm，而分布函數(shù)是p，那么使用pnorm（1）就得出了結(jié)果0.84；若計算扔10次硬幣實(shí)驗(yàn)中有3次正面向上的概率，類似的dbinom（x=3，size=10，prob=0.5）得出0.11
三、抽樣函數(shù)
我們想從1到10中隨機(jī)抽取5個數(shù)字，那么這樣來做：首先產(chǎn)生一個序列，然后用sample函數(shù)進(jìn)行無放回抽取。
x=1:10
sample（x，size=5）
有放回抽取則是
sample（x，size=5，replace=T）
sample函數(shù)在建模中經(jīng)常用來對樣本數(shù)據(jù)進(jìn)行隨機(jī)的劃分，一部分作為訓(xùn)練數(shù)據(jù)，另一部分作為檢驗(yàn)數(shù)據(jù)。
#e#
R語言入門教程四：常用的統(tǒng)計推斷
通常一個研究項(xiàng)目能夠獲得的數(shù)據(jù)是有限的，以有限的樣本特征來推斷總體特征就稱為統(tǒng)計推斷。推斷又可細(xì)分為區(qū)間估計和假設(shè)檢驗(yàn)，二者雖有區(qū)別，但卻是一枚硬幣的兩面，之間有著緊密的關(guān)聯(lián)。
1 對總體均值進(jìn)行區(qū)間估計
假設(shè)我們從總體中抽得一個樣本，希望根據(jù)樣本均值判斷總體均值的置信區(qū)間，如下例所示：
x=rnorm（50，mean=10，sd=5） #隨機(jī)生成50個均值為10，標(biāo)準(zhǔn)差為5的隨機(jī)數(shù)為作為研究對象
mean（x）-qt（0.975，49）*sd（x）/sqrt（50） #根據(jù)統(tǒng)計學(xué)區(qū)間估計公式，得到95%置信度下的區(qū)間下界
mean（x）+qt（0.975，49）*sd（x）/sqrt（50） #95%置信度下的區(qū)間上界
也可以直接利用R語言內(nèi)置函數(shù)t.test
t.test（x，conf.level=0.95）
從如下結(jié)果可得95%置信區(qū)間為（9.56，12.36）
One Sample t-test
data： x
t = 15.7301， df = 49， p-value 《 2.2e-16
alternative hypothesis： true mean is not equal to 0
95 percent confidence interval：
9.563346 12.364729
sample estimates：
mean of x
10.96404
2 對總體均值進(jìn)行假設(shè)檢驗(yàn)
還是以上面的X數(shù)據(jù)作為對象，來檢驗(yàn)總體均值是否為10
t.test（x，mu=10，alternative=‘two.sided’） #這里的原假設(shè)是總體均值（mu）為10，使用雙側(cè)檢驗(yàn)，得到P值為0.17，可見P值不夠小，不能夠拒絕原假設(shè)。
T檢驗(yàn)是極為常用的檢驗(yàn)方法，除了單樣本推斷之外，t.test命令還可以實(shí)現(xiàn)兩樣本推斷和配對樣本推斷。如果要對總體比率或總體方差進(jìn)行推斷，可以使用prop.test和var.test。
3 正態(tài)分布檢驗(yàn)
T檢驗(yàn)的前提條件是總體服從正態(tài)分布，因此我們有必要先檢驗(yàn)正態(tài)性。而且在評價回歸模型時，對殘差也需要檢驗(yàn)正態(tài)性。檢驗(yàn)正態(tài)性的函數(shù)是shapiro.test
shapiro.test（x）
結(jié)果所下：
Shapiro-Wilk normality test
data： x
W = 0.9863， p-value = 0.8265
該檢驗(yàn)的原假設(shè)是服從正態(tài)分布，由P值為0.82可判斷不能拒絕總體服從正態(tài)的假設(shè)
4 非參數(shù)檢驗(yàn)
如果總體不服從正態(tài)分布，那么T檢驗(yàn)就不再適用，此時我們可以利用非參數(shù)方法推斷中位數(shù)。wilcoxon.test函數(shù)可實(shí)現(xiàn)符號秩檢驗(yàn)。
wilcox.test（x，conf.int=T） #指定conf.int讓函數(shù)返回中位數(shù)的置信區(qū)間
wilcox.test（x，mu=1） #指定mu讓函數(shù)返回中位數(shù)為10的檢驗(yàn)結(jié)果
5 獨(dú)立性檢驗(yàn)（聯(lián)列表檢驗(yàn)）
卡方分布有一個重要應(yīng)用就是根據(jù)樣本數(shù)據(jù)來檢驗(yàn)兩個分類變量的獨(dú)立性，我們以CO2數(shù)據(jù)為例來說明chisq.test函數(shù)的使用，help（CO2）可以了解更多信息。
data（CO2） #讀入內(nèi)置的數(shù)據(jù)包，其中Type和Treatmen是其中兩個分類變量。
chisq.test（table（CO2$Type，CO2$Treatment）） #使用卡方檢驗(yàn)函數(shù)來檢驗(yàn)這兩個因子之間是否獨(dú)立
結(jié)果顯示P值為0.82，因此可以認(rèn)為兩因子之間獨(dú)立。在樣本較小的情況下，還可以使用fisher精確檢驗(yàn)，對應(yīng)的函數(shù)是fisher.test。
R語言入門教程五：簡單線性回歸
線性回歸可能是數(shù)據(jù)分析中最為常用的工具了，如果你認(rèn)為手上的數(shù)據(jù)存在著線性定量關(guān)系，不妨先畫個散點(diǎn)圖觀察一下，然后用線性回歸加以分析。下面簡單介紹一下如何在R中進(jìn)行線性回歸。
1 回歸建模
我們利用R語言中內(nèi)置的trees數(shù)據(jù)，其中包含了Volume（體積）、Girth（樹圍）、Height（樹高）這三個變量，我們希望以體積為因變量，樹圍為自變量進(jìn)行線性回歸。
plot（Volume~Girth，data=trees，pch=16，col=‘red’）
model=lm（Volume~Girth，data=trees）
abline（model，lty=2）
summary（model）
首先繪制了兩變量的散點(diǎn)圖，然后用lm函數(shù)建立線性回歸模型，并將回歸直線加在原圖上，最后用summary將模型結(jié)果進(jìn)行了展示，從變量P值和F統(tǒng)計量可得回歸模型是顯著的。但截距項(xiàng)不應(yīng)該為負(fù)數(shù)，所以也可以用下面方法將截距強(qiáng)制為0。
model2=lm（Volume~Girth-1，data=trees）
2 模型診斷
在模型建立后會利用各種方式來檢驗(yàn)?zāi)Ｐ偷恼_性，對殘差進(jìn)行分析是常見的方法，下面我們來生成四種用于模型診斷的圖形。
par（mfrow=c（2，2））
plot（model）
par（mfrow=c（1，1））
這里左上圖是殘差對擬合值作圖，整體呈現(xiàn)出一種先下降后下升的模式，顯示殘差中可能還存在未提煉出來的影響因素。右上圖殘差QQ圖，用以觀察殘差是否符合正態(tài)分布。左下圖是標(biāo)準(zhǔn)化殘差對擬合值，用于判斷模型殘差是否等方差。右下圖是標(biāo)準(zhǔn)化殘差對杠桿值，虛線表示的cooks距離等高線。我們發(fā)現(xiàn)31號樣本有較大的影響。
3 變量變換
因?yàn)?1號樣本有著高影響力，為了降低其影響，一種方法就是將變量進(jìn)行開方變換來改善回歸結(jié)果，從殘差標(biāo)準(zhǔn)誤到殘差圖，各項(xiàng)觀察都說明變換是有效的。
plot（sqrt（Volume）~Girth，data=trees，pch=16，col=‘red’）
model2=lm（sqrt（Volume）~Girth，data=trees）
abline（model2，lty=2）
summary（model2）
4 模型預(yù)測
下面根據(jù)上述模型計算預(yù)測值以及置信區(qū)間，predict函數(shù)可以獲得模型的預(yù)測值，加入?yún)?shù)可以得到預(yù)測區(qū)間
plot（sqrt（Volume）~Girth，data=trees，pch=16，col=‘red’）
model2=lm（sqrt（Volume）~Girth，data=trees）
data.pre=data.frame（predict（model2，interval=‘prediction’））
lines（data.pre$lwr~trees$Girth，col=‘blue’，lty=2）
lines（data.pre$upr~trees$Girth，col=‘blue’，lty=2）
我們還可以將樹圍和樹高都加入到模型中去，進(jìn)行多元回歸。如果要考慮的變量很多，可以用step函數(shù)進(jìn)行變量篩選，它是以AIC作為評價指標(biāo)來判斷一個變量是否應(yīng)該加入模型，建議使用這種自動判斷函數(shù)時要謹(jǐn)慎。對于嵌套模型，還可以使用anova建立方差分析表來比較模型。對于變量變換的形式，則可以使用MASS擴(kuò)展包中的boxcox函數(shù)來進(jìn)行COX變換。
R語言入門教程六（完）：Logistic回歸
讓我們用logistic回歸來結(jié)束本系列的內(nèi)容吧，本文用例來自于John Maindonald所著的《Data Analysis and Graphics Using R》一書，其中所用的數(shù)據(jù)集是anesthetic，數(shù)據(jù)集來自于一組醫(yī)學(xué)數(shù)據(jù)，其中變量conc表示麻醉劑的用量，move則表示手術(shù)病人是否有所移動，而我們用nomove做為因變量，因?yàn)檠芯康闹攸c(diǎn)在于conc的增加是否會使nomove的概率增加。
首先載入數(shù)據(jù)集并讀取部分文件，為了觀察兩個變量之間關(guān)系，我們可以利cdplot函數(shù)來繪制條件密度圖。
library（DAAG）
head（anesthetic）
cdplot（factor（nomove）~conc，data=anesthetic，main=‘條件密度圖’，ylab=‘病人移動’，xlab=‘麻醉劑量’）
從圖中可見，隨著麻醉劑量加大，手術(shù)病人傾向于靜止。下面利用logistic回歸進(jìn)行建模，得到intercept和conc的系數(shù)為-6.47和5.57，由此可見麻醉劑量超過1.16（6.47/5.57）時，病人靜止概率超過50%。
anes1=glm（nomove~conc，family=binomial（link=‘logit’），data=anesthetic）
summary（anes1）
上面的方法是使用原始的0-1數(shù)據(jù)進(jìn)行建模，即每一行數(shù)據(jù)均表示一個個體，另一種是使用匯總數(shù)據(jù)進(jìn)行建模，先將原始數(shù)據(jù)按下面步驟進(jìn)行匯總
anestot=aggregate（anesthetic［，c（‘move’，‘nomove’）］，by=list（conc=anesthetic$conc），F(xiàn)UN=sum）
anestot$conc=as.numeric（as.character（anestot$conc））
anestot$total=apply（anestot［，c（‘move’，‘nomove’）］，1，sum）
anestot$prop=anestot$nomove/anestot$total
得到匯總數(shù)據(jù)anestot如下所示
conc move nomove total prop
1 0.8 6 1 7 0.1428571
2 1.0 4 1 5 0.2000000
3 1.2 2 4 6 0.6666667
4 1.4 2 4 6 0.6666667
5 1.6 0 4 4 1.0000000
6 2.5 0 2 2 1.0000000
對于匯總數(shù)據(jù)，有兩種方法可以得到同樣的結(jié)果，一種是將兩種結(jié)果的向量合并做為因變量，如anes2模型。另一種是將比率做為因變量，總量做為權(quán)重進(jìn)行建模，如anes3模型。這兩種建模結(jié)果是一樣的。
anes2=glm（cbind（nomove，move）~conc，family=binomial（link=‘logit’），data=anestot）
anes3=glm（prop~conc，family=binomial（link=‘logit’），weights=total，data=anestot）
根據(jù)logistic模型，我們可以使用predict函數(shù)來預(yù)測結(jié)果，下面根據(jù)上述模型來繪圖
x=seq（from=0，to=3，length.out=30）
y=predict（anes1，data.frame（conc=x），type=‘response’）
plot（prop~conc，pch=16，col=‘red’，data=anestot，xlim=c（0.5，3），main=‘Logistic回歸曲線圖’，ylab=‘病人靜止概率’，xlab=‘麻醉劑量’）
lines（y~x，lty=2，col=‘blue’）
工商網(wǎng)監(jiān)
評論