##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')
Prof Brian Ripley
2008-Dec-07 16:28 UTC
[R] confint.glm(...) fails for binomial count data format
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