Regression without intercept in R and Stata

user_4266314 0
Viewed 850

Recently, I stumbled upon the fact that Stata and R handle regressions without intercept differently. I'm not a statistician, so please be kind if my vocabulary is not ideal.

I tried to make the example somewhat reproducible. This is my example in R:

> set.seed(20210211)
> df <- data.frame(y = runif(50), x = runif(50))
> df$d <- df$x > 0.5
> 
> (tmp <- tempfile("data", fileext = ".csv"))
[1] "C:\\Users\\s1504gl\\AppData\\Local\\Temp\\1\\RtmpYtS6uk\\data1b2c1c4a96.csv"
> write.csv(df, tmp, row.names = FALSE)
> 
> summary(lm(y ~ x + d, data = df))

Call:
lm(formula = y ~ x + d, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.48651 -0.27449  0.03828  0.22119  0.53347 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.4375     0.1038   4.214 0.000113 ***
x            -0.1026     0.3168  -0.324 0.747521    
dTRUE         0.1513     0.1787   0.847 0.401353    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2997 on 47 degrees of freedom
Multiple R-squared:  0.03103,   Adjusted R-squared:  -0.0102 
F-statistic: 0.7526 on 2 and 47 DF,  p-value: 0.4767

> summary(lm(y ~ x + d + 0, data = df))

Call:
lm(formula = y ~ x + d + 0, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.48651 -0.27449  0.03828  0.22119  0.53347 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
x       -0.1026     0.3168  -0.324 0.747521    
dFALSE   0.4375     0.1038   4.214 0.000113 ***
dTRUE    0.5888     0.2482   2.372 0.021813 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2997 on 47 degrees of freedom
Multiple R-squared:  0.7196,    Adjusted R-squared:  0.7017 
F-statistic: 40.21 on 3 and 47 DF,  p-value: 4.996e-13

And here is what I have in Stata (please note that I have copied the filename from R to Stata):


. import delimited "C:\Users\s1504gl\AppData\Local\Temp\1\RtmpYtS6uk\data1b2c1c4a96.csv"
(3 vars, 50 obs)

. encode d, generate(d_enc)

. 
. regress y x i.d_enc

      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 47)        =      0.75
       Model |  .135181652         2  .067590826   Prob > F        =    0.4767
    Residual |  4.22088995        47  .089806169   R-squared       =    0.0310
-------------+----------------------------------   Adj R-squared   =   -0.0102
       Total |   4.3560716        49   .08889942   Root MSE        =    .29968

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |  -.1025954   .3168411    -0.32   0.748    -.7399975    .5348067
             |
       d_enc |
       TRUE  |   .1512977   .1786527     0.85   0.401    -.2081052    .5107007
       _cons |   .4375371    .103837     4.21   0.000     .2286441    .6464301
------------------------------------------------------------------------------

. regress y x i.d_enc, noconstant

      Source |       SS           df       MS      Number of obs   =        50
-------------+----------------------------------   F(2, 48)        =     38.13
       Model |  9.23913703         2  4.61956852   Prob > F        =    0.0000
    Residual |  5.81541777        48  .121154537   R-squared       =    0.6137
-------------+----------------------------------   Adj R-squared   =    0.5976
       Total |  15.0545548        50  .301091096   Root MSE        =    .34807

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |    .976214   .2167973     4.50   0.000     .5403139    1.412114
             |
       d_enc |
       TRUE  |  -.2322011   .1785587    -1.30   0.200    -.5912174    .1268151
------------------------------------------------------------------------------

As you can see, the results of the regression with intercept are identical. But if I omit the intercept (+ 0 in R, , noconstant in Stata), the results differ. In R, the intercept is now captured in dFALSE, which is reasonable from what I understand. I don't understand what Stata is doing here. Also the degrees of freedom differ.

My questions:

  1. Can anyone explain to me how Stata is handling this?
  2. How can I replicate Stata's behavior in R?
stata lm linear-regression r
1 Answers

I believe bas pointed in the right direction, but I am still unsure why both results differ.

I am not attempting to answer the question, but provdide deeper understanding of what stata is doing (by digging into the source of R's lm() function. In the following lines I replicate what lm() does, but jumping over sanity checks and options such as weights, contrasts, etc...

(I cannot yet fully understand why in the second regression (with NO CONSTANT) the dFALSE coefficient captures the effect of the intercept in the default regression (with constant)

set.seed(20210211)
df <- data.frame(y = runif(50), x = runif(50))
df$d <- df$x > 0.5

lm() With Constant

form_default <- as.formula(y ~ x + d)
mod_frame_def <- model.frame(form_default, df)
mod_matrix_def <- model.matrix(object = attr(mod_frame_def, "terms"), mod_frame_def)
head(mod_matrix_def)
#>   (Intercept)         x dTRUE
#> 1           1 0.7861162     1
#> 2           1 0.2059603     0
#> 3           1 0.9793946     1
#> 4           1 0.8569093     1
#> 5           1 0.8124811     1
#> 6           1 0.7769280     1

stats:::lm.fit(
  y = model.response(mod_frame_def),
  x = mod_matrix_def
)$coefficients
#> (Intercept)           x       dTRUE 
#>   0.4375371  -0.1025954   0.1512977

lm() No Constant

form_nocon  <- as.formula(y ~ x + d + 0)
mod_frame_nocon <- model.frame(form_nocon, df)
mod_matrix_nocon <- model.matrix(object = attr(mod_frame_nocon, "terms"), mod_frame_nocon)
head(mod_matrix_nocon)
#>           x dFALSE dTRUE
#> 1 0.7861162      0     1
#> 2 0.2059603      1     0
#> 3 0.9793946      0     1
#> 4 0.8569093      0     1
#> 5 0.8124811      0     1
#> 6 0.7769280      0     1

stats:::lm.fit(
  y = model.response(mod_frame_nocon),
  x = mod_matrix_nocon
)$coefficients
#>          x     dFALSE      dTRUE 
#> -0.1025954  0.4375371  0.5888348

lm() with as.numeric()

[as indicated in the comments by bas]

form_asnum  <- as.formula(y ~ x + as.numeric(d) + 0)
mod_frame_asnum <- model.frame(form_asnum, df)
mod_matrix_asnum <- model.matrix(object = attr(mod_frame_asnum, "terms"), mod_frame_asnum)
head(mod_matrix_asnum)
#>           x as.numeric(d)
#> 1 0.7861162             1
#> 2 0.2059603             0
#> 3 0.9793946             1
#> 4 0.8569093             1
#> 5 0.8124811             1
#> 6 0.7769280             1

stats:::lm.fit(
  y = model.response(mod_frame_asnum),
  x = mod_matrix_asnum
)$coefficients
#>             x as.numeric(d) 
#>     0.9762140    -0.2322012

Created on 2021-03-18 by the reprex package (v1.0.0)

user_82371860
Related