{"word_count": 232, "character_count": 3034}
Most of our intuitions about regression models come from the Gaussian standard linear model. One interesting feature is that, when we have a factor explanatory variable, the sum of predictions per class is the sum of observations of the endogeneous variable, per class. To be more specific, consider some factor variable , and a regression modelUse ordinary least squares to fit that modelThen for all &#62; n=200&#62; X1=rep(0:1,each=n/2)&#62; set.seed(1)&#62; X2=runif(2*n)&#62; L=X1-X2&#62; B=data.frame(Y=rnorm(n,L),X1=as.factor(X1),X2=X2)&#62; pd=aggregate(x=B$Y,by=list(B$X1),mean)$x&#62; pd[1] -0.4881735  0.5341301&#62; fit=lm(Y~X1+X2,data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] -0.4881735  0.5341301Actually, this result is still wherever the data come from. For instance, if we consider some Poisson based observations,&#62; L=exp(X1-X2)&#62; B=data.frame(Y=rpois(n,L),X1=as.factor(X1),X2=X2)&#62; pd=aggregate(x=B$Y,by=list(B$X1),mean)$x&#62; pd[1] 0.68 1.80&#62; fit=lm(Y~X1+X2,data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.68 1.80This result is still valid here. But this is not the case with Generalized Linear Models. I mean, not everytime. The thing is that it is stil valid with a logistic regression,&#62; B=data.frame(Y=(rpois(n,L)&#62;1)*1,X1=as.factor(X1),X2=X2)&#62; pd=aggregate(x=B$Y,by=list(B$X1),mean)$x&#62; pd[1] 0.13 0.38&#62; fit=glm(Y~X1+X2,data=B,family=binomial(link="logit"))&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.13 0.38But if we move away from the logistic regression, this result is no longer valid. The probit regression for instance&#62; fit=glm(Y~X1+X2,family=binomial(link="probit"),data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.1301143 0.3799545or a cloglog link function&#62; fit=glm(Y~X1+X2,family=binomial(link="cloglog"),data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.1299614 0.3800653or a log link function&#62; fit=glm(Y~X1+X2,family=binomial(link="log"),data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.1299618 0.3800345Actually, we can try any kind of power link function&#62; fit=glm(Y~X1+X2,family=binomial(link=power(.5)),data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.1301992 0.3800619&#62; fit=glm(Y~X1+X2,family=binomial(link=power(.8)),data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.1306485 0.3798376Even if we will soon see a similar property with the Poisson regression,&#62; fit=glm(Y~X1+X2,family=poisson(link="log"),data=B)&#62; B2=data.frame(x=B$X1,y=predict(fit,type="response"))&#62; aggregate(x=B2$y,by=list(B2$x),mean)$x[1] 0.13 0.38in general, for any 