Tag Archives: categorical

Regression

Combining automatically factor levels with trees

Leave a comment

Last year, in a post, I discussed how to merge levels of factor variables, using combinatorial techniques (it was for my STT5100 cours, and trees are not in the syllabus), with an extension on trees at the end of the post.

consider the following (simulated dataset)

n=200
set.seed(1)
x1=runif(n)
x2=runif(n)
y=1+2*x1-x2+rnorm(n,0,.2)
LB=sample(LETTERS[1:10])
b=data.frame(y=y,x1=x1,
  x2=cut(x2,breaks=
  c(-1,.05,.1,.2,.35,.4,.55,.65,.8,.9,2),
  labels=LB))
str(b)
'data.frame':	200 obs. of  3 variables:
 $ y : num  1.345 1.863 1.946 2.481 0.765 ...
 $ x1: num  0.266 0.372 0.573 0.908 0.202 ...
 $ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ...
table(b$x2)[LETTERS[1:10]]
 
 A  B  C  D  E  F  G  H  I  J 
11 12 23 34 23 36 12 32  3 14

Just by looking at the data (see the previous post), we could easily get the feeling that 10 levels was too much.

Following my post, Przemyslaw sent a comment suggesting to use

library(factorMerger)

It is indeed a nice package (unless you have really really big datasets with a lot of categories in your factor variables – as I experienced recently), and you can get great graphs

MF = mergeFactors(response = b$y, 
             factor = b$x2, 
             family = "gaussian")
plot(MF)

Here is suggests to create three categories. Recall that with student t-tests (changing the reference), we got

Another interesting package, by Piro Polo, is

library(tree.bins)

To use it, we simply call the following function, and we transform automatically our dataset : the continuous variables remain unchanged, and (possibly) categories of categorical variables are merged

b.bins = tree.bins(data=b, y=y)
str(b.bins)
Classes ‘data.table’ and 'data.frame':	200 obs. of  3 variables:
 $ y : num  1.345 1.863 1.946 2.481 0.765 ...
 $ x1: num  0.266 0.372 0.573 0.908 0.202 ...
 $ x2: chr  "Group.4" "Group.4" "Group.4" "Group.4" ...
 - attr(*, ".internal.selfref")= 
table(b.bins$x2)

Group.1 Group.2 Group.3 Group.4 
     23      35      26     116

here in four groups. To get the correspondance, use

tree.bins(data=b, y=y, return = "lkup.list")
[[1]]
   x2 Categories
1   E    Group.1
2   G    Group.2
3   C    Group.2
4   B    Group.3
5   J    Group.3
6   I    Group.4
7   A    Group.4
8   H    Group.4
9   F    Group.4
10  D    Group.4

(we have a list with one element, one dataframe, since there is only one factor variable). Cool, isn’t it ? I miss Przemyslaw’s plot, but this is rather quick, and efficient..

 

2018 Automne, STT5100

Régression sur une variable qualitative et ANOVA

Leave a comment

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

Davis=read.table(
  "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]
Davis=data.frame(Y=Davis$weight * 2.204622,
                 X1=Davis$sex)

On voulait estimer le modèle y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_iOn avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') 
 Y=Davis$Y

car la matrice \mathbf{X}^T\mathbf{X} est inversible (une fois que l’on enlève la constante)

 solve(t(X)%*%X)
            [,1]       [,2]
[1,] 0.008928571 0.00000000
[2,] 0.000000000 0.01136364

et donc l’estimateur par moindres carrés est (classiquement)\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

 solve(t(X)%*%X) %*% (t(X)%*%Y)
         [,1]
[1,] 125.4272
[2,] 167.3258

ce qui correspond effectivement à la sortie de R,

 reg=lm(Y~0+X1,data=Davis)
 summary(reg)
 
Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X1F  125.427      1.960   64.00   <2e-16 ***
X1M  167.326      2.211   75.68   <2e-16 ***

Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes

x=Y[X[,1]==1]
y=Y[X[,2]==1]
nx=length(x)
ny=length(y)

On avait vu en cours que les \widehat{\mathbf{\beta}} avaient une interprétation très simple, puisque\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_iautrement dit \widehat{{\beta}}_M   est le poids moyen des hommes. Et en effet

 mean(y)
[1] 167.3258

C’est finalement très naturel, ou intuitif.

On peut maintenant s’interroger sur l’écart-type de l’estimateur de \widehat{{\beta}}_M . Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici

 sqrt(var(y)/ny)
[1] 2.794391
 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)
[1] 2.794391

car pour rappel\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de \mathbf{\beta} est proportionnel à \sigma^2 , la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de \sigma^2 

 s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2))
 sqrt(s2/ny)
[1] 2.210863

et en effet, on retombe sur la valeur donnée dans le tableau de régresion

 sqrt(s2/nx)
[1] 1.959721

(pareil pour l’autre coefficient).

On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).

 X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y)
          [,1]
[1,] 125.42724
[2,]  41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)
[1] 41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  125.427      1.960   64.00   <2e-16 ***
X1M           41.899      2.954   14.18   <2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -6.4475, df = 286, p-value = 4.826e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -30.62603 -16.30035
sample estimates:
mean of x mean of y 
 143.8626  167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable \mathbf{y} ), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Call:
lm(formula = Y ~ X1, data = Davis)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-57.458 -22.184  -5.512  17.809 118.912 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 143.4382     2.7820   51.56   <2e-16 ***
X1M           0.9645     4.1940    0.23    0.818    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 29.44 on 198 degrees of freedom
Multiple R-squared:  0.000267,	Adjusted R-squared:  -0.004782 
F-statistic: 0.05289 on 1 and 198 DF,  p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y
 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -0.22998, df = 198, p-value = 0.8183
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.235209  7.306165
sample estimates:
mean of x mean of y 
 143.4382  144.4027

et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents.

Regression

Visualizing effects of a categorical explanatory variable in a regression

Leave a comment

Recently, I’ve been working on two problems that might be related to semiotic issues in predictive modeling (i.e. instead of a standard regression table, how can we plot coefficient values in a regression model). To be more specific, I have a variable of interest Y that is observed for several individuals i, with explanatory variables \mathbf{x}_i, year t, in a specific region z_i\in\{A,B,C,D,E\}. Suppose that we have a simple (standard) linear model (forget about time here) y_i=\beta_0+\beta_1x_{1,i}+\cdots+\beta_kx_{k,i}+\sum_j \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i

Let us forget the temporal effect to focus on the spatial effect today. And consider some simulated dataset. There will be only one (continuous) explanatory variable. And I will generate correlated covariates, just to be more realistic.

n=1000
library(mnormt)
r=.5
Sigma=matrix(c(1,r,r,1), 2, 2)
set.seed(1)
X=rmnorm(n,c(0,0),Sigma)
X1=cut(X[,1],c(-100,quantile(X[,1],c(.1,.4,.7,.85)),
100),labels=LETTERS[1:5])
X2=X[,2]
Y=5+X[,1]-X[,2]+rnorm(n)/2
db=data.frame(Y,X1,X2)

Here we have y_i=\beta_0+\beta_1x_{1,i}+\sum_{j\in\{A,B,C,D,E\}} \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i The goal here is to get to graph to visualize the vector \hat\alpha=(\hat\alpha_A,\cdots,\hat\alpha_E). Let us run the linear regression

reg1=lm(Y~X1+X2,data=db)
idx=which(substr(names(reg1$coefficients), 1,2)=="X1")
v1=reg1$coefficients[idx]
names(v1)=LETTERS[2:5]
barplot(v1,col=rgb(0,0,1,.4))

Note that it is possible to add some sort of “confidence interval” to discuss significance (or to avoid to spend hours discussing differences in bar heights that are not significantly different)

library(Hmisc)
sv1=summary(reg1)$coefficients[idx,2]
(bp1=barplot(v1,ylim=range(c(0,v1+2*sv1))))
errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

My main concern here is the “reference” that is considered. Should A be the reference? Why not B

db$X1=relevel(db$X1,"B")
reg1=lm(Y~X1+X2,data=db)
idx=which(substr(names(reg1$coefficients),1,2)=="X1")
v1=reg1$coefficients[idx]
names(v1)=LETTERS[c(1,3:5)]
library(Hmisc)
sv1=summary(reg1)$coefficients[idx,2]
(bp1=barplot(v1)
errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

Why not the smallest one? Why not the largest one?… What if there is no simple way to choose. Furthermore, let us get back to the original point, which is that there might be some temporal aspects. More precisely, we can have \hat\alpha^{(t)}=(\hat\alpha_A^{(t)},\cdots,\hat\alpha_E^{(t)}). If we have also \hat\alpha^{(t+1)} and we get another plot, how do we interpret it. If for E the bar is taller, it means that relative to A, the difference has increased. I have the feeling that the interpretation is more complicated because we do not see, on that graph, changes in \hat\alpha^{(t)}_A.

Let us try something else. First, let us get back to the original setting

db$X1=relevel(db$X1,"A")

Consider here the regression without the intercept, so that all values remain

reg1=lm(Y~0+X1+X2,data=db)
idx=which(substr(names(reg1$coefficients),1,2)=="X1")
v1=reg2$coefficients[idx]
names(v1)=LETTERS[1:5]
barplot(v1)

It can be hard to read, especially if Y takes (very) large values, and you think that barplots should start at 0. But still, having those 5 values is nice. Why not rescale that graph?

A natural idea my be to consider the case where no spatial component is considered, and to look at the difference with that reference.

reg1=lm(Y~1+X2,data=db)
reg2=lm(Y~0+X1+X2,data=db)
idx=which(substr(names(reg2$coefficients),1,2)=="X1")
v1=reg2$coefficients[idx]
v2=v1-reg1$coefficients["(Intercept)"]
barplot(v2,col=rgb(0,0,1,.4))
sv2=summary(reg2)$coefficients[idx,2]
(bp2=barplot(v2,ylim=range(c(v2-2*sv2,v2+2*sv2))))
errbar(bp2[,1],v2,v2-2*sv2,v2+2*sv2,add=TRUE)

I like that graph, I should admit it. Now, I still have some remaining questions. For instance, can we insure that when only the intercept is considered, the value of \hat\beta_0 is somewhere between \hat\beta_A,\cdots,\hat\beta_E? Is it possible that \hat\beta_A-\hat\beta_0,\cdots,\hat\beta_E-\hat\beta_0 are all positive? In that case, I would find that hard to interpret.

Actually, if I really want values that can be seen as compared to some average, why not consider a (weighted) average of \hat\beta_A,\cdots,\hat\beta_E? (weights being here proportion in each class, in each region)

w=table(db$X1)
v3=v1-sum(w*v1)/sum(w)
(bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3))))
errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE)

I like that one. But what if, instead of normalizing at the end, we normalize the original dependent variable. By “normalize”, I mean “rescale”, to have a centered variable.

db$Y0=db$Y-mean(db$Y)
reg3=lm(Y0~0+X1+X2,data=db)
sv3=summary(reg3)$coefficients[idx,2]
(bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3))))
errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE)

This one is nice, because it is extremely simple to explain. But what if instead of a linear regression, we add a logistic one (with Y\in\{0,1\})? or a Poisson regression…

So maybe it cannot be the best solution here. Let us try something else… In insurance ratemaking, people like to use “zonier“. It is a two-stage regression. The idea is to run a regression without any spatial components, first. Then, consider the regression of residuals on spatial variables. Here, it would be something like

reg1=lm(Y~1+X2,data=db)
reg4=lm(residuals(reg1)~0+db$X1)

Since we focus on residuals, those are centered, and we have an easy interpretation of respective values

sv4=summary(reg4)$coefficients[idx,2]
v4=reg4$coefficients
(bp4=barplot(v4,names.arg=LETTERS[1:5])))
errbar(bp4[,1],v4,v4-2*sv4,v4+2*sv4,add=TRUE)

I guess that it can also be use in generalized linear models, with Pearson (or deviance) residuals.

Another possible idea can be the following. Again, the goal is not to have the true values, but to visualize on a graph how regions can be different. Here, all of them are significantly different. And in region A, Y is smaller, ceteris paribus (other things equal in the sense that we have taken into account x_1). And in region E it is larger. Here, the graph helps to “see” those differences.

Why not consider a completely different graph. What if we plot vector a instead of \alpha, where a_A can be interpreted as the value of the coefficient if we consider region A against “not region A“. What if we consider 5 regressions where dichotomous versions of Z are considered : Z_j=\mathbf{1}_{Z=j}.

v5=sv5=rep(NA,5)
names(v5)=LETTERS[1:5]
for(k in 1:5){
reg=lm(Y~I(X1==LETTERS[k])+X2,data=db)
v5[k]=reg$coefficients[2]
sv5[k]=summary(reg)$coefficients[2,2]}

We can plot that sequence of values, including some confidence intervals (that would be related to significance with respect to all other regions)

(bp5=barplot(v5,ylim=range(c(v5-2*sv5,v5+2*sv5))))
errbar(bp5[,1],v5,v5-2*sv5,v5+2*sv5,add=TRUE)

Looking at values does not give intuitive results, but I have the feeling that it is easy to explain what we plot (we compare each region to “the rest of the world”), and the ordering of a seems to be consistent with \alpha (but I could not prove it).

Here are some ideas I got. I should be able to provide other graphs, but I would love to discuss with anyone on that topics, to find a proper and nice way to visualize effects of a categorical explanatory variable in a regression model (that can be a logistic one). Comments are open…

Statistics

How Could Classification Trees Be So Fast on Categorical Variables?

1 Comment

I think that over the past months, I have been saying non-correct things about classification with categorical covariates. Because I never took time to look at it carefuly. Consider some simulated dataset, with a logistic regression,

> n=1e3
> set.seed(1)
> X1=runif(n)
> q=quantile(X1,(0:26)/26)
> q[1]=0
> X2=cut(X1,q,labels=LETTERS[1:26])
> p=exp(-.1+qnorm(2*(abs(.5-X1))))/(1+exp(-.1+qnorm(2*(abs(.5-X1)))))
> Y=rbinom(n,size=1,p)
> df=data.frame(X1=X1,X2=X2,p=p,Y=Y)

Here, we use some continuous covariate, except that is considered as not-observed. Instead, we have a categorical covariate with 26 categories. The (theoretical) relationship between the covariate and the probability is given below,

> vx1=seq(0,1,by=.001)
> vp=exp(-.1+qnorm(2*(abs(.5-vx1))))/(1+exp(-.1+qnorm(2*(abs(.5-vx1)))))
> plot(vx1,vp,type="l")

and the empirical probability, for each modality is

If we run a classification tree, we get

> library(rpart)
> tree=rpart(Y~X2,data=df)
> library(rpart.plot)
> prp(tree, type=2, extra=1)

To be more specific, the output is here

> tree
1) root 1000 249.90000 0.4900000  
  2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.3 0.302
    4) X2=J,K,L,M,N,O,P,Q,R 346  65.12 0.25144  *
    5) X2=F,G,H,I 153  37.22876 0.4183007       *
  3) X2=A,B,C,D,E,S,T,U,V,W,X,Y,Z 501 109.61 0.67
    6) X2=B,C,D,E,S,T,U,V,W,X 385  90.38 0.623  *
    7) X2=A,Y,Z 116  14.50862 0.8534483         *

 

Note that it takes less than a second to get that output. So clearly, we did not look for all combinations between modalities. For the first node, there are like  possible groups, i.e.

> 67108864

It is big… not huge, but too big to try all combinations, since that’s only the first node, and we have to do it again on the two leaves, etc. Antoine (aka @ly_antoine) told me – while we were having a coffee after lunch today – the trick to get a fast algorithm, on categories. And as usual, the idea is very clever…

First, we need a function to compute Gini index

> gini=function(y,classe){
+    T=table(y,classe)
+    nx=apply(T,2,sum)
+    n=sum(T)
+    pxy=T/matrix(rep(nx,each=2),nrow=2)
+    omega=matrix(rep(nx,each=2),nrow=2)/n
+    g=-sum(omega*pxy*(1-pxy))
+    return(g)}

For the first node, the idea is very simple:

  • Compute empirical averages 
> cond_prob=aggregate(df$Y,by=list(df$X2),mean)
  • Then sort those values, ,
  • Based on that ordering, consider 
> Group_Letters=cond_prob[order(cond_prob$x),2]

  • Then consider (only)  possible partitions,

against 

> v_gini=rep(NA,26)
> for(v in 1:26){
+   CLASSE=df$X2 %in% Group_Letters[1:v]
+   v_gini[v]=gini(y=df$Y,classe=CLASSE)
+ }

If we plot them, we get

> plot(1:26,v_gini,type="b)

As for continuous variables, we seek for the maximum value, and then, we have our two groups,

> sort(Group_Letters[1:which.max(v_gini)])
 [1] F G H I J K L M N O P Q R

That’s exactly what we got with the tree function in R,

1) root 1000 249.90000 0.4900000  
  2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30

Now, consider the leaf on the left (for instance)

> sub_df=df[df$X2 %in% sort(Group_Letters[1:which.max(v_gini)]),]

Then use the same algorithm as before: sort the conditional means,

> cond_prob=aggregate(sub_df$Y,by=
+ list(sub_df$X2),mean)
> s_Group_Letters=cond_prob[order(cond_prob$x),2]

Then compute Gini indices based on groups obtained from that ordering,

> v_gini=rep(NA,length(sub_Group_Letters))
> for(v in 1:length(sub_Group_Letters)){
+   CLASSE=sub_df$X2 %in% s_Group_Letters[1:v]
+   v_gini[v]=gini(y=sub_df$Y,classe=CLASSE)
+ }

If we plot it, we get our two groups,

> plot(1:length(s_Group_Letters),v_gini,type="b")

And the first group is here

> sort(sub_Group_Letters[1:which.max(v_gini)])
[1] J K L M N O P Q R

Again, that’s exactly what we got with the R function

1) root 1000 249.90000 0.4900000  
  2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30
    4) X2=J,K,L,M,N,O,P,Q,R 346  65.12 0.25144  *

Clever, isn’t?

2013, ACT2040, Actuarial science, Automne, Course, GLM, Ratemaking

Regression on variables, or on categories?

Leave a comment

I admit it, the title sounds weird. The problem I want to address this evening is related to the use of the stepwise procedure on a regression model, and to discuss the use of categorical variables (and possible misinterpreations). Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db2.txt",header=TRUE,sep=";")

First, let us change the reference in our categorical variable  (just to get an easier interpretation later on)

> db$X3=relevel(as.factor(db$X3),ref="E")

If we run a logistic regression on the three variables (two continuous, one categorical), we get

> reg=glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0758   0.1226   0.2805   0.4798   2.0345  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.39528    0.86649  -6.227 4.77e-10 ***
X1           0.51618    0.09163   5.633 1.77e-08 ***
X2           0.24665    0.05911   4.173 3.01e-05 ***
X3A         -0.09142    0.32970  -0.277   0.7816    
X3B         -0.10558    0.32526  -0.325   0.7455    
X3C          0.63829    0.37838   1.687   0.0916 .  
X3D         -0.02776    0.33070  -0.084   0.9331    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 806.29  on 999  degrees of freedom
Residual deviance: 582.29  on 993  degrees of freedom
AIC: 596.29

Number of Fisher Scoring iterations: 6

Now, if we use a stepwise procedure, to select variables in the model, we get

> step(reg)
Start:  AIC=596.29
Y ~ X1 + X2 + X3

       Df Deviance    AIC
- X3    4   587.81 593.81
<none>      582.29 596.29
- X2    1   600.56 612.56
- X1    1   617.25 629.25

Step:  AIC=593.81
Y ~ X1 + X2

       Df Deviance    AIC
<none>      587.81 593.81
- X2    1   606.90 610.90
- X1    1   622.44 626.44

So clearly, we should remove the categorical variable if our starting point was the regression on the three variables.

Now, what if we consider the same model, but slightly different: on the five categories,

> X3complete = model.matrix(~0+X3,data=db)
> db2 = data.frame(db,X3complete)
> head(db2)
  Y       X1       X2 X3 X3A X3B X3C X3D X3E
1 1 3.297569 16.25411  B   0   1   0   0   0
2 1 6.418031 18.45130  D   0   0   0   1   0
3 1 5.279068 16.61806  B   0   1   0   0   0
4 1 5.539834 19.72158  C   0   0   1   0   0
5 1 4.123464 18.38634  C   0   0   1   0   0
6 1 7.778443 19.58338  C   0   0   1   0   0

From a technical point of view, it is exactly the same as before, if we look at the regression,

> reg = glm(Y~X1+X2+X3A+X3B+X3C+X3D+X3E,family=binomial,data=db2)
> summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3A + X3B + X3C + X3D + X3E, family = binomial, 
    data = db2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0758   0.1226   0.2805   0.4798   2.0345  

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.39528    0.86649  -6.227 4.77e-10 ***
X1           0.51618    0.09163   5.633 1.77e-08 ***
X2           0.24665    0.05911   4.173 3.01e-05 ***
X3A         -0.09142    0.32970  -0.277   0.7816    
X3B         -0.10558    0.32526  -0.325   0.7455    
X3C          0.63829    0.37838   1.687   0.0916 .  
X3D         -0.02776    0.33070  -0.084   0.9331    
X3E               NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 806.29  on 999  degrees of freedom
Residual deviance: 582.29  on 993  degrees of freedom
AIC: 596.29

Number of Fisher Scoring iterations: 6

Both regressions are equivalent. Now, what about a stepwise selection on this new model?

> step(reg)
Start:  AIC=596.29
Y ~ X1 + X2 + X3A + X3B + X3C + X3D + X3E

Step:  AIC=596.29
Y ~ X1 + X2 + X3A + X3B + X3C + X3D

       Df Deviance    AIC
- X3D   1   582.30 594.30
- X3A   1   582.37 594.37
- X3B   1   582.40 594.40
<none>      582.29 596.29
- X3C   1   585.21 597.21
- X2    1   600.56 612.56
- X1    1   617.25 629.25

Step:  AIC=594.3
Y ~ X1 + X2 + X3A + X3B + X3C

       Df Deviance    AIC
- X3A   1   582.38 592.38
- X3B   1   582.41 592.41
<none>      582.30 594.30
- X3C   1   586.30 596.30
- X2    1   600.58 610.58
- X1    1   617.27 627.27

Step:  AIC=592.38
Y ~ X1 + X2 + X3B + X3C

       Df Deviance    AIC
- X3B   1   582.44 590.44
<none>      582.38 592.38
- X3C   1   587.20 595.20
- X2    1   600.59 608.59
- X1    1   617.64 625.64

Step:  AIC=590.44
Y ~ X1 + X2 + X3C

       Df Deviance    AIC
<none>      582.44 590.44
- X3C   1   587.81 593.81
- X2    1   600.73 606.73
- X1    1   617.66 623.66

What do we get now? This time, the stepwise procedure recommends that we keep one category (namely C). So my point is simple: when running a stepwise procedure with factors, either we keep the factor as it is, or we drop it. If it is necessary to change the design, by pooling together some categories, and we forgot to do it, then it will be suggested to remove that variable, because having 4 categories meaning the same thing will cost us too much if we use the Akaike criteria. Because this is exactly what happens here

> library(car)
> reg = glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)
> linearHypothesis(reg,c("X3A=X3B","X3A=X3D","X3A=0"))
Linear hypothesis test

Hypothesis:
X3A - X3B = 0
X3A - X3D = 0
X3A = 0

Model 1: restricted model
Model 2: Y ~ X1 + X2 + X3

  Res.Df Df  Chisq Pr(>Chisq)
1    996                     
2    993  3 0.1446      0.986

So here, we should pool together categories A, B, D and E (which was here the reference). As mentioned in a previous post, it is necessary to pool together categories that should be pulled together as soon as possible. If not, the stepwise procedure might yield to some misinterpretations.

2013, ACT2040, Actuarial science, Automne, Course, GLM, Tests

Logistic regression and categorical covariates

8 Comments

A short post to get back – for my nonlife insurance course – on the interpretation of the output of a regression when there is a categorical covariate. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)
> tail(db)
     Y       X1       X2 X3
995  1 4.801836 20.82947  A
996  1 9.867854 24.39920  C
997  1 5.390730 21.25119  D
998  1 6.556160 20.79811  D
999  1 4.710276 21.15373  A
1000 1 6.631786 19.38083  A

Let us run a logistic regression on that dataset

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.45885    1.04646  -4.261 2.04e-05 ***
X1           0.51664    0.11178   4.622 3.80e-06 ***
X2           0.21008    0.07247   2.899 0.003745 ** 
X3B          1.74496    0.49952   3.493 0.000477 ***
X3C         -0.03470    0.35691  -0.097 0.922543    
X3D          0.08004    0.34916   0.229 0.818672    
X3E          2.21966    0.56475   3.930 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 397.69  on 993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

Here, the reference is modality . Which means that for someone with characteristics , we predict the following probability

where  denotes the cumulative distribution function of the logistic distribution

For someone with characteristics , we predict the following probability

For someone with characteristics , we predict the following probability

(etc.) Here, if we accept  (against ), it means that modality  cannot be considerd as different from .

A natural idea can be to change the reference modality, and to look at the -values. If we consider the following loop, we get

> M = matrix(NA,5,5)
> rownames(M)=colnames(M)=LETTERS[1:5]
> for(k in 1:5){
+ db$X3 = relevel(X3,LETTERS[k])
+ reg = glm(Y~X1+X2+X3,family=binomial,data=db)
+ M[levels(db$X3)[-1],k] = summary(reg)$coefficients[4:7,4]
+ } 
> M
             A            B            C            D            E
A           NA 0.0004771853 9.225428e-01 0.8186723647 8.482647e-05
B 4.771853e-04           NA 4.841204e-04 0.0009474491 4.743636e-01
C 9.225428e-01 0.0004841204           NA 0.7506242347 9.194193e-05
D 8.186724e-01 0.0009474491 7.506242e-01           NA 1.730589e-04
E 8.482647e-05 0.4743636442 9.194193e-05 0.0001730589           NA

and if we simply want to know if the -value exceeds – or not – 5%, we get the following,

> M.TF = M>.05
> M.TF
      A     B     C     D     E
A    NA FALSE  TRUE  TRUE FALSE
B FALSE    NA FALSE FALSE  TRUE
C  TRUE FALSE    NA  TRUE FALSE
D  TRUE FALSE  TRUE    NA FALSE
E FALSE  TRUE FALSE FALSE    NA

The first column is obtained when  is the reference, and then, we see which parameter should be considered as null. The interpretation is the following:

  •  and  are not different from 
  •  is not different from 
  •  and  are not different from 
  •  and  are not different from 
  •  is not different from 

Note that we only have, here, some kind of intuition. So, let us run a more formal test. Let us consider the following regression (we remove the intercept to get a model easier to understand)

> library(car)
> db$X3=relevel(X3,"A")
> reg=glm(Y~0+X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
    Estimate Std. Error z value Pr(>|z|)    
X1   0.51664    0.11178   4.622 3.80e-06 ***
X2   0.21008    0.07247   2.899  0.00374 ** 
X3A -4.45885    1.04646  -4.261 2.04e-05 ***
X3E -2.23919    1.06666  -2.099  0.03580 *  
X3D -4.37881    1.04887  -4.175 2.98e-05 ***
X3C -4.49355    1.06266  -4.229 2.35e-05 ***
X3B -2.71389    1.07274  -2.530  0.01141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1386.29  on 1000  degrees of freedom
Residual deviance:  397.69  on  993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

It is possible to use Fisher test to test if some coefficients are equal, or not (more generally if some linear constraints are satisfied)

> linearHypothesis(reg,c("X3A=X3C","X3A=X3D","X3B=X3E"))
Linear hypothesis test

Hypothesis:
X3A - X3C = 0
X3A - X3D = 0
- X3E  + X3B = 0

Model 1: restricted model
Model 2: Y ~ 0 + X1 + X2 + X3

  Res.Df Df  Chisq Pr(>Chisq)
1    996                     
2    993  3 0.6191      0.892

Here, we clearly accept the assumption that the first three factors are equal, as well as the last two. What is the next step? Well, if we believe that there are mainly two categories,  and , let us create that factor,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.39439    1.02791  -4.275 1.91e-05 ***
X1           0.51378    0.11138   4.613 3.97e-06 ***
X2           0.20807    0.07234   2.876  0.00402 ** 
X3bisBE      1.94905    0.36852   5.289 1.23e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 398.31  on 996  degrees of freedom
AIC: 406.31

Number of Fisher Scoring iterations: 7

Here, all the categories are significant. So we do have a proper model.