電子發燒友App
R語言入門教程四：常用的統計推斷
通常一個研究項目能夠獲得的數據是有限的，以有限的樣本特征來推斷總體特征就稱為統計推斷。推斷又可細分為區間估計和假設檢驗，二者雖有區別，但卻是一枚硬幣的兩面，之間有著緊密的關聯。
1 對總體均值進行區間估計
假設我們從總體中抽得一個樣本，希望根據樣本均值判斷總體均值的置信區間，如下例所示：
x=rnorm（50，mean=10，sd=5） #隨機生成50個均值為10，標準差為5的隨機數為作為研究對象
mean（x）-qt（0.975，49）*sd（x）/sqrt（50） #根據統計學區間估計公式，得到95%置信度下的區間下界
mean（x）+qt（0.975，49）*sd（x）/sqrt（50） #95%置信度下的區間上界
也可以直接利用R語言內置函數t.test
t.test（x，conf.level=0.95）
從如下結果可得95%置信區間為（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 對總體均值進行假設檢驗
還是以上面的X數據作為對象，來檢驗總體均值是否為10
t.test（x，mu=10，alternative=‘two.sided’） #這里的原假設是總體均值（mu）為10，使用雙側檢驗，得到P值為0.17，可見P值不夠小，不能夠拒絕原假設。
T檢驗是極為常用的檢驗方法，除了單樣本推斷之外，t.test命令還可以實現兩樣本推斷和配對樣本推斷。如果要對總體比率或總體方差進行推斷，可以使用prop.test和var.test。
3 正態分布檢驗
T檢驗的前提條件是總體服從正態分布，因此我們有必要先檢驗正態性。而且在評價回歸模型時，對殘差也需要檢驗正態性。檢驗正態性的函數是shapiro.test
shapiro.test（x）
結果所下：
Shapiro-Wilk normality test
data： x
W = 0.9863， p-value = 0.8265
該檢驗的原假設是服從正態分布，由P值為0.82可判斷不能拒絕總體服從正態的假設
4 非參數檢驗
如果總體不服從正態分布，那么T檢驗就不再適用，此時我們可以利用非參數方法推斷中位數。wilcoxon.test函數可實現符號秩檢驗。
wilcox.test（x，conf.int=T） #指定conf.int讓函數返回中位數的置信區間
wilcox.test（x，mu=1） #指定mu讓函數返回中位數為10的檢驗結果
5 獨立性檢驗（聯列表檢驗）
卡方分布有一個重要應用就是根據樣本數據來檢驗兩個分類變量的獨立性，我們以CO2數據為例來說明chisq.test函數的使用，help（CO2）可以了解更多信息。
data（CO2） #讀入內置的數據包，其中Type和Treatmen是其中兩個分類變量。
chisq.test（table（CO2$Type，CO2$Treatment）） #使用卡方檢驗函數來檢驗這兩個因子之間是否獨立
結果顯示P值為0.82，因此可以認為兩因子之間獨立。在樣本較小的情況下，還可以使用fisher精確檢驗，對應的函數是fisher.test。
R語言入門教程五：簡單線性回歸
線性回歸可能是數據分析中最為常用的工具了，如果你認為手上的數據存在著線性定量關系，不妨先畫個散點圖觀察一下，然后用線性回歸加以分析。下面簡單介紹一下如何在R中進行線性回歸。
1 回歸建模
我們利用R語言中內置的trees數據，其中包含了Volume（體積）、Girth（樹圍）、Height（樹高）這三個變量，我們希望以體積為因變量，樹圍為自變量進行線性回歸。
plot（Volume~Girth，data=trees，pch=16，col=‘red’）
model=lm（Volume~Girth，data=trees）
abline（model，lty=2）
summary（model）
首先繪制了兩變量的散點圖，然后用lm函數建立線性回歸模型，并將回歸直線加在原圖上，最后用summary將模型結果進行了展示，從變量P值和F統計量可得回歸模型是顯著的。但截距項不應該為負數，所以也可以用下面方法將截距強制為0。
model2=lm（Volume~Girth-1，data=trees）
2 模型診斷
在模型建立后會利用各種方式來檢驗模型的正確性，對殘差進行分析是常見的方法，下面我們來生成四種用于模型診斷的圖形。
par（mfrow=c（2，2））
plot（model）
par（mfrow=c（1，1））
這里左上圖是殘差對擬合值作圖，整體呈現出一種先下降后下升的模式，顯示殘差中可能還存在未提煉出來的影響因素。右上圖殘差QQ圖，用以觀察殘差是否符合正態分布。左下圖是標準化殘差對擬合值，用于判斷模型殘差是否等方差。右下圖是標準化殘差對杠桿值，虛線表示的cooks距離等高線。我們發現31號樣本有較大的影響。
3 變量變換
因為31號樣本有著高影響力，為了降低其影響，一種方法就是將變量進行開方變換來改善回歸結果，從殘差標準誤到殘差圖，各項觀察都說明變換是有效的。
plot（sqrt（Volume）~Girth，data=trees，pch=16，col=‘red’）
model2=lm（sqrt（Volume）~Girth，data=trees）
abline（model2，lty=2）
summary（model2）
4 模型預測
下面根據上述模型計算預測值以及置信區間，predict函數可以獲得模型的預測值，加入參數可以得到預測區間
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）
我們還可以將樹圍和樹高都加入到模型中去，進行多元回歸。如果要考慮的變量很多，可以用step函數進行變量篩選，它是以AIC作為評價指標來判斷一個變量是否應該加入模型，建議使用這種自動判斷函數時要謹慎。對于嵌套模型，還可以使用anova建立方差分析表來比較模型。對于變量變換的形式，則可以使用MASS擴展包中的boxcox函數來進行COX變換。
R語言入門教程六（完）：Logistic回歸
讓我們用logistic回歸來結束本系列的內容吧，本文用例來自于John Maindonald所著的《Data Analysis and Graphics Using R》一書，其中所用的數據集是anesthetic，數據集來自于一組醫學數據，其中變量conc表示麻醉劑的用量，move則表示手術病人是否有所移動，而我們用nomove做為因變量，因為研究的重點在于conc的增加是否會使nomove的概率增加。
首先載入數據集并讀取部分文件，為了觀察兩個變量之間關系，我們可以利cdplot函數來繪制條件密度圖。
library（DAAG）
head（anesthetic）
cdplot（factor（nomove）~conc，data=anesthetic，main=‘條件密度圖’，ylab=‘病人移動’，xlab=‘麻醉劑量’）
從圖中可見，隨著麻醉劑量加大，手術病人傾向于靜止。下面利用logistic回歸進行建模，得到intercept和conc的系數為-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數據進行建模，即每一行數據均表示一個個體，另一種是使用匯總數據進行建模，先將原始數據按下面步驟進行匯總
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
得到匯總數據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
對于匯總數據，有兩種方法可以得到同樣的結果，一種是將兩種結果的向量合并做為因變量，如anes2模型。另一種是將比率做為因變量，總量做為權重進行建模，如anes3模型。這兩種建模結果是一樣的。
anes2=glm（cbind（nomove，move）~conc，family=binomial（link=‘logit’），data=anestot）
anes3=glm（prop~conc，family=binomial（link=‘logit’），weights=total，data=anestot）
根據logistic模型，我們可以使用predict函數來預測結果，下面根據上述模型來繪圖
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’）
工商網監
評論