[R] constraint matrices in vglm (VGAM package)
Thomas Yee
t.yee at auckland.ac.nz
Sun Nov 18 21:15:24 CET 2007
Hi Lin,
try
multi2 = vglm(case123con ~ SNP_A1+SNP_A2+age,
multinomial(parallel = TRUE ~ SNP_A1+SNP_A2 - 1),
work.analy)
or
multi3 = vglm(case123con ~ SNP_A1+SNP_A2+age,
multinomial(parallel = FALSE ~ 1 + age),
work.analy)
After your fit, typing coef(multi2, mat=TRUE) is a way of checking.
Also constraints(multi2).
Further info is at
http://www.stat.auckland.ac.nz/~yee/VGAM/doc/constraints.pdf
More generally, one can use the constraints argument to pass in a
list with explicit 'constraint matrices'; the parallel argument
here is just a more convenient way to set them up.
cheers
Thomas
> Hello R users,
>
> I am performing a multinomial logit regression and would like to
> constrain a few model coefficients to be equal.
>
> Here is my model:
>
> multi <- vglm(case123con ~ SNP_A1+SNP_A2+age, multinomial, work.analy)
>
> where case123con is a four level categorical variable (case 1, case 2,
> case 3, control) and SNP_A1 and SNP_A2 are indicator functions (yes/no).
>
> The output of this model is:
>
> Coefficients:
>
> Value Std. Error t value
>
> (Intercept):1 6.9798044 0.8145521 8.5688866
> (Intercept):2 2.5729346 0.8733182 2.9461592
> (Intercept):3 0.0129745 1.0966651 0.0118308
> age:1 -0.0468339 0.0087766 -5.3362053
> age:2 -0.0306066 0.0091644 -3.3397314
> age:3 -0.0058718 0.0114316 -0.5136443
> SNP_A1:1 0.0284303 0.1590333 0.1787698
> SNP_A1:2 0.1685160 0.1674925 1.0061104
> SNP_A1:3 0.0137997 0.2052559 0.0672315
> SNP_A2:1 -0.6265717 0.2694583 -2.3253010
> SNP_A2:2 -0.2528873 0.2790632 -0.9062006
> SNP_A2:3 -0.0789181 0.3339153 -0.2363418
>
> I would like to constrain the model coefficients so that SNP_A1:1 =
> SNP_A1:2 = SNP_A1:3 and SNP_A2:1 = SNP_A2:2 = SNP_A2:3. The reason I
> want to do this is because I want to test the null hypothesis that these
> coefficients are equal through a likelihood ratio test. Can I do this
> by defining the constraint matrices and if so, how do I go about doing
> so?
>
> Any help or tips would be greatly appreciated! Thank you in advance!!
>
> Lin
>
More information about the R-help
mailing list