{"word_count": 575, "character_count": 6515}
This morning, Stéphane asked me tricky question about extracting coefficients from a regression with categorical explanatory variates. More precisely, he asked me if it was possible to store the coefficients in a nice table, with information on the variable and the modality (those two information being in two different columns). Here is some code I did to produce the table he was looking for, but I guess that some (much) smarter techniques can be used (comments - see below - are open). Consider the following dataset&#62; base   x sex   hair1  1   H  Black2  4   F  Brown3  6   F  Black4  6   H  Black5 10   H  Brown6  5   H Blondewith two factors,&#62; levels(base$hair)[1] "Black"  "Blonde" "Brown" &#62; levels(base$sex)[1] "F" "H"Let us run a (standard linear) regression,&#62; reg=lm(x~hair+sex,data=base)which is here&#62; summary(reg)Call:lm(formula = x ~ hair + sex, data = base)Residuals:         1          2          3          4          5          6 -3.714e+00 -2.429e+00  2.429e+00  1.286e+00  2.429e+00 -2.220e-16 Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)(Intercept)   3.5714     3.4405   1.038    0.408hairBlonde    0.2857     4.8655   0.059    0.959hairBrown     2.8571     3.7688   0.758    0.528sexH          1.1429     3.7688   0.303    0.790Residual standard error: 4.071 on 2 degrees of freedomMultiple R-squared: 0.2352,Adjusted R-squared: -0.9121 F-statistic: 0.205 on 3 and 2 DF,  p-value: 0.886If we want to extract the names of the factors (assuming here that there are no numbers in the name of the factor), and the values of the associated modality, one can use&#62; VARIABLE=c("",gsub("[-^0-9]", "", names(unlist(reg$xlevels))))&#62; MODALITY=c("",as.character(unlist(reg$xlevels)))&#62; names=data.frame(VARIABLE,MODALITY,NOMVAR=c(+ "(Intercept)",paste(VARIABLE,MODALITY,sep="")[-1]))&#62; regression=data.frame(NOMVAR=names(coefficients(reg)),+ COEF=as.numeric(coefficients(reg)))&#62; merge(names,regression,all.x=TRUE)       NOMVAR VARIABLE MODALITE      COEF1 (Intercept)                   3.57142862   hairBlack     hair    Black        NA3  hairBlonde     hair   Blonde 0.28571434   hairBrown     hair    Brown 2.85714295        sexF      sex        F        NA6        sexH      sex        H 1.1428571or, if we want modalities exluding references,&#62; merge(names,regression)       NOMVAR VARIABLE MODALITE      COEF1 (Intercept)                   3.57142862  hairBlonde     hair   Blonde 0.28571433   hairBrown     hair    Brown 2.85714294        sexH      sex        H 1.1428571In order to reproduce the table Stéphane sent me, let us use the following code to produce an html table,&#62; library(xtable)&#62; htlmtable &#60;- xtable(merge(names,regression))&#62; print(htlmtable,type="html")NOMVARVARIABLEMODALITYCOEF1(Intercept)3.572hairBlondehairBlonde0.293hairBrownhairBrown2.864sexHsexH1.14So yes, it is possible to build a table with the variable, modalities, and coefficients. This function can be interesting on prospective mortality, when we do have a large number of modalities per factor (years, ages and year of birth). Consider the following datasets&#62; DEATH=read.table(+ "http://freakonometrics.free.fr/DeathsSwitzerland.txt",+ header=TRUE,skip=2)&#62; EXPOSURE=read.table(+ "http://freakonometrics.free.fr/ExposuresSwitzerland.txt",+ header=TRUE,skip=2)&#62; DEATH$Age=as.numeric(as.character(DEATH$Age))&#62; DEATH=DEATH[-which(is.na(DEATH$Age)),]&#62; EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age))&#62; EXPOSURE=EXPOSURE[-which(is.na(EXPOSURE$Age)),]&#62; base=data.frame(y=as.factor(DEATH$Year),a=as.factor(DEATH$Age),+ c=as.factor(DEATH$Year-DEATH$Age),D=DEATH$Total,E= EXPOSURE$Total)&#62; base=base[base$E&#62;0,]and the following nonlinear model, based on Lee-Carter model (including a cohort effect),can be estimated using&#62; library(gnm)&#62; reg=gnm(D~a+Mult(a,y)+Mult(a,c),offset=log(E),family=poisson,data=base)In order to extract the 671 coefficients from the regresssion,&#62; length(coefficients(reg))[1] 671(as properly as possible) we have to be careful: names of coefficients are not that simple to handle. For instance, we can see things like&#62; coefficients(reg)[200]Mult(., year).age98          0.04203519In order to extract them, define&#62; na=length((reg$xlevels)$age)&#62; ny=length((reg$xlevels)$year)&#62; nc=length((reg$xlevels)$cohort)&#62; VARIABLElong=c("",rep("age",na),rep("Mult(., year).age",na),+ rep("Mult(a, .).y",ny),+ rep("Mult(., cohort).age",na),rep("Mult(age, .).cohort",nc))&#62; VARIABLEshort=c("",rep("age",na),rep("age",na),rep("year",ny),+ rep("age",na),rep("cohort",nc))&#62; MODALITY=c("",(reg$xlevels)$age,(reg$xlevels)$age,+ (reg$xlevels)$year,(reg$xlevels)$age,(reg$xlevels)$cohort)&#62; names=data.frame(VARIABLElong,VARIABLEshort,+ MODALITY,NOMVAR=c("(Intercept)",paste(VARIABLElong,MODALITY,sep="")[-1]))&#62; regression=data.frame(NOMVAR=names(coefficients(reg)),+ COEF=as.numeric(coefficients(reg)))Here we go, now we have the coefficients from the regression in a nice table,&#62; outputreg=merge(names,regression)&#62; outputreg[1:10,]        NOMVAR VARIABLElong VARIABLEshort MODALITY        COEF1  (Intercept)                                     -8.222254582         age1          age           age        1 -0.874954513        age10          age           age       10 -1.671457044       age100          age           age      100  4.910416505        age11          age           age       11 -1.001869906        age12          age           age       12 -1.059534977        age13          age           age       13 -0.909528598        age14          age           age       14  0.028806689        age15          age           age       15  0.4283073810       age16          age           age       16  1.35961403It is now possible to plot all the coefficients, as functions of the age, the year of observation, or the year of birth. For instance, for the standard average age effect (namely  as a function of ), we can use&#62; typevariable=as.character(unique(outputreg$VARIABLElong))&#62; basegraph=outputreg[outputreg$VARIABLElong==typevariable[2],]&#62; x=as.numeric(as.character(basegraph$MODALITY))&#62; y=basegraph$COEF&#62; plot(x,y,type="p",col="blue",xlab="Age")while the cohort effect ( as a function of ) is obtained using&#62; basegraph=outputreg[outputreg$VARIABLElong==typevariable[5],]&#62; x=as.numeric(as.character(basegraph$MODALITY))&#62; y=basegraph$COEF&#62; plot(x,y,type="p",col="blue",xlab="Cohort (year of birth)",ylim=c(0,10))