From cberry at tajo.ucsd.edu  Sat May  3 20:44:45 2008
From: cberry at tajo.ucsd.edu (Charles C. Berry)
Date: Sat, 3 May 2008 11:44:45 -0700
Subject: [Rd] zero variance in part of a glm (PR#11355)
In-Reply-To: <20080501072509.561C22834169@mail.pubhealth.ku.dk>
References: <20080501072509.561C22834169@mail.pubhealth.ku.dk>
Message-ID: <Pine.LNX.4.64.0805031124490.24535@tajo.ucsd.edu>



NOT a bug.

The estimate of the probability of mi is precisely 1.0 for one cell in 
your setup.

The information contributed by that cell has a factor of 
p*(1-p) in it, which works out to zero.

The bad behavior of asymptotic methods for binary regressions in such 
settings is well known.

HTH,

Chuck

On Thu, 1 May 2008, m.crawley at imperial.ac.uk wrote:

> In this real example (below), all four of the replicates in one
> treatment combination had zero failures, and this produced a very high
> standard error in the summary.lm.
> =20
> Just adding one failure to one of the replicates produced a well-behaved
> standard error.
> =20
> I don't know if this is a bug, but it is certainly hard for users to
> understand.
> =20
> I would value your comments=20
> =20
> Thanks
> =20
> Best wishes,
> =20
> Mick
> =20
> Prof  M.J. Crawley
> =20
> Imperial College London
> Silwood Park
> Ascot
> Berks
> SL5 7PY
> UK
> =20
> Phone (0) 207 5942 216
> Fax     (0) 207 5942 339
> =20
> The data are from a bioassay in which a factorial experiment with two
> factors (antibio and toxin) each with three levels was replicated four
> times. The response is "mi" and "n-mi"
> =20
> Note that lines 17 to 20 in the dataframe have no failures (24 dead out
> of 24 individuals)
>
> data<-read.table("c:\\temp\\ab1.txt",header=3DT)
> attach(data)
> names(data)
> =20
> [1] "antibio" "toxin"   "rep"     "alive"   "af"      "dead"    "mi"
>
> [8] "n"=20=20=20=20
> =20
> data
> =20
> =20
> =20
>   antibio   toxin rep alive af dead mi  n
> 1     camp control   1    24  0    0  0 24
> 2     camp control   2    23  0    1  1 24
> 3     camp control   3    23  0    1  1 24
> 4     camp control   4    21  0    3  3 24
> 5     camp  Cry1Ab   1    21  4    3  7 24
> 6     camp  Cry1Ab   2    20  3    4  7 24
> 7     camp  Cry1Ab   3    20  4    4  8 24
> 8     camp  Cry1Ab   4    18  7    6 13 24
> 9     camp   Vip3A   1     7  3   17 20 24
> 10    camp   Vip3A   2    10  6   14 20 24
> 11    camp   Vip3A   3    11  5   13 18 24
> 12    camp   Vip3A   4    10  2   14 16 24
> 13   bcock control   1    21  0    3  3 24
> 14   bcock control   2    24  0    0  0 24
> 15   bcock control   3    24  3    0  3 24
> 16   bcock control   4    23  1    1  2 24
> 17   bcock  Cry1Ab   1     4  4   20 24 24
> 18   bcock  Cry1Ab   2     4  4   20 24 24
> 19   bcock  Cry1Ab   3     1  1   23 24 24
> 20   bcock  Cry1Ab   4     2  2   22 24 24
> 21   bcock   Vip3A   1    11  4   13 17 24
> 22   bcock   Vip3A   2    15  5    9 14 24
> 23   bcock   Vip3A   3     5  3   19 22 24
> 24   bcock   Vip3A   4    10  6   14 20 24
> 25     ana control   1    23  0    1  1 24
> 26     ana control   2    23  0    1  1 24
> 27     ana control   3    22  0    2  2 24
> 28     ana control   4    23  0    1  1 24
> 29     ana  Cry1Ab   1    18  1    6  7 24
> 30     ana  Cry1Ab   2    17  4    7 11 24
> 31     ana  Cry1Ab   3    13  4   11 15 24
> 32     ana  Cry1Ab   4    14  3   10 13 24
> 33     ana   Vip3A   1    21 12    3 15 24
> 34     ana   Vip3A   2    12  6   12 18 24
> 35     ana   Vip3A   3     9  5   15 20 24
> 36     ana   Vip3A   4     9  1   15 16 24
> =20
> y<-cbind(mi,n-mi)
> model<-glm(y~antibio*toxin,quasibinomial)
> summary(model)
> =20
>
> Call:
> glm(formula =3D y ~ antibio * toxin, family =3D quasibinomial)
> =20
> Deviance Residuals:=20
>    Min       1Q   Median       3Q      Max=20=20
> -2.0437  -0.5645  -0.1022   0.6921   1.9996=20=20
> =20
> Coefficients:
>                           Estimate Std. Error  t value Pr(>|t|)=20=20=20=
> =20
> (Intercept)              -2.901e+00  5.020e-01   -5.780 3.79e-06 ***
> antibiobcock              5.035e-01  6.441e-01    0.782    0.441=20=20=20=
> =20
> antibiocamp               2.100e-15  7.099e-01 2.96e-15    1.000=20=20=20=
> =20
> toxinCry1Ab               2.818e+00  5.494e-01    5.129 2.15e-05 ***
> toxinVip3A                3.840e+00  5.600e-01    6.857 2.29e-07 ***
> antibiobcock:toxinCry1Ab  2.050e+01  2.365e+03    0.009    0.993=20=20=20=
> =20
> antibiocamp:toxinCry1Ab  -4.721e-01  7.795e-01   -0.606    0.550=20=20=20=
> =20
> antibiobcock:toxinVip3A  -2.868e-01  7.381e-01   -0.389    0.701=20=20=20=
> =20
> antibiocamp:toxinVip3A    2.748e-01  7.975e-01    0.345    0.733=20=20=20=
> =20
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1=20
> =20
> (Dispersion parameter for quasibinomial family taken to be 1.19442)
> =20
>    Null deviance: 515.115  on 35  degrees of freedom
> Residual deviance:  35.047  on 27  degrees of freedom
> AIC: NA
> =20
> Number of Fisher Scoring iterations: 17
> =20
> # note the standard error of the firast interaction term =3D 2365)
> =20
> # add a single failure to one replicate
> =20
> y2<-y
> y2[17,]<-c(23,1)
> model2<-glm(y2~antibio*toxin,quasibinomial)
> summary(model2)
> =20
> Call:
> glm(formula =3D y2 ~ antibio * toxin, family =3D quasibinomial)
> =20
> Deviance Residuals:=20
>   Min      1Q  Median      3Q     Max=20=20
> -2.044  -0.627  -0.221   0.709   2.000=20=20
> =20
> Coefficients:
>                           Estimate Std. Error   t value Pr(>|t|)=20=20=20=
> =20
> (Intercept)              -2.901e+00  5.251e-01    -5.526 7.44e-06 ***
> antibiobcock              5.035e-01  6.737e-01     0.747    0.461=20=20=20=
> =20
> antibiocamp              -4.058e-18  7.426e-01 -5.47e-18    1.000=20=20=20=
> =20
> toxinCry1Ab               2.818e+00  5.747e-01     4.904 3.94e-05 ***
> toxinVip3A                3.840e+00  5.857e-01     6.556 4.96e-07 ***
> antibiobcock:toxinCry1Ab  4.134e+00  1.352e+00     3.057    0.005 **=20
> antibiocamp:toxinCry1Ab  -4.721e-01  8.153e-01    -0.579    0.567=20=20=20=
> =20
> antibiobcock:toxinVip3A  -2.868e-01  7.720e-01    -0.372    0.713=20=20=20=
> =20
> antibiocamp:toxinVip3A    2.748e-01  8.341e-01     0.329    0.744=20=20=20=
> =20
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1=20
> =20
> (Dispersion parameter for quasibinomial family taken to be 1.306703)
> =20
>    Null deviance: 506.602  on 35  degrees of freedom
> Residual deviance:  37.851  on 27  degrees of freedom
> AIC: NA
> =20
> Number of Fisher Scoring iterations: 5
> =20
> #   Now the standard errors are all well-behaved
> =20
>
> =20
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901



[Rd] zero variance in part of a glm (PR#11355)
Charles C. Berry