VGAM and vglm

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

VGAM and vglm

Ted.Harding-2
Hi Folks,
I wonderif someone who is familiar with the details
of vglm in the VGAM package can assist me. I'm new
to using it, and there doesn;t seem much in the
documentation that's relevant to the question below.

Say I have a vector x of 0/1 responses and another
vector y of 0/1 responses, these in fact being a
bivariate set of 0/1 responses equivalent to
cbind(x,y).

E.g.
set.seed(12345)
x<-sample(c(0,1),50,replace=TRUE)
y<-sample(c(0,1),50,replace=TRUE)
table(x,y)
##    y
##    0  1
##x 0 10 17
##  1 11 12

Now, it seems to me that if I use vglm to fit cbind(x,y)
it simply produces separate fits for x and y, without
taking account of their association. Thus, first (some
lines of the output removed):

summary(vglm(cbind(x,y)~1,family=binomialff(mv=TRUE)))
##Coefficients:
##                 Value Std. Error  t value
##(Intercept):1 -0.16034    0.28375 -0.56508
##(Intercept):2  0.32277    0.28652  1.12651
##
##Number of linear predictors:  2
##Names of linear predictors: logit(E[x]), logit(E[y])

Now if I shuffle say y:

y1<-sample(y,50)
table(x,y1)
##    y1
##    0  1
##x 0 11 16
##  1 10 13

summary(vglm(cbind(x,y1)~1,family=binomialff(mv=TRUE)))
##Coefficients:
##                 Value Std. Error  t value
##(Intercept):1 -0.16034    0.28375 -0.56508
##(Intercept):2  0.32277    0.28652  1.12651
##
##Number of linear predictors:  2
##Names of linear predictors: logit(E[x]), logit(E[y1])

Thus getting exactly the same coefficients; indeed the
results are exactly what I would get if I did separate
glm fits in the first place:

summary(glm(x~1,family=binomial))
##Call:
##glm(formula = x ~ 1, family = binomial)
##Coefficients:
##            Estimate Std. Error z value Pr(>|z|)
##(Intercept)  -0.1603     0.2838  -0.565    0.572

##summary(glm(y~1,family=binomial))
##Coefficients:
##            Estimate Std. Error z value Pr(>|z|)
##(Intercept)   0.3228     0.2865   1.126     0.26

Plus I get the P-values thrown in for good measure!

I had been hoping that vglm, with the "mv=TRUE" option
set in the family=binomialff(mv=TRUE) option, would
treat cbind(x,y) as multivariate data and take account
of their association in preforming the fit (i.e., in
the above example where there are no covariates, would
fit the 2x2 table). This issue becomes more interesting
when there are covariates as well. In fact, it seems to
simply fit the marginals.

So I'm wondering what is to be gained by using vglm
instead of simple glm?

Or is there something I've overlooked, which would
cause a true multivaruate fit to to be done by vglm?

If course, I can use brute force to fit the 2x2 table:

summary(vglm(cbind(x*y,x*(1-y),(1-x)*y,(1-x)*(1-y))~1,
        family=binomialff(mv=TRUE)))
##Coefficients:
##                Value Std. Error t value
##(Intercept):1 -1.1527    0.33113 -3.4810
##(Intercept):2 -1.2657    0.34139 -3.7073
##(Intercept):3 -0.6633    0.29854 -2.2218
##(Intercept):4 -1.3863    0.35355 -3.9210
##
##Number of linear predictors:  4
##Names of linear predictors: logit(mu1),
##  logit(mu2), logit(mu3), logit(mu4)

but that seems to be going to greater lengths than
one should need to!

Any comments and suggestions will be welcome!

Best wishes to all,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <[hidden email]>
Fax-to-email: +44 (0)870 094 0861
Date: 29-Oct-07                                       Time: 23:42:08
------------------------------ XFMail ------------------------------

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: VGAM and vglm

Thomas Yee
Hi Ted,

yes, vglm(..., family=binomialff(mv=TRUE)) does treat each response as
independent. You can see that with coef(fit, matrix=TRUE) where fit is
the object.  If you want to model two dependent binary responses then
you can try the VGAM family functions binom2.or(), loglinb2(),
binom2.rho().

 > So I'm wondering what is to be gained by using vglm instead of
 > simple glm?

Some VGAM family functions handle a multivariate response.  Whether
they are treated as independent or dependent depends on the specific
family function. Handling multivariate responses becomes important
when the modelling function is rrvglm(), cqo(), cao() etc.  Being able
to handle a multivariate response is a good thing, e.g., one can
easily plot several response curves simultaneously, as in the last
example in the vgam() help file.

 > I wonderif someone who is familiar with the details of vglm in the
 > VGAM package can assist me. I'm new to using it, and there doesn;t
 > seem much in the documentation that's relevant to the question
 > below.

Yes, the documentation
(http://www.stat.auckland.ac.nz/~yee/VGAM/download.shtml) is rather
patchy. Sorry! Things are developing there, albeit at a slow pace.

cheers

Thomas

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.