2011-08-09

順序なしカテゴリカルデータをvglm()で

前のデータで順序なしカテゴリであるx4を被説明変数にしてy.2を順序ありカテゴリの説明変数にする

# 最終形質が順序なしカテゴリの場合

# 順序なしカテゴリ x4を被説明変数にしてやる

fat.out0 <- vgam(x4 ~ 1,family = cumulative(parallel = TRUE, reverse = TRUE))
head(coef(fat.out0,matrix=TRUE),10)

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


fat.out1 <- vglm(x4 ~ x1+x2+x3+y.2,family = cumulative(parallel = TRUE, reverse = TRUE))

colSums(fat.out1@y)

head(coef(fat.out1,matrix=TRUE),10)

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


pchisq(deviance(fat.out0)-deviance(fat.out1),df=df.residual(fat.out0)-df.residual(fat.out1),lower.tail=FALSE)


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

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

出力

> # 順序なしカテゴリ x4を被説明変数にしてやる
> 
> fat.out0 <- vgam(x4 ~ 1,family = cumulative(parallel = TRUE, reverse = TRUE))
 警告メッセージ： 
In model.matrix.default(mt, mf, contrasts) :
  variable 'x4' converted to a factor
> head(coef(fat.out0,matrix=TRUE),10)
            logit(P[Y>=2]) logit(P[Y>=3]) logit(P[Y>=4]) logit(P[Y>=5])
(Intercept)       1.463058      0.4599402     -0.4599402      -1.482832
> 
> head(coef(summary(fat.out0)),10)
(Intercept):1 (Intercept):2 (Intercept):3 (Intercept):4 
    1.4630584     0.4599402    -0.4599402    -1.4828323 
> 
> 
> fat.out1 <- vglm(x4 ~ x1+x2+x3+y.2,family = cumulative(parallel = TRUE, reverse = TRUE))
 警告メッセージ： 
In model.matrix.default(mt, mf, contrasts) :
  variable 'x4' converted to a factor
> 
> colSums(fat.out1@y)
  A   B   C   D   E 
188 199 226 202 185 
> 
> head(coef(fat.out1,matrix=TRUE),10)
            logit(P[Y>=2]) logit(P[Y>=3]) logit(P[Y>=4]) logit(P[Y>=5])
(Intercept)   1.8816656095   0.8754619270  -0.0475793864  -1.0724638363
x1            0.4572981994   0.4572981994   0.4572981994   0.4572981994
x2           -0.0064020316  -0.0064020316  -0.0064020316  -0.0064020316
x3           -0.0003846138  -0.0003846138  -0.0003846138  -0.0003846138
y.2           0.0009281906   0.0009281906   0.0009281906   0.0009281906
> 
> head(coef(summary(fat.out1)),10)
                      Value  Std. Error     t value
(Intercept):1  1.8816656095 0.307879495  6.11169513
(Intercept):2  0.8754619270 0.302953928  2.88975268
(Intercept):3 -0.0475793864 0.301653599 -0.15772856
(Intercept):4 -1.0724638363 0.304455951 -3.52255829
x1             0.4572981994 0.300156547  1.52353231
x2            -0.0064020316 0.009480932 -0.67525340
x3            -0.0003846138 0.000325561 -1.18138791
y.2            0.0009281906 0.061208467  0.01516441
> 
> 
> pchisq(deviance(fat.out0)-deviance(fat.out1),df=df.residual(fat.out0)-df.residual(fat.out1),lower.tail=FALSE)
[1] 0.3965315
> 
> 
> fat.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(fat.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(fat.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(fat.out1)-deviance(fat.out2),df=df.residual(fat.out1)-df.residual(fat.out2),lower.tail=FALSE)
[1] 0.0031684
>