用R語言進(jìn)行數(shù)據(jù)分析：常規(guī)和廣義線性模型

線性模型

對于常規(guī)的多重模型（multiple model）擬合，最基本的函數(shù)是lm()。 下面是調(diào)用它的方式的一種改進(jìn)版：

     >fitted.model<- lm(formula, data =data.frame)

例如

     > fm2 <- lm(y ~ x1 + x2, data = production)

將會擬合 y 對 x1 和 x2 的多重回歸模型(和一個隱式的截距項(xiàng))。

一個重要的(技術(shù)上可選)參數(shù)是data = production。它指定任何構(gòu)建這個模型的參數(shù)首先必須來自 數(shù)據(jù)框 production。 這里不需要考慮數(shù)據(jù)框 production 是否被綁定在搜索路徑中。

廣義線性模型

廣義線性建模是線性模型在研究響應(yīng)值的非正態(tài)分布以及非線性模型的簡潔直接的線性轉(zhuǎn)化 時的一種發(fā)展。 廣義線性模型 是基于下面一系列 假設(shè)前提的：

• 有一個響應(yīng)變量 y和一系列有趣的刺激變量（stimulus variable） x_1, x_2,…。 這些刺激變量決定響應(yīng)變量的最終分布。
• 刺激變量僅僅通過一個線性函數(shù)影響響應(yīng)值y 的分布。 該線性函數(shù)稱為 線性預(yù)測器（linear predictor），常常寫成

            eta = beta_1 x_1 + beta_2 x_2 +...+ beta_p x_p,

  因此 x_i 當(dāng)且僅當(dāng) beta_i 等于0時對 y 的分布沒有影響。

• y 分布的形式為

            f_Y(y; mu, phi)
              = exp((A/phi) * (y lambda(mu) - gamma(lambda(mu))) + tau(y, phi))

  其中 phi 是度量參數(shù)（scale parameter）(可能已知)，對所有觀測 恒定；A 是一個先驗(yàn)的權(quán)重，假定知道但是 可能隨觀測不同有所不同；mu是 y 的均值。 也就是說假定 y 的分布是由 均值和一個可能的度量參數(shù)決定的。

• 均值 mu 是線性預(yù)測器的平滑可逆函數(shù)（smooth invertible function）：

            mu = m(eta),    eta = m^{-1}(mu) = ell(mu)

  該可逆函數(shù) ell() 稱為 關(guān)聯(lián)函數(shù)（link function）。

這些假定比較寬松，足以包括統(tǒng)計(jì)實(shí)踐中大多數(shù)有用的統(tǒng)計(jì)模型， 但也足夠嚴(yán)謹(jǐn)，使得可以發(fā)展計(jì)算和推論中一致的方法（ 至少可以近似一致）。 讀者如果想了解這方面最新的進(jìn)展，可以 參考 McCullagh & Nelder (1989) 或者 Dobson (1990)。

族

R 提供了一系列廣義線性建模工具，從類型上來說包括 gaussian, 二項(xiàng)式, poisson, 反 gaussian 和 gamma 模型的響應(yīng)變量分布以及 在響應(yīng)變量分布沒有明確給定時的逆似然（quasi-likelihood）模型。 在后者，方差函數(shù)（variance function） 可以認(rèn)為是均值的函數(shù)，但是在另外一些情況下， 該函數(shù)可以由響應(yīng)變量的分布得到。

每一種響應(yīng)分布允許各種關(guān)聯(lián)函數(shù)將均值和線性預(yù)測器關(guān)聯(lián)起來。 這些自動關(guān)聯(lián)函數(shù)如 下表所示：

這些用于模型構(gòu)建過程中的響應(yīng)分布，關(guān)聯(lián)函數(shù)和各種 其他必要的信息統(tǒng)稱為 廣義線性模型的族（family）。

glm()函數(shù)

既然響應(yīng)的分布僅僅通過單一的一個線性函數(shù)依賴于 刺激變量，那么用于線性模型的機(jī)制同樣 可以用于指定一個廣義模型的線性部分。 但是族必須以一種不同的方式指定。

R 用于廣義線性回歸的函數(shù)是glm()， 它的使用形式為

     >fitted.model<- glm(formula, family=family.generator, data=data.frame)

和lm()相比，唯一的一個新特性就是描述族的參數(shù)family.generator。 它是產(chǎn)生函數(shù)和表達(dá)式列表的函數(shù)名字。這些函數(shù) 用于定義和控制模型的構(gòu)建與計(jì)算過程。 盡管開始看起來有點(diǎn)復(fù)雜， 但它們非常容易 使用。

這些名字是標(biāo)準(zhǔn)的。程序給定的族生成器可以參見 Families 列表中 的“族名”。當(dāng)選擇一個關(guān)聯(lián)函數(shù)時， 該關(guān)聯(lián)函數(shù)名和族名可以同時在括弧里面作為 參數(shù)設(shè)定。在擬（quasi） 家族里面，方差函數(shù)也是以這種方式設(shè)定。

一些例子可能會使這個過程更清楚。

gaussian族

命令

     > fm <- glm(y ~ x1 + x2, family = gaussian, data = sales)

和下面的命令結(jié)果一致。

     > fm <- lm(y ~ x1+x2, data=sales)

但是效率上，前者差一點(diǎn)。注意，gaussian 族沒有相關(guān)參數(shù)， 因此它不提供關(guān)聯(lián)函數(shù)的。 如一個問題需要用非標(biāo)準(zhǔn)關(guān)聯(lián)函數(shù)的 gaussian 族， 那么只能采用我們后面討論的擬族。

二項(xiàng)式族

考慮 Silvey (1970) 提供的一個小的例子。

在 Kalythos 的 Aegean 島上，男性居民常?；加?一種先天的眼科疾病，并且隨著年齡的增長而變的愈顯著。 現(xiàn)在搜集了各種年齡段島上男性居民的樣本，同時記錄了盲眼的數(shù)目。 數(shù)據(jù)顯示如下：

我們想知道的是這些數(shù)據(jù)是否吻合 logistic 和 probit 模型， 并且分別估計(jì)各個模型的 LD50，也就是一個男性居民盲眼的概率 為50%時候的年齡。

如果 y 和 n 是年齡為 x 時的盲眼數(shù)目和檢測 樣本數(shù)目，兩種模型的形式都為 y ~ B(n, F(beta_0 + beta_1 x))， 其中在 probit 模型中， F(z) = Phi(z) 是標(biāo)準(zhǔn)的正態(tài)分布函數(shù)，而在 logit 模型 (默認(rèn))中， F(z) = e^z/(1+e^z)。 這兩種模型中， LD50 = – beta_0/beta_1 ，即分布函數(shù)的參數(shù)為0時 所在的點(diǎn)。

第一步是把數(shù)據(jù)轉(zhuǎn)換成數(shù)據(jù)框。

     > kalythos <- data.frame(x = c(20,35,45,55,70), n = rep(50,5),
                              y = c(6,17,26,37,44))

在glm()擬合二項(xiàng)式模型時，響應(yīng)變量 有三種可能性：

• 如果響應(yīng)變量是向量， 則假定操作二元（binary） 數(shù)據(jù)，因此要求是0/1向量。
• 如果響應(yīng)變量是雙列矩陣，則假定第一列 為試驗(yàn)成功的次數(shù)第二列 為試驗(yàn)失敗的次數(shù)。
• 如果響應(yīng)變量是因子，則第一水平作為失敗 (0) 考慮而其他的作為`成功'(1) 考慮。

我們采用的是第二種慣例。我們在數(shù)據(jù)框中 增加了一個矩陣：

     > kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)

為了擬合這些模型，我們采用

     > fmp <- glm(Ymat ~ x, family = binomial(link=probit), data = kalythos)
     > fml <- glm(Ymat ~ x, family = binomial, data = kalythos)

既然 logit 的關(guān)聯(lián)函數(shù)是默認(rèn)的，因此我們可以在第二條命令中省略該參數(shù)。 為了查看擬合結(jié)果，我們使用

     > summary(fmp)
     > summary(fml)

兩種模型都擬合的很好。為了計(jì)算 LD50，我們可以 利用一個簡單的函數(shù)：

     > ld50 <- function(b) -b[1]/b[2]
     > ldp <- ld50(coef(fmp)); ldl <- ld50(coef(fml)); c(ldp, ldl)

從這些數(shù)據(jù)中得到的年齡分別是43.663年和 43.601年。

Poisson 模型

在 Poisson 族中，默認(rèn)的關(guān)聯(lián)函數(shù)是log。在實(shí)際操作中， 這一族常常用于擬合計(jì)數(shù)資料的 Poisson 對數(shù)線性模型。 這些計(jì)數(shù)資料的實(shí)際分布往往符合二項(xiàng)式分布。 這是一個非常重要而又龐大的話題，我們不想在這里深入展開。 它構(gòu)成了非-gaussian 廣義模型內(nèi)容 的很大一部分。

有時候，實(shí)踐中產(chǎn)生的 Poisson 數(shù)據(jù)在對數(shù)或者平方根 轉(zhuǎn)化后可當(dāng)作正態(tài)數(shù)據(jù)處理。 作為后者的另一種選擇是，一個 Poisson 廣義線性模型可以通過下面的例子 擬合：

     > fmod <- glm(y ~ A + B + x, family = poisson(link=sqrt),
                   data = worm.counts)

擬似然模型

對于所有的族，響應(yīng)變量的方差依賴于均值并且擁有 作為系數(shù)（multiplier）的尺度參數(shù)。 方差對均值的依賴方式是響應(yīng)分布的一個特性； 例如對于poisson分布 Var(y) = mu。

對于擬似然估計(jì)和推斷，我們不是設(shè)定精確的響應(yīng)分布而是 設(shè)定關(guān)聯(lián)函數(shù)和方差函數(shù)的形式。因?yàn)殛P(guān)聯(lián)函數(shù)和方差函數(shù)都依賴于均值。 既然擬似然估計(jì) 和 gaussian 分布使用的技術(shù)非常相似， 因此這一族順帶提供了一種用非標(biāo)準(zhǔn)關(guān)聯(lián)函數(shù)或者方差函數(shù) 擬合gaussian模型的 方法。

例如，考慮非線性回歸的擬合 y = theta_1 z_1 / (z_2 – theta_2) + e 同樣還可以寫成 y = 1 / (beta_1 x_1 + beta_2 x_2) + e 其中 x_1 = z_2/z_1, x_2 = -1/x_1, beta_1 = 1/theta_1, and beta_2 = theta_2/theta_1。 假如有適合的數(shù)據(jù)框，我們可以如下 進(jìn)行非線性擬合

     > nlfit <- glm(y ~ x1 + x2 - 1,
                    family = quasi(link=inverse, variance=constant),
                    data = biochem)