# 実例）比率データの一般化線形モデル 

# 集団内の性別の個体数に関する比率データ
numbers <- read.table("../data/sexratio.txt", header=T)
numbers

##   density females males
## 1       1       1     0
## 2       4       3     1
## 3      10       7     3
## 4      22      18     4
## 5      55      22    33
## 6     121      41    80
## 7     210      52   158
## 8     444      79   365

attach(numbers)

# 性比データの描画
par(mfrow=c(1,2))
plot(density, males/density)
plot(log(density), males/density)

par(mfrow=c(1,1))

# ロジット変換
# オスの比率を求める
p <- males/density 

par(mfrow=c(1,2))
# ロジットの描画
plot(density, log(p/(1-p))) 
# ロジットの描画（対数密度）
plot(log(density), log(p/(1-p))) 

par(mfrow=c(1,1))

# 一般化線形モデルのあてはめ
# オスとメスの対のデータを合体
y <- cbind(males, females) 
# 対数密度に対してロジット♂比を一般化線形モデルで説明する
output <- glm(y ~ log(density), family=binomial(logit))
# 解析結果のモデル評価（AIC）
summary(output) 

## 
## Call:
## glm(formula = y ~ log(density), family = binomial(logit))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9697  -0.3411   0.1499   0.4019   1.0372  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.65927    0.48758  -5.454 4.92e-08 ***
## log(density)  0.69410    0.09056   7.665 1.80e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 71.1593  on 7  degrees of freedom
## Residual deviance:  5.6739  on 6  degrees of freedom
## AIC: 38.201
## 
## Number of Fisher Scoring iterations: 4

# 逸脱度分析（analysis of deviance）
library(car)

## Loading required package: carData

# Exponentiated coefficients ("odds ratios")
exp(coef(output))  

##  (Intercept) log(density) 
##   0.06999959   2.00190680

Anova(output, type="II", test="LR")

## Analysis of Deviance Table (Type II tests)
## 
## Response: y
##              LR Chisq Df Pr(>Chisq)    
## log(density)   65.485  1  5.855e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 最適モデルの描画
# 区間0 ≦ x ≦ 6 の数列を入力
xx <- seq(0, 6, 0.1)
# オスの比率を対数密度に対してプロットする
plot(log(density), p, ylab="Proportion of Males")
# オスの予測ロジット比率を元の比率に逆変換して描画する
lines(xx,predict(output, list(density=exp(xx)), type="response"))

detach(numbers)