R語言入門教程四：常用的統(tǒng)計(jì)推斷

通常一個(gè)研究項(xiàng)目能夠獲得的數(shù)據(jù)是有限的，以有限的樣本特征來推斷總體特征就稱為統(tǒng)計(jì)推斷。推斷又可細(xì)分為區(qū)間估計(jì)和假設(shè)檢驗(yàn)，二者雖有區(qū)別，但卻是一枚硬幣的兩面，之間有著緊密的關(guān)聯(lián)。

1 對總體均值進(jìn)行區(qū)間估計(jì)

假設(shè)我們從總體中抽得一個(gè)樣本，希望根據(jù)樣本均值判斷總體均值的置信區(qū)間，如下例所示：

x=rnorm（50，mean=10，sd=5） #隨機(jī)生成50個(gè)均值為10，標(biāo)準(zhǔn)差為5的隨機(jī)數(shù)為作為研究對象

mean（x）-qt（0.975，49）*sd（x）/sqrt（50） #根據(jù)統(tǒng)計(jì)學(xué)區(qū)間估計(jì)公式，得到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)性。而且在評價(jià)回歸模型時(shí)，對殘差也需要檢驗(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ù)方法推斷中位數(shù)。wilcoxon.test函數(shù)可實(shí)現(xiàn)符號(hào)秩檢驗(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)）

卡方分布有一個(gè)重要應(yīng)用就是根據(jù)樣本數(shù)據(jù)來檢驗(yàn)兩個(gè)分類變量的獨(dú)立性，我們以CO2數(shù)據(jù)為例來說明chisq.test函數(shù)的使用，help（CO2）可以了解更多信息。

data（CO2） #讀入內(nèi)置的數(shù)據(jù)包，其中Type和Treatmen是其中兩個(gè)分類變量。

chisq.test（table（CO2$Type，CO2$Treatment）） #使用卡方檢驗(yàn)函數(shù)來檢驗(yàn)這兩個(gè)因子之間是否獨(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)系，不妨先畫個(gè)散點(diǎn)圖觀察一下，然后用線性回歸加以分析。下面簡單介紹一下如何在R中進(jìn)行線性回歸。

1 回歸建模

我們利用R語言中內(nèi)置的trees數(shù)據(jù)，其中包含了Volume（體積）、Girth（樹圍）、Height（樹高）這三個(gè)變量，我們希望以體積為因變量，樹圍為自變量進(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)計(jì)量可得回歸模型是顯著的。但截距項(xiàng)不應(yīng)該為負(fù)數(shù)，所以也可以用下面方法將截距強(qiáng)制為0。

model2=lm（Volume~Girth-1，data=trees）

2 模型診斷

在模型建立后會(huì)利用各種方式來檢驗(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號(hào)樣本有較大的影響。

3 變量變換

因?yàn)?1號(hào)樣本有著高影響力，為了降低其影響，一種方法就是將變量進(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ù)上述模型計(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作為評價(jià)指標(biāo)來判斷一個(gè)變量是否應(yīng)該加入模型，建議使用這種自動(dòng)判斷函數(shù)時(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ù)病人是否有所移動(dòng)，而我們用nomove做為因變量，因?yàn)檠芯康闹攸c(diǎn)在于conc的增加是否會(huì)使nomove的概率增加。

首先載入數(shù)據(jù)集并讀取部分文件，為了觀察兩個(gè)變量之間關(guān)系，我們可以利cdplot函數(shù)來繪制條件密度圖。

library（DAAG）

head（anesthetic）

cdplot（factor（nomove）~conc，data=anesthetic，main=‘條件密度圖’，ylab=‘病人移動(dòng)’，xlab=‘麻醉劑量’）

從圖中可見，隨著麻醉劑量加大，手術(shù)病人傾向于靜止。下面利用logistic回歸進(jìn)行建模，得到intercept和conc的系數(shù)為-6.47和5.57，由此可見麻醉劑量超過1.16（6.47/5.57）時(shí)，病人靜止概率超過50%。

anes1=glm（nomove~conc，family=binomial（link=‘logit’），data=anesthetic）

summary（anes1）

上面的方法是使用原始的0-1數(shù)據(jù)進(jìn)行建模，即每一行數(shù)據(jù)均表示一個(gè)個(gè)體，另一種是使用匯總數(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’）