一般化線形モデル（Generalized Linear Model: GLM）のメモ。

【目次】

・線形回帰モデルでは、誤差が互いに独立に正規分布に従うことを前提としている（不遍性・等分散性・独立性・正規性）。

$\qquad Y=X' \beta + \varepsilon , \quad \varepsilon \sim N( 0,\sigma ^{2})$

・そのため誤差構造が正規分布ではない場合、モデルの当てはまりが悪い。
・GLMでは誤差構造が何の分布に従うかを指定できる。

計算式

GLMは分布族・リンク関数・線形予測子でモデリングする。
「リンク関数 (期待値) ＝線形予測子」で期待値を $(-∞, ∞)$ のスケールに写像し、この空間で線形モデルを作成する。
「逆リンク関数（リンク関数の逆数）」で線形予測子を元のスケールに写像し、期待値を得る。

$\quad g(θ)=X' \beta$

$\quad E(Y|X) = θ = g^{-1}(X' \beta )$

■ ロジスティック回帰での例
・ロジスティック回帰では誤差構造が二項分布に従うと仮定する。
・ロジット関数（リンク関数）でロジット値に変換し、線形モデルを作成する。確率値は(0, 1)の範囲をとるが、ロジット値は $(-∞, ∞)$ の範囲となる。

$\qquad g(θ)=ln (\dfrac{θ}{1-θ})=X' \beta$

・ロジスティック関数（ロジット関数の逆数、シグモイド関数とも言う）を使って線形予測子から確率値を得る。

$\qquad E(Y|X) = θ = g^{-1}( X' \beta )=\dfrac{exp( X' \beta )}{1+exp( X' \beta )}$

・下図は1%刻みの確率値とそのロジット値の対応図。

（Rコードを見る）

data1 <- data.frame(p=seq(0,1,0.01)) %>% mutate(logit=log(p/(1-p)), type="確率")
data2 <- data1 %>% mutate(type="ロジット")
data3 <- union_all(data1,data2) %>% 
  mutate( val =if_else(type=="確率", p, logit)
         ,grp =p
         ,type=factor(type, level=c("確率","ロジット")) )

ggplot(data=data3, mapping=aes(type, val, group=grp)) + 
  geom_point() +
  geom_line() +
  labs(x="", y="", title="確率値とロジット値", caption="🄫Cochineal19") +
  scale_y_continuous(limits=c(-6,6), breaks=c(seq(-5,5,1),0.5)) +
  theme(plot.title=element_text(size=30)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20)
        ,legend.title=element_text(size=20)
        ,legend.text=element_text(size=20))

GLMの例

Rのfamilyオブジェクトの分布族とリンク関数（デフォルト）：

binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")

上記から４モデルを記載。

モデル	目的変数	分布族	リンク関数 $g(θ)$	期待値 $E(Y\|X)=θ$
線形回帰	連続値 $(-∞ ～ ∞)$	正規分布	恒等リンク $θ=X' \beta$	$θ=X' \beta$
ロジスティック回帰	二値の離散値 $(0, 1)$	二項分布	ロジットリンク $ln (\dfrac{θ}{1-θ})=X' \beta$	$θ=\dfrac{exp(X' \beta )}{1+exp(X' \beta )}$
ポアソン回帰	非負の離散値 $(0 ～ ∞)$	ポアソン分布	対数リンク $ln (θ)=X' \beta$	$θ=exp(X' \beta )$
ガンマ回帰	非負の連続値 $(0 ～ ∞)$	ガンマ分布	逆数リンク $-θ^{-1}=X' \beta$	$θ=-(X' \beta )^{-1}$

glm(y~x, data=ads, family=gaussian(link="identity"))  #-- 線形回帰
glm(y~x, data=ads, family=binomial(link="logit"))  #-- ロジスティック回帰
glm(y~x, data=ads, family=poisson(link="log"))  #-- ポアソン回帰
glm(y~x, data=ads, family=Gamma(link="inverse"))  #-- ガンマ回帰
  
geeglm(y~x, data=ads,family=poisson(link="log"), id=id)  #-- 修正ポアソン回帰

#-- ロジスティック回帰の例
proc genmod data=ads;
    class effects(ref="control") ;
    model response(event="1") = effects / dist=binomial link=logit lrci type3 ;
    ods output ParameterEstimates=_Est;
run;

SAS Help Center

ちなみに
ロジスティック回帰の回帰係数を指数変換⇒オッズ比
ポアソン回帰の回帰係数を指数関数⇒リスク比
（コックス回帰の回帰係数を指数変換⇒ハザード比）

> library(tidyverse)
> library(geepack)
> n <- 100
> ads <- data.frame(
+    id=c(seq(1,n*2))
+   ,x=c(rep(0,n),rep(1,n))
+   ,y=c(rbinom(n,1,0.7),rbinom(n,1,0.3)) )
> (mtx1 <- matrix(xtabs(~x+y,data=ads),2))
     [,1] [,2]
[1,]   22   78
[2,]   67   33
> glm1 <- glm(y~x, data=ads, family=binomial(link="logit"))  #-- ロジスティック回帰
> glm2 <- glm(y~x, data=ads, family=poisson(link="log"))  #-- ポアソン回帰
> exp(coef(glm1))  #-- オッズ比
(Intercept)           x 
  3.5454545   0.1389208 
> exp(coef(glm2))  #-- リスク比
(Intercept)           x 
  0.7800000   0.4230769 
> 
> glm3 <- geeglm(y~x, data=ads,family=poisson(link="log"), id=id)  #-- 修正ポアソン回帰
> exp(coef(glm3))  #-- リスク比
(Intercept)           x 
  0.7800000   0.4230769