# Regression on factors

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 $x_1\in\{0,1\}$, and a regression model

$y_i=\beta_0+\beta_1 \boldsymbol{1}(x_1=1)+\beta_2 x_2+\varepsilon_i$

Use ordinary least squares to fit that model

$\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 \boldsymbol{1}(x_1=1)+\widehat{\beta}_2 x_2$

Then for all $x\in\{0,1\}$

$\sum_{i:x_i=x} y_i = \sum_{i:x_i=x} \widehat{y}_i$

> n=200 > X1=rep(0:1,each=n/2) > set.seed(1) > X2=runif(2*n) > L=X1-X2 > B=data.frame(Y=rnorm(n,L),X1=as.factor(X1),X2=X2) > pd=aggregate(x=B$Y,by=list(B$X1),mean)$x > pd [1] -0.4881735 0.5341301 > fit=lm(Y~X1+X2,data=B) > B2=data.frame(x=B$X1,y=predict(fit)) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] -0.4881735 0.5341301 Actually, this result is still wherever the data come from. For instance, if we consider some Poisson based observations, > L=exp(X1-X2) > B=data.frame(Y=rpois(n,L),X1=as.factor(X1),X2=X2) > pd=aggregate(x=B$Y,by=list(B$X1),mean)$x > pd [1] 0.68 1.80 > fit=lm(Y~X1+X2,data=B) > B2=data.frame(x=B$X1,y=predict(fit)) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.68 1.80

This 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,

> B=data.frame(Y=(rpois(n,L)>1)*1,X1=as.factor(X1),X2=X2) > pd=aggregate(x=B$Y,by=list(B$X1),mean)$x > pd [1] 0.13 0.38 > fit=glm(Y~X1+X2,data=B,family=binomial(link="logit")) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.13 0.38 But if we move away from the logistic regression, this result is no longer valid. The probit regression for instance > fit=glm(Y~X1+X2,family=binomial(link="probit"),data=B) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.1301143 0.3799545 or a cloglog link function > fit=glm(Y~X1+X2,family=binomial(link="cloglog"),data=B) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.1299614 0.3800653 or a log link function > fit=glm(Y~X1+X2,family=binomial(link="log"),data=B) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.1299618 0.3800345 Actually, we can try any kind of power link function > fit=glm(Y~X1+X2,family=binomial(link=power(.5)),data=B) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.1301992 0.3800619 > fit=glm(Y~X1+X2,family=binomial(link=power(.8)),data=B) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] 0.1306485 0.3798376 Even if we will soon see a similar property with the Poisson regression, > fit=glm(Y~X1+X2,family=poisson(link="log"),data=B) > B2=data.frame(x=B$X1,y=predict(fit,type="response")) > aggregate(x=B2$y,by=list(B2$x),mean)\$x [1] 0.13 0.38

in general, for any $x\in\{0,1\}$

$\sum_{i:x_i=x} y_i \neq \sum_{i:x_i=x} \widehat{y}_i$