{"word_count": 605, "character_count": 4078}
This morning, Benoit sent me an email, about an exercise he found in an econometric textbook, about linear regression. Consider the following dataset,Here, variable X denotes the income, and Y the expenses. The goal was to fit a linear regression (actually, in the email, it was mentioned that we should try to fit an heteroscedastic model, but let us skip this part). So Benoit's question was more or less: how do you fit a linear regression from a contingency table?Usually, when I got an email on Saturday morning, I try to postpone. But the kids had their circus class, so I had some time to answer. And this did not look like a complex puzzle... Let us import this dataset in R, so that we can start playing&#62; df=read.table("http://freakonometrics.free.fr/baseexo.csv",sep=";",header=TRUE)&#62; M=as.matrix(df[,2:ncol(df)])&#62; M[is.na(M)]&#60;-0&#62; M      X14 X19 X21 X23 X25 X27 X29 X31 X33 X35 [1,]  74  13   7   1   0   0   0   0   0   0 [2,]   6   4   2   7   4   0   0   0   0   0 [3,]   2   3   2   2   4   0   0   0   0   0 [4,]   1   1   2   3   3   2   0   0   0   0 [5,]   2   0   1   3   2   0   6   0   0   0 [6,]   2   0   2   1   0   0   1   2   1   0 [7,]   0   0   0   2   0   0   1   1   3   0 [8,]   0   1   0   1   0   0   0   0   2   0 [9,]   0   0   0   0   1   1   0   1   0   1The first idea I had was to use those counts as weights. Weighted least squares should be perfect. The dataset is built from this matrix,&#62; W=as.vector(M)&#62; x=df[,1]&#62; X=rep(x,ncol(M))&#62; y=as.numeric(substr(names(df)[-1],2,3))&#62; Y=rep(y,each=nrow(M))&#62; base1=data.frame(X1=X,Y1=Y,W1=W)Here we have&#62; head(base1,10)   X1 Y1 W11  16 14 742  23 14  63  25 14  24  27 14  15  29 14  26  31 14  27  33 14  08  35 14  09  37 14  010 16 19 13The regression is the following,&#62; summary(reg1)Call:lm(formula = Y1 ~ X1, data = base1, weights = W1)Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)    (Intercept)  4.35569    2.03022   2.145    0.038 *  X1           0.68263    0.09016   7.572 3.04e-09 ***---Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 7.892 on 40 degrees of freedomMultiple R-squared:  0.589,Adjusted R-squared:  0.5787 F-statistic: 57.33 on 1 and 40 DF,  p-value: 3.038e-09It looks like the output is the same as what Benoit found, so we should be happy. Now, I had a second thought. Why not create the implied dataset. Using replicates, we should be able to create the dataset that was used to get this contingency table,&#62; vX=vY=rep(NA,sum(W))&#62; sumW=c(0,cumsum(W))&#62; for(i in 1:length(W)){+ if(W[i]&#62;0){+ vX[(1+sumW[i]):sumW[i+1]]=X[i]+ vY[(1+sumW[i]):sumW[i+1]]=Y[i]+ }}&#62; base2=data.frame(X2=vX,Y2=vY)Here, the dataset is much larger, and there is no weight,&#62; tail(base2,10)    X2 Y2172 31 31173 33 31174 37 31175 31 33176 33 33177 33 33178 33 33179 35 33180 35 33181 37 35If we run a linear regression on this dataset, we obtain&#62; summary(reg2)Call:lm(formula = Y2 ~ X2, data = base2)Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)    (Intercept)  4.35569    0.95972   4.538 1.04e-05 ***X2           0.68263    0.04262  16.017  &#60; 2e-16 ***---Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 3.731 on 179 degrees of freedomMultiple R-squared:  0.589,Adjusted R-squared:  0.5867 F-statistic: 256.5 on 1 and 179 DF,  p-value: &#60; 2.2e-16If we compare the two regressions, we have&#62; rbind(coefficients(summary(reg1)),+ coefficients(summary(reg2)))             Estimate Std. Error   t value     Pr(&#62;|t|)(Intercept) 4.3556857 2.03021637  2.145429 3.804237e-02X1          0.6826296 0.09015771  7.571506 3.038443e-09(Intercept) 4.3556857 0.95972279  4.538483 1.036711e-05X2          0.6826296 0.04261930 16.016913 2.115373e-36The estimators are exactly the same (which does not surprise me), but standard deviation (ans significance levels) are quite different. And to be honest, I find that surprising. Which approach here is the most legitimate (since they are finally not equivalent)?