[R] confint.glm(...) fails for binomial count data format

[Prof Brian Ripley]

For the record, this is because that example has a useless row (row 24 has 
no respondents and so adds nothing).  confint() works if you remove the 
pointless row.

We'll add a precautionary check in due course, but such datasets are 
unsurprisingly rare.

On Sun, 16 Nov 2008, Xiaoxu LI wrote:

> ##Q1. confint.glm(...) fails for an example of HSAUR
>
> data("womensrole", package = "HSAUR");
> ## summary(womensrole);
> womensrole_glm_2 <- glm(fm2, data = womensrole,family = binomial())
> ## summary(womensrole_glm_2);
> confint(womensrole_glm_2);
> ## -------Fail---------
> # Waiting for profiling to be done...
> # Error in if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1") :
> #  missing value where TRUE/FALSE needed
> ###############################
>
> ##Q2. Any quick function to transform a count/weight data.frame  into
> a simple factor data.frame?  Dislike "for" routine.
>
> (womensrole.factor <- womensrole[c(),1:2] )
> k=0;
> for (i in as.integer(rownames(womensrole))){
> 	if (womensrole$agree[i] > 0)
> 		for (j in 1:womensrole$agree[i]){
> 			k=k+1;
> 			womensrole.factor[k,1:2]<-womensrole[i,1:2];
> 			womensrole.factor[k,3]<-TRUE;
> 		}
> 	if (womensrole$disagree[i] > 0)
> 		for (j in 1:womensrole$disagree[i]){
> 			k=k+1;
> 			womensrole.factor[k,1:2]<-womensrole[i,1:2];
> 			womensrole.factor[k,3]<-FALSE;
> 		}
> }
> colnames(womensrole.factor)[3]<-'agree';
> ## summary(womensrole.factor)
> ## sum(womensrole$agree)
> ## sum(womensrole$disagree)
>
> ##Two dataset will report same prediction, Chisq and different sample
> size, residual deviance, ...
>
> fm2 <- cbind(agree,disagree) ~ sex * education;
> womensrole_glm_2 <- glm(fm2, data = womensrole, family = binomial());
> womensrole.factor_glm_2 <- glm(agree~sex*education, data =
> womensrole.factor, family = binomial());
> ## Same prediction
> myplot <- function(role.fitted) {
> 	f <- womensrole$sex == "Female"
> 	plot(womensrole$education, role.fitted, type = "n",
> 	ylab = "Probability of agreeing",
> 	xlab = "Education", ylim = c(0,1))
> 	lines(womensrole$education[!f], role.fitted[!f], lty = 1)
> 	lines(womensrole$education[f], role.fitted[f], lty = 2)
> 	lgtxt <- c("Fitted (Males)", "Fitted (Females)")
> 	legend("topright", lgtxt, lty = 1:2, bty = "n")
> 	y <- womensrole$agree / (womensrole$agree +
> 	womensrole$disagree)
> 	text(womensrole$education, y, ifelse(f, "\\VE", "\\MA"),
> 	family = "HersheySerif", cex = 1.25)
> }
> role.fitted2 <- predict(womensrole_glm_2, type = "response");
> myplot(role.fitted2);
> role.fitted2.factor <-
> predict(womensrole.factor_glm_2,newdata=womensrole[,1:2], type =
> "response");
> f <- womensrole$sex == "Female"
> lines(womensrole$education[!f], role.fitted2.factor[!f], lty = 1,col='red');
> lines(womensrole$education[f], role.fitted2.factor[f], lty = 2,col='red');
> ##  Same Chisq, different sample size and residual deviance, AIC
> anova(womensrole_glm_2,test='Chisq')
> anova(womensrole.factor_glm_2,test='Chisq')
>
> ______________________________________________
> R-help at r-project.org 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.
>

-- 
Brian D. Ripley,                  ripley at 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 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595