R lm versus Python sklearn linear_model
<p>When I study Python SKlearn, the first example that I come across is the linear regression.</p>

<p><a href="http://scikit-learn.org/stable/modules/linear_model.html" rel="nofollow">http://scikit-learn.org/stable/modules/linear_model.html</a></p>

<p>Code of its very first example:</p>

<pre><code>from sklearn import linear_model
reg = linear_model.LinearRegression()
reg.fit([[0, 0], [1, 1], [2,2]], [0, 1,2])
reg.fit
reg.coef_
array([ 0.5,  0.5])
</code></pre>

<p>Here I assume <code>[[0, 0], [1, 1], [2,2]]</code> represents a data.frame containing <code>x1 = c(0,1,2)</code> and <code>x2 = c(0,1,2)</code> and <code>y = c(0,1,2)</code> as well.</p>

<p>Immediately, I begin to think that <code>array([ 0.5,  0.5])</code> are the coeffs for <code>x1</code> and <code>x2</code>. </p>

<p>But, are there standard errors for those estimates? What about t tests p values, R2 and other figures?</p>

<p>Then I try to do the same thing in R.</p>

<pre><code>X = data.frame(x1 = c(0,1,2),x2 = c(0,1,2),y = c(0,1,2))
lm(data=X, y~x1+x2)
Call:
lm(formula = y ~ x1 + x2, data = X)

#Coefficients:
#(Intercept)           x1           x2  
#  1.282e-16    1.000e+00           NA  
</code></pre>

<p>Obviously <code>x1</code> and <code>x2</code> are completely linearly dependent so the OLS will fail. Why the SKlearn still works and gives this results? Am I getting sklearn in a wrong way? Thanks.</p>

<p>Both solutions are correct (assuming that NA behaves like a zero). Which solution is favored depends on the numerical solver used by the OLS estimator.</p>

<p><code>sklearn.linear_model.LinearRegression</code> is based on <code>scipy.linalg.lstsq</code> which in turn calls the LAPACK <code>gelsd</code> routine which is described here:</p>

<p><a href="http://www.netlib.org/lapack/lug/node27.html" rel="nofollow">http://www.netlib.org/lapack/lug/node27.html</a></p>

<p>In particular it says that when the problem is rank deficient it seeks the minimum norm least squares solution.</p>

<p>If you want to favor the other solution, you can use a coordinate descent solver with a tiny bit of L1 penalty as implemented in th Lasso class:</p>

<pre><code>&gt;&gt;&gt; from sklearn.linear_model import Lasso
&gt;&gt;&gt; reg = Lasso(alpha=1e-8)
&gt;&gt;&gt; reg.fit([[0, 0], [1, 1], [2, 2]], [0, 1, 2])

Lasso(alpha=1e-08, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
&gt;&gt;&gt; reg.coef_
array([  9.99999985e-01,   3.97204719e-17])
</code></pre>

python scikit-learn regression linear