##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@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.
