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

Xiaoxu LI lixiaoxu at gmail.com
Sun Nov 16 09:09:05 CET 2008


##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')



More information about the R-help mailing list