用R語言進行數據分析：常規和廣義線性模型

線性模型

對于常規的多重模型（multiple model）擬合，最基本的函數是lm()。 下面是調用它的方式的一種改進版：

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

例如

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

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

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

廣義線性模型

廣義線性建模是線性模型在研究響應值的非正態分布以及非線性模型的簡潔直接的線性轉化 時的一種發展。 廣義線性模型 是基于下面一系列 假設前提的：

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

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

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

• y 分布的形式為

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

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

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

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

  該可逆函數 ell() 稱為 關聯函數（link function）。

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

族

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

每一種響應分布允許各種關聯函數將均值和線性預測器關聯起來。 這些自動關聯函數如 下表所示：

這些用于模型構建過程中的響應分布，關聯函數和各種 其他必要的信息統稱為 廣義線性模型的族（family）。

glm()函數

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

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

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

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

這些名字是標準的。程序給定的族生成器可以參見 Families 列表中 的“族名”。當選擇一個關聯函數時， 該關聯函數名和族名可以同時在括弧里面作為 參數設定。在擬（quasi） 家族里面，方差函數也是以這種方式設定。

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

gaussian族

命令

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

和下面的命令結果一致。

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

但是效率上，前者差一點。注意，gaussian 族沒有相關參數， 因此它不提供關聯函數的。 如一個問題需要用非標準關聯函數的 gaussian 族， 那么只能采用我們后面討論的擬族。

二項式族

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

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

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

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

第一步是把數據轉換成數據框。

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

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

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

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

     > 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 的關聯函數是默認的，因此我們可以在第二條命令中省略該參數。 為了查看擬合結果，我們使用

     > summary(fmp)
     > summary(fml)

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

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

從這些數據中得到的年齡分別是43.663年和 43.601年。

Poisson 模型

在 Poisson 族中，默認的關聯函數是log。在實際操作中， 這一族常常用于擬合計數資料的 Poisson 對數線性模型。 這些計數資料的實際分布往往符合二項式分布。 這是一個非常重要而又龐大的話題，我們不想在這里深入展開。 它構成了非-gaussian 廣義模型內容 的很大一部分。

有時候，實踐中產生的 Poisson 數據在對數或者平方根 轉化后可當作正態數據處理。 作為后者的另一種選擇是，一個 Poisson 廣義線性模型可以通過下面的例子 擬合：

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

擬似然模型

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

對于擬似然估計和推斷，我們不是設定精確的響應分布而是 設定關聯函數和方差函數的形式。因為關聯函數和方差函數都依賴于均值。 既然擬似然估計 和 gaussian 分布使用的技術非常相似， 因此這一族順帶提供了一種用非標準關聯函數或者方差函數 擬合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。 假如有適合的數據框，我們可以如下 進行非線性擬合

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