Re: [S] 'glim' and 'glm' covariances

Prof Brian Ripley (ripley@stats.ox.ac.uk)
Tue, 14 Jul 1998 17:23:55 +0100 (BST)


Your code does not even work for me. Try the following recommended
approaches:

> glim(X, y = r, error = "binomial", link = "probit", n = n, intercept=F)
$coef:
X1 X2 X3 X4
-4.221782 -2.004999 -3.694616 1.937577

$var:
[,1] [,2] [,3] [,4]
[1,] 0.4392274 0.19038032 0.3478578 -0.17621315
[2,] 0.1903803 0.16359623 0.1804440 -0.09140688
[3,] 0.3478578 0.18044397 0.3950234 -0.16701621
[4,] -0.1762132 -0.09140688 -0.1670162 0.08460483

$deviance:
[1] 104.609067 3.327818

$df:
[1] 12 8

$scale:
[1] 1

$intercept:
[1] F

or, better,

> fm <- glm(cbind(r, n-r) ~ X -1, binomial(link = probit))
> library(MASS)
> vcov(fm)
X1 X2 X3 X4
X1 0.4390593 0.19030009 0.3477222 -0.17614027
X2 0.1903001 0.16354522 0.1803751 -0.09136983
X3 0.3477222 0.18037510 0.3948984 -0.16695378
X4 -0.1761403 -0.09136983 -0.1669538 0.08457119

[Note, the right formula is to multiply cov.unscaled by the dispersion,
here 1. That is what the glm method for vcov does.]

In this example the deviance, 3.32 on 8 df, is rather small, so so
people would suspect underdispersion and reduce the estimate of the
dispersion.

> Date: Tue, 14 Jul 1998 16:51:16 +0200
> From: Manfred Wilhelm <wilma@pei.de>
> To: s-news@wubios.wustl.edu
> Subject: [S] 'glim' and 'glm' covariances
>
> Hallo,
>
> in order to compare three preparations of a quantal bioassy I perform a probit
analysis. The matrix X below is designed to yield three different intercepts and
a common slope (X[,4] = dilution).
> To obtain Fieller's confidence limits for the relative potency I need the
covariance parameters of the estimated coefficients.
>
> I used two methods,'glm' and 'glim' in S-PLUS 4.5 for Windows, but got
extremley different covariance matrices, while only moderate differences in the
coefficients occur. In both cases the dispersion or scale parameter is taken to
be 1 (as usual for binomial error or family). Which one is the correct way?
> (For perfect confusion: the data are from a worked example in the literature
resulting in a third covariance matrix lying in between those of 'glim' and
'glm'.)
>
> Differences seem to be due to a common mulitplicity constant for all
covariance parameters. I also recognized that covariance estimates differ by
changing the 'wt' or 'weights' arguments. Are the chosen ones not sensible? Do I
have to use any scale or dispersion parameter and which one?
>
> > X <-
cbind(c(rep(1,4),rep(0,8)),c(rep(0,4),rep(1,4),rep(0,4)),c(rep(0,8),rep(1,4)),lo
g(c(15,7.5,3.75,1.875,10,5,2.5,1.25,32,16,8,4)))
>
> > X
> [,1] [,2] [,3] [,4]
> [1,] 1 0 0 2.7080502
> [2,] 1 0 0 2.0149030
> [3,] 1 0 0 1.3217558
> [4,] 1 0 0 0.6286087
> [5,] 0 1 0 2.3025851
> [6,] 0 1 0 1.6094379
> [7,] 0 1 0 0.9162907
> [8,] 0 1 0 0.2231436
> [9,] 0 0 1 3.4657359
> [10,] 0 0 1 2.7725887
> [11,] 0 0 1 2.0794415
> [12,] 0 0 1 1.3862944
>
> > r <- c(9,3,1,0,12,10,6,0,12,11,8,2) # no. of reactive animals
> > n <- c(10,11,12,11,12,12,12,10,12,12,12,12) # no. of challenged animals
> > p <- r/n
>
>
> > model.glim_glim(X, y = r, error = "binomial", link = "probit", n = n, wt =
n, intercept = F)
>
> > model.glim$coef
> X1 X2 X3 X4
> -4.167826 -1.961318 -3.641401 1.910795
>
> > model.glim$var
> [,1] [,2] [,3] [,4]
> [1,] 0.03698645 0.016329518 0.02935748 -0.014847039
> [2,] 0.01632952 0.014264447 0.01576007 -0.007970386
> [3,] 0.02935748 0.015760073 0.03372016 -0.014329291
> [4,] -0.01484704 -0.007970386 -0.01432929 0.007246792
>
>
> > model.glm_glm(p ~ -1 + X, family = binomial(link = probit), weights = n)
>
> > model.glm$coef
> X1 X2 X3 X4
> -4.221851 -2.005094 -3.694732 1.937629
>
> > summary(model.glm)$cov.unscaled
> X1 X2 X3 X4
> X1 0.4397946 0.1906967 0.3483506 -0.17647000
> X2 0.1906967 0.1638209 0.1807388 -0.09156000
> X3 0.3483506 0.1807388 0.3955179 -0.16725506
> X4 -0.1764700 -0.0915600 -0.1672551 0.08472929
>
>
> Thanks for any help,
>
> Manfred Wilhelm
>
>
>
> _______________________________________
>
> Dr. Manfred Wilhelm
> Paul-Ehrlich-Institut (PEI)
> Federal Agency for Sera and Vaccines
> Paul-Ehrlich-Str. 51-59
> D-63225 Langen
> Tel: +49 6103 77 2064
> Fax: +49 6103 77 1253
> Email: wilma@pei.de
> Web: www.pei.de
> _______________________________________
>
>
>
> -----------------------------------------------------------------------
> This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
> send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
> message: unsubscribe s-news

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

----------------------------------------------------------------------- This message was distributed by s-news@wubios.wustl.edu. To unsubscribe send e-mail to s-news-request@wubios.wustl.edu with the BODY of the message: unsubscribe s-news