In [1]:
k = 2 # variables independientes
N = 10 # observaciones
Out[1]:
In [2]:
X = rand(N, k) # simulación: valores observados para las variables independientes
Out[2]:
In [3]:
α = rand() # simulación: valor real de la intersección
Out[3]:
In [4]:
β = rand(k) # simulación: valores reales de los coeficientes
Out[4]:
In [5]:
E = 0.1 .* rand(N) # simulación: agrega errores al azar al modelo
Out[5]:
In [6]:
Y = α + X * β + E # valores para la variable dependiente con el modelo lineal que acabamos de armar
Out[6]:
La biblioteca MultivariateStats de Julia
implementa la regresión lineal por cuadrados mínimos en su función llsq. Este método toma la matriz de valores de la variable dependiente y el vector de la independiente. Retorna un vector de la forma $[\beta_{1}, …, \beta_{k}, \alpha]$.
In [7]:
using MultivariateStats
predicción = llsq(X, Y)
Out[7]:
In [8]:
αₚ = predicción[end] # valor de la intersección del ajuste
Out[8]:
In [9]:
βₚ = predicción[1:end-1] # valores de los coeficientes del ajuste
Out[9]:
In [10]:
Yₚ = αₚ + X*βₚ # valores predichos para la variable dependiente
Out[10]:
Diferencias entre los valores reales y los ajustados
In [11]:
α - αₚ
Out[11]:
In [12]:
β - βₚ
Out[12]:
In [13]:
Eₚ = Y - Yₚ # son los valores residuales de la regresión
Out[13]:
In [14]:
using Plots
pyplot(size=(500,250))
plot(
scatter(Y, Yₚ, legend=false,
xlab="Valores", ylab="Predicciones"),
histogram(Eₚ, label="Residuales", bins=5)
)
Out[14]:
Supuestos, muchos de los cuales aplican a los residuos:
In [15]:
using GLM
In [16]:
Xₐ = hcat(ones(N), X) # Y = 1*α + β₁X₁ + β₂X₂ + ...
Out[16]:
In [17]:
modelo = lm(Xₐ,Y)
Out[17]:
In [18]:
residuos = residuals(modelo)
Out[18]:
En el caso de no cumplirse con la homocedacia, es posible darle a cada observación un peso inversamente proporcional a la varianza para ese valor o grupo de valores de la/s variable/s independiente/s en una regresion por minimos cuadrados ponderados. Lo óptimo es conocer de antemano cual es la varianza de cada medición. Por ejemplo si medimos cada dato con una herramienta diferente y conocemos su error. Sino existen maneras de determinar los posibles valores de la varianza (o de un estimador robusto correlacionado con ésta) con respecto a las variables independientes.
In [19]:
using DataFrames
In [20]:
N = 100
X = rand(N)
β = rand()
α = 5
E = rand(Normal(1,0.5),N) .* X # heterocedasticidad
Y = α + X * β + E
datos = DataFrame(hcat(X, Y))
names!(datos, [:X, :Y])
Out[20]:
In [25]:
modelo = glm(@formula(Y ~ X), datos, Normal(), IdentityLink())
Out[25]:
In [26]:
coeficientes = coef(modelo)
Out[26]:
In [29]:
using StatPlots
pyplot(size=(300,300))
scatter(datos, :X, :Y, legend=false)
Plots.abline!(coeficientes[2], coeficientes[1])
Out[29]:
In [30]:
residuos = Y - predict(modelo)
Out[30]:
In [31]:
scatter(X, residuos, legend=false)
Out[31]:
In [32]:
residuos_modelo = DataFrame(
x = collect(0:0.05:0.8).+0.1, # inicio + paso/2 = centro de la ventana
y = Float64[ mad( residuos[i .<= X .<= i+0.2] ) for i in collect(0:0.05:0.8) ] # paso = 0.2
)
Out[32]:
In [35]:
modelo = glm(@formula(y ~ x), residuos_modelo, Normal(), IdentityLink())
Out[35]:
In [36]:
plot(
scatter(residuos_modelo, :x,:y),
scatter(residuos_modelo[:x],residuos_modelo[:y]-predict(modelo)),
legend=false, size=(400,200))
Out[36]:
In [37]:
C = coef(modelo)
Out[37]:
In [38]:
pesos = Float64[ 1/(C[1] + C[2]*x) for x in datos[:X] ]
Out[38]:
In [39]:
lineal_ponderada = glm(@formula(Y ~ X), datos, Normal(), IdentityLink(), wts=pesos)
Out[39]:
In [40]:
reg_plot = scatter(datos, :X, :Y)
Plots.abline!(reg_plot, coef(lineal_ponderada)[2], coef(lineal_ponderada)[1])
residuos = ( datos[:Y] - predict(lineal_ponderada) ) .* pesos
normal = Normal(0,std(residuos))
res_plot = histogram(residuos, bins=10, normed=true)
plot!(res_plot, res -> pdf(normal, res), minimum(residuos), maximum(residuos))
res_scatter = scatter(datos[:X], residuos)
plot(reg_plot, res_scatter, res_plot, legend=false, size=(800,200), layout=grid(1,3))
Out[40]:
En general es preferible usar pesos cuando éstos se conocen de antemano. Por ejemplo si se conoce el error experimental de los instrumentos, B factors o resolución en estructuras cristalográficas, etc. Tratar de encontrar pesos puede ser complicado. En este caso, los pesos hacen que los residuales ponderados cumplan homocedacia, pero quizás a costa de la normalidad.
Ridge regression en MultivariateStats y en ScikitLearn.
Esta regresión es mucho más robusta a la presencia de variables explicativas correlacionadas (colineales).
Ridge Regression usando ScikitLearn.jl, el valor de α es seleccionado usando validación cruzada:
In [41]:
k = 3 # variables independientes
@assert k > 2
N = 10 # observaciones
X = rand(N, k) # simulación: valores observados para las variables independientes
X[:,1] = X[:,1].*X[:,2] + X[:,1] # colineales
α = rand() # simulación: valor real de la intersección
β = rand(k) # simulación: valores reales de los coeficientes
E = 0.1 .* rand(N) # simulación: agrega errores al azar al modelo
Y = α + X * β + E # valores para la variable dependiente con el modelo lineal que acabamos de armar
Out[41]:
In [42]:
using ScikitLearn
@sk_import linear_model: RidgeCV
Out[42]:
In [43]:
rcv = RidgeCV(alphas=0.01:0.01:10.0)
Out[43]:
In [44]:
ScikitLearn.fit!(rcv, X, Y)
Out[44]:
In [45]:
rcv[:coef_]
Out[45]:
In [46]:
rcv[:intercept_]
Out[46]:
In [47]:
α = rcv[:alpha_]
Out[47]:
Ridge Regression usando ScikitLearn.jl:
In [48]:
@sk_import linear_model: Ridge
Out[48]:
In [49]:
r = Ridge(alpha = α)
Out[49]:
In [50]:
ScikitLearn.fit!(r, X, Y)
Out[50]:
In [51]:
r[:coef_]
Out[51]:
In [52]:
rcv[:intercept_]
Out[52]:
Ridge Regression usando MultivariateStats:
In [53]:
MultivariateStats.ridge(X, Y, α)
Out[53]:
Ridge Regression usando R and RCall:
In [57]:
using RCall
R"""
library(MASS)
lmr <- lm.ridge($Y ~ $X, lambda=$α)
"""
Out[57]: