2011-08-09

順序ありカテゴリカルデータをvglm()で

本日の1つ前の記事の続き
データを適当に作ってvglm()にかけてみる

# VGAMを使ってみよう
# 最終形質が順序ありカテゴリの場合
# y 量的形質
# y.2が直接観測できずに、順序ありカテゴリ値が得られる
# たとえば {0,1,2,3} 重症度分類
# 説明変数が複数ある
# x1,x2,x3,x4
# x1 性別
# x2 発症時年齢
# x3 発症時検査項目t1 (量的因子)
# x4 発症時検査項目t2 (順序なしカテゴリ)
# サンプル数 N

N<-1000
f.frac<-0.8
Nf<-N*f.frac
Nm<-N-Nf
x1<-c(rep(0,Nf),rep(1,Nm))
x2<-c(rnorm(Nf,mean=30,sd=5),rnorm(Nm,mean=60,sd=10))
hist(x2)

x3<-rnorm(N,10)*x2+rnorm(N,10,sd=5)*mean(x2)
hist(x3)
x4<-sample(c("A","B","C","D","E"),N,replace=TRUE)
# yがx1,x2,x3に依存するようにyを計算する
y<-x1*10+x2*x3/1000 + rnorm(N,10,20)
# x4は順序なしカテゴリ。カテゴリ"C"だけ、寄与させる
y[which(x4=="C")]<-y[which(x4=="C")]*2
y.2<-rep(0,length(y))

qy<-quantile(y)
y2<-which(y>qy[2])
y3<-which(y>qy[3])
y4<-which(y>qy[4])

y.2[y2]<-1
y.2[y3]<-2
y.2[y4]<-3


plot(x3,y.2)

# 説明変数なしのモデル
fit.out0 <- vglm(y.2 ~ 1,family = cumulative(parallel = TRUE, reverse = TRUE))
head(coef(fit.out0,matrix=TRUE),10)

head(coef(summary(fit.out0)),10)

# すべての変数x1,x2,x3,x4にリニアな係数(parallel=TRUE)を設定したモデル
fit.out1 <- vglm(y.2 ~ x1+x2+x3+x4,family = cumulative(parallel = TRUE, reverse = TRUE))
# 説明変数は、カテゴリの数の列に分けて格納されている
colSums(fit.out1@y)
# 推定された係数が、y.2の段階によらず同じである(parallel=TRUE)であることに注意
head(coef(fit.out1,matrix=TRUE),10)

head(coef(summary(fit.out1)),10)

# 検定は「モデルなし」と「モデルあり」との間で尤度比検定
# 尤度比検定は、カイ自乗分布を仮定した統計量と自由度で決まる(自由度は増やした変数の個数)
pchisq(deviance(fit.out0)-deviance(fit.out1),df=df.residual(fit.out0)-df.residual(fit.out1),lower.tail=FALSE)

# x4だけ、paralleleの条件を外す(parallel=FALSE~1+x4)と、変数が増えるから自由度も増える
# そのような変更が有意かどうかも尤度比検定でみてやる

fit.out2 <- vglm(y.2 ~ x1+x2+x3+x4,family = cumulative(parallel = FALSE ~ 1+x4, reverse = TRUE))
head(coef(fit.out2,matrix=TRUE),10)

head(coef(summary(fit.out2)),10)
pchisq(deviance(fit.out1)-deviance(fit.out2),df=df.residual(fit.out1)-df.residual(fit.out2),lower.tail=FALSE)

出力

> # 説明変数なしのモデル
> fit.out0 <- vglm(y.2 ~ 1,family = cumulative(parallel = TRUE, reverse = TRUE))
> head(coef(fit.out0,matrix=TRUE),10)
            logit(P[Y>=2]) logit(P[Y>=3]) logit(P[Y>=4])
(Intercept)       1.098612   1.274734e-16      -1.098612
> 
> head(coef(summary(fit.out0)),10)
                      Value Std. Error       t value
(Intercept):1  1.098612e+00 0.07302967  1.504337e+01
(Intercept):2  1.274734e-16 0.06324555  2.015531e-15
(Intercept):3 -1.098612e+00 0.07302967 -1.504337e+01
> 
> # すべての変数x1,x2,x3,x4にリニアな係数(parallel=TRUE)を設定したモデル
> fit.out1 <- vglm(y.2 ~ x1+x2+x3+x4,family = cumulative(parallel = TRUE, reverse = TRUE))
 警告メッセージ： 
In model.matrix.default(mt, mf, contrasts) :
  variable 'x4' converted to a factor
> # 説明変数は、カテゴリの数の列に分けて格納されている
> colSums(fit.out1@y)
  0   1   2   3 
250 250 250 250 
> # 推定された係数が、y.2の段階によらず同じである(parallel=TRUE)であることに注意
> head(coef(fit.out1,matrix=TRUE),10)
            logit(P[Y>=2]) logit(P[Y>=3]) logit(P[Y>=4])
(Intercept)    -3.37798590    -4.89396871    -6.85838071
x1              0.80068004     0.80068004     0.80068004
x2              0.06732731     0.06732731     0.06732731
x3              0.00282544     0.00282544     0.00282544
x4B             0.17519237     0.17519237     0.17519237
x4C             2.51468084     2.51468084     2.51468084
x4D            -0.02673022    -0.02673022    -0.02673022
x4E             0.07796965     0.07796965     0.07796965
> 
> head(coef(summary(fit.out1)),10)
                    Value   Std. Error     t value
(Intercept):1 -3.37798590 0.4013065486  -8.4174702
(Intercept):2 -4.89396871 0.4156017513 -11.7756210
(Intercept):3 -6.85838071 0.4456535110 -15.3894910
x1             0.80068004 0.3462284361   2.3125773
x2             0.06732731 0.0114422468   5.8840985
x3             0.00282544 0.0003674152   7.6900467
x4B            0.17519237 0.1943464748   0.9014435
x4C            2.51468084 0.2117535184  11.8755091
x4D           -0.02673022 0.1949736865  -0.1370966
x4E            0.07796965 0.1984772896   0.3928392
> 
> # 検定は「モデルなし」と「モデルあり」との間で尤度比検定
> # 尤度比検定は、カイ自乗分布を仮定した統計量と自由度で決まる(自由度は増やした変数の個数)
> pchisq(deviance(fit.out0)-deviance(fit.out1),df=df.residual(fit.out0)-df.residual(fit.out1),lower.tail=FALSE)
[1] 8.493108e-134
> 
> # x4だけ、paralleleの条件を外す(parallel=FALSE~1+x4)と、変数が増えるから自由度も増える
> # そのような変更が有意かどうかも尤度比検定でみてやる
> 
> fit.out2 <- vglm(y.2 ~ x1+x2+x3+x4,family = cumulative(parallel = FALSE ~ 1+x4, reverse = TRUE))
 警告メッセージ： 
In model.matrix.default(mt, mf, contrasts) :
  variable 'x4' converted to a factor
> head(coef(fit.out2,matrix=TRUE),10)
            logit(P[Y>=2]) logit(P[Y>=3]) logit(P[Y>=4])
(Intercept)   -3.468009222   -5.057181562   -7.770261183
x1             1.083793033    1.083793033    1.083793033
x2             0.072781576    0.072781576    0.072781576
x3             0.002867438    0.002867438    0.002867438
x4B            0.145736083    0.104751298    0.454973687
x4C            1.200523387    1.974227617    3.579792821
x4D            0.017035933   -0.259397806    0.243186015
x4E            0.038567630    0.172056690   -0.045528562
> 
> head(coef(summary(fit.out2)),10)
                     Value   Std. Error     t value
(Intercept):1 -3.468009222 0.4182957827  -8.2908061
(Intercept):2 -5.057181562 0.4327272985 -11.6867634
(Intercept):3 -7.770261183 0.5281333590 -14.7126877
x1             1.083793033 0.3663289848   2.9585238
x2             0.072781576 0.0115986302   6.2750148
x3             0.002867438 0.0003716743   7.7149203
x4B:1          0.145736083 0.2432247429   0.5991828
x4B:2          0.104751298 0.2419113476   0.4330152
x4B:3          0.454973687 0.4021672687   1.1313046
x4C:1          1.200523387 0.2759482980   4.3505374
> pchisq(deviance(fit.out1)-deviance(fit.out2),df=df.residual(fit.out1)-df.residual(fit.out2),lower.tail=FALSE)
[1] 3.443547e-10