R语言基础应用
http://my.oschina.net/u/236959/blog/93731
1 什么是R語言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)計分析,如果感興趣的話可以到 官方網(wǎng)站了解關(guān)于其功能的更多信息。
至于R語言名稱的由來則是根據(jù)兩位主要作者的首字母(Robert Gentleman and Ross Ihaka),但過于簡短的關(guān)鍵詞也造成在搜索引擎中很不容易找到相關(guān)的資料。不過這個專門的 搜索網(wǎng)站可以幫到你。
2 為什么要學(xué)習(xí)R語言
可能你想說,“我已經(jīng)學(xué)會了spss/sas/stata...,為什么還要去學(xué)習(xí)R呢?”如下幾方面可能會吸引到你:
- R是免費開源軟件:現(xiàn)在很多學(xué)術(shù)期刊都對分析軟件有版權(quán)要求,而免費的分析工具可以使你在這方面不會有什么擔(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版可直接點擊這個 連接
在ubuntu下面安裝R則更容易,在終端里頭運行如下命令即可
sudo apt-get update
sudo apt-get install r-base
此外,學(xué)習(xí)R語言時強(qiáng)烈推薦安裝Rstudio做為R的圖形界面,關(guān)于Rstudio之前的 博文有過簡單介紹,點 這里可能轉(zhuǎn)到它的官方網(wǎng)站。
4 R語言的學(xué)習(xí)方法
學(xué)習(xí)R并不是一件非常輕松的事情,初學(xué)者需要記住的就是:
- 親手鍵入代碼并理解其意義
- 在筆記里記下一些重點或心得(個人推薦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子論壇?(免費資料也不少)
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?(需)
如果有一些簡單的入門問題,也可以在推特上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)計方法;方差分析;線性回歸和一般線性模型。
R語言基礎(chǔ)入門之二:數(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)計表格,每一行都代表一個樣本點,而每一列則代表了樣本的不同屬性或特征。初學(xué)者需要掌握的基本操作方法就是dataframe的編輯、抽取和運算。
盡管建議初學(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ù),因為Sepal.Width變量是第2個變量,所以此時鍵入下面的命令即可:
newdata=data[20:30,2] 如果需要抽取所有數(shù)據(jù)的Sepal.Width變量,那么下面兩個命令是等價的:
newdata=data[,2]
newdata=data$Sepal.Width 第三種情況是需要對數(shù)據(jù)進(jìn)行一些運算,例如需要將所有樣本的Sepal.Width變量都放大10倍,我們先將原數(shù)據(jù)進(jìn)行一個復(fù)制,再用$符號來提取運算對象即可:
newdata=data
newdata$Sepal.Width=newdata$Sepal.Width*10 3 描述統(tǒng)計
描述統(tǒng)計是一種從大量數(shù)據(jù)中壓縮提取信息的工具,最常用的就是 summary命令,運行summary(data)得到結(jié)果如下:對于數(shù)值變量計算了五個分位點和均值,對于分類變量則計算了頻數(shù)。
也可以單獨計算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ù),則可以通過散點圖來觀察規(guī)律
plot(data$Sepal.Width,Sepal.Length) 如果需要保存繪圖結(jié)果,建議使用Rstudio中的plot菜單命令,選擇save plot as image
R語言基礎(chǔ)入門之三:常用統(tǒng)計函數(shù)運算
在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,FUN=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,FUN=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)
舉例來講,我們求標(biāo)準(zhǔn)正態(tài)分布曲線下小于1的面積p(x<1),正態(tài)分布是norm,而分布函數(shù)是p,那么使用pnorm(1)就得出了結(jié)果0.84;若計算扔10次硬幣實驗中有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ù),另一部分作為檢驗數(shù)據(jù)。
R語言基礎(chǔ)入門之四:常用的統(tǒng)計推斷
通常一個研究項目能夠獲得的數(shù)據(jù)是有限的,以有限的樣本特征來推斷總體特征就稱為統(tǒng)計推斷。推斷又可細(xì)分為區(qū)間估計和假設(shè)檢驗,二者雖有區(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è)檢驗
還是以上面的X數(shù)據(jù)作為對象,來檢驗總體均值是否為10
t.test(x,mu=10,alternative='two.sided') ?#這里的原假設(shè)是總體均值(mu)為10,使用雙側(cè)檢驗,得到P值為0.17,可見P值不夠小,不能夠拒絕原假設(shè)。 T檢驗是極為常用的檢驗方法,除了單樣本推斷之外,t.test命令還可以實現(xiàn)兩樣本推斷和配對樣本推斷。如果要對總體比率或總體方差進(jìn)行推斷,可以使用 prop.test和 var.test。
3 正態(tài)分布檢驗
T檢驗的前提條件是總體服從正態(tài)分布,因此我們有必要先檢驗正態(tài)性。而且在評價回歸模型時,對殘差也需要檢驗正態(tài)性。檢驗正態(tài)性的函數(shù)是 shapiro.test
shapiro.test(x) 結(jié)果所下:
Shapiro-Wilk normality test
data: ?x
W = 0.9863, p-value = 0.8265 該檢驗的原假設(shè)是服從正態(tài)分布,由P值為0.82可判斷不能拒絕總體服從正態(tài)的假設(shè)
4 非參數(shù)檢驗
如果總體不服從正態(tài)分布,那么T檢驗就不再適用,此時我們可以利用非參數(shù)方法推斷中位數(shù)。 wilcoxon.test函數(shù)可實現(xiàn)符號秩檢驗。
wilcox.test(x,conf.int=T) ? #指定conf.int讓函數(shù)返回中位數(shù)的置信區(qū)間
wilcox.test(x,mu=1) ?#指定mu讓函數(shù)返回中位數(shù)為10的檢驗結(jié)果 5 獨立性檢驗(聯(lián)列表檢驗)
卡方分布有一個重要應(yīng)用就是根據(jù)樣本數(shù)據(jù)來檢驗兩個分類變量的獨立性,我們以CO2數(shù)據(jù)為例來說明 chisq.test函數(shù)的使用,help(CO2)可以了解更多信息。
data(CO2) ? #讀入內(nèi)置的數(shù)據(jù)包,其中Type和Treatmen是其中兩個分類變量。
chisq.test(table(CO2$Type,CO2$Treatment)) ? #使用卡方檢驗函數(shù)來檢驗這兩個因子之間是否獨立 結(jié)果顯示P值為0.82,因此可以認(rèn)為兩因子之間獨立。在樣本較小的情況下,還可以使用fisher精確檢驗,對應(yīng)的函數(shù)是 fisher.test。
R語言基礎(chǔ)入門之五:簡單線性回歸
線性回歸可能是數(shù)據(jù)分析中最為常用的工具了,如果你認(rèn)為手上的數(shù)據(jù)存在著線性定量關(guā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) 首先繪制了兩變量的散點圖,然后用lm函數(shù)建立線性回歸模型,并將回歸直線加在原圖上,最后用summary將模型結(jié)果進(jìn)行了展示,從變量P值和F統(tǒng)計量可得回歸模型是顯著的。但截距項不應(yīng)該為負(fù)數(shù),所以也可以用下面方法將截距強(qiáng)制為0。
model2=lm(Volume~Girth-1,data=trees) 2 模型診斷
在模型建立后會利用各種方式來檢驗?zāi)P偷恼_性,對殘差進(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 變量變換
因為31號樣本有著高影響力,為了降低其影響,一種方法就是將變量進(jìn)行開方變換來改善回歸結(jié)果,從殘差標(biāo)準(zhǔn)誤到殘差圖,各項觀察都說明變換是有效的。
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語言基礎(chǔ)入門之六(完):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做為因變量,因為研究的重點在于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),FUN=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')
?
與50位技術(shù)專家面對面20年技術(shù)見證,附贈技術(shù)全景圖總結(jié)
- 上一篇: ftpClient的连接超时设置(set
- 下一篇: R语言笔记:快速入门