Regresión [Regression]

Iremos completando y actualizando este amplio tema con el tiempo (no dejes de revisar la página!).Pincha aquí para acceder al archivo de regresión lineal con cuestiones teóricas, resumen sobre los contraste de ajuste en el modelo y análisis de datos.Para más información mira aquí (regresión lineal simple) y aquí (script R ANOVA).

Objetivo: modelar asociaciones entre variables.
El análisis de regresión se utiliza para modelar la relación entre una variable Y (v.respuesta o dependiente) y una o más variables x1, x2,..., xp-1 (v.predictora o independiente). La v.respuesta debe ser continua, pero las v.predictoras pueden ser contínuas, discretas o categóricas.
En su forma más simple, la regresión describe una relación lineal entre una variable predictora o independiente (eje-x) y una variable respuesta o dependiente (eje-y). Para ajustar una regersión lineal a los datos, se utiliza el método de mínimos cuadrados y a fin de evaluar los parámetros del modelo ajustado se usan los test de hipótesis.
Se deben tener en cuenta los supuestos del modelo y los test diagnósticos para evaluar el ajuste de los datos al modelo y las predicciones realizadas con él.
Como tópicos avanzados discutiremos: regresión logística, regersión múltiple, regersión no-lineal, regresión robusta y análisis de vías (path analysis).

Índice

• Regresión lineal simple
• Regresión lineal múltiple
• Mínimos cuadrados ordinarios (LSQ)
• Máxima probabilidad (ML)
• Test de hipótesis para los parámetros β
• Aproximación ANOVA a la regresión
• Validación y selección del modelo adecuado
• Diagnóstico
• Transformaciones
• Colinealidad

1. Regresión lineal simple

Y_i=β_0 + β_1*x_i + ε_i

donde Y_i es el valor de la variable respuesta para la observación i, β_0 y β_1 son parámetros, x_i es una constante conocida para la observación i y ε_i es un término de error con distribución N(0,σ^2) con varianza σ^2 desconocida.
Si tenemos n observaciones podemos expresar el sistema de ecuaciones en notación matricial: Y=X*β+ε, donde Y es un vector nx1, X una matriz nx2, β un vector 2x1 y ε un vector nx1 con ε~N(0, σ^2I) donde σ^2I es la matriz de varianza-covarianza del vector de errores.

Ejemplo: veremos una serie de scripts en R, muy sencillos.Estudiamos la relación entre grade point average (GPA, de estudiantes de primer sementre en una Universidad) y scholastic aptitude test (sat, de estos estudiantes).

Graficamos las variables y estimamos por mínimos cuadrados los valores de beta0 y beta1 (mediante la función lm()).

attach(Grades)                                                       
plot(sat,gpa) #el gráfico sugiere relación lineal entre gpa y sat     
Y<-gpa; x<-sat                                                         
model<-lm(Y~x); model$coefficients #LSQ utilizando la función lm()

Created by Pretty R at inside-R.org

Calculamos la matriz de varianza-covarianza de beta, probamos si existe una relación lineal a nivel alfa=.10 (test de hipótesis) y construimos intervalos de confianza al 90% para beta0 y beta1.

simple.model<-lm(Y~x)       #cálculo directo de la matriz var-cov vcoc(simple.model)            #el test estadístico estandarizado tiene distribución t198 y H1 es una hipótesis two-sided, por lo cual la región de rechazo es |tobs|>t.95,198. Aquí tobs=beta/SEbeta                                                                    
tobs<-summary(simple.model)$coef[2,3]          #intervalo de confianza al 90% para beta0 y beta1                                                                       
confint(simple.model,level=.9)

Construimos la tabla ANOVA y probamos si existe una relación lineal entre las variables a nivel alfa=.05.
anova(simple.model)




2. Regresión lineal múltiple

Y_i=β_0+β_1 x_i1+⋯+β_((p-1)) x_(i(p-1))+ε_i, i=1,…,n

Que podemos expresar en notación matricial como Y=Xβ+ε, donde Y es un vector nx1, X una matriz nxp, β un vector px1 y ε un vector nx1 con ε~N(0,σ^2 I) donde σ^2 I es la matriz de varianza-covarianza del vector de errores. Cuando se asume ε~N(0,σ^2 I), el modelo se llama modelo de error normal y se tiene que Y~N(Xβ,σ^2 I).
El vector β es un vector de constantes desconocidas que son estimadas a partir de los datos. Cada β_j j=1,…,p-1 indica el cambio en E[Y|x_ij] para un i fijo cuando x_ij aumenta una unidad y los demás predictores permanecen constantes.
Asumiendo que no existe error de medición en los valores de x_ij, podemos estimar los parámetros β_j mediante uno de las dos técnicas siguientes: mínimos cuadrados ordinarios (ordinary lest squares LSQ) o método de máxima probabilidad (máximum likelihood ML).


3. Mínimos cuadrados ordinarios (LSQ)Este criterio minimiza la suma de las desviaciones de los cuadrados de los Y_i respecto a los valores esperados. La desviación i (error) es: ε_i=Y_i-E(Y_i ), y las estimaciones de los β_j se calculan minimizando la cantidad Q: Q=∑ε_i^2 =ϵ^' ϵ=(Y-Xβ)^' (Y-Xβ).
Para el modelo de regresión lineal simple: ε_i=Y_i-(β_0+β_1 x_i). Resolviendo para (β_0 ) ̂ y (β_1 ) ̂, tenemos:
(β_0 ) ̂=Y ̅-(β_1 ) ̂x ̅ y (β_1 ) ̂=(∑(x_i-x ̅ ) (Y_i-Y))/(∑(x_i-x ̅ )^2 )
Cuando calculamos los valores β, la línea de regresión ajustada es (Y_i ) ̂=(β_0 ) ̂+(β_1 ) ̂x_i, donde los errores estimados (predichos), llamados residuos, son: (ε_i ) ̂=Y_i-(Y_(i.) ) ̂

Ejercicio en R
Investigamos la relación entre la puntuación media de estudiantes en el primer semestre de una universidad (GPA) y su puntuación de aptitud escolar (SPA).

attach(Grades)                  
plot(sat,gpa)                                                                
Y<-gpa;x<-sat
model<-lm(Y~x)                                            
model$coeff

4. Máxima probabilidad (ML)
Otro método común para estimar β es utilizar el estimador de máxima probabilidad (MLE) de β cuando ε~N(0,σ^2 I) mediante la construcción de una función de probabilidad. Sea la función probabilidad para β y σ^2 cuando X es dado:
L(β,σ^2│X)=∏[1/√(2πσ^2 )]* exp⁡[-(∑ε_i^2 )/(2σ^2 ) ]=∏[1/√(2πσ^2 )]* exp⁡[-((Y-Xβ)^' (Y-Xβ))/(2σ^2 ) ]
entonces la función log-likelihood es lnL(β,σ^2│X), entonces β ̂=(X'X)^(-1) X^' Y, (σ^2 ) ̂=((Y-Xβ ̂ )^' (Y-Xβ ̂ ))/n. Aunque el estimador (σ^2 ) ̂ es un estimador sesgado de σ^2.

5. Test de hipótesis para los parámetros β
Sea Ho:β_k=β_k0 vs. H1: β_k =! β_k0, utilizamos un t-estadístico estándar:
(Estimador insesgado-valor hipotetizado)/(error estándar del estimador)
esto es: (β ̂_k-β_k0)/s_(β_k ) ~t_(n-p) para k=0,…,p-1, donde s_(β_k )=√(matriz varianza-covarianza)
Sea Ho:β_k=0 vs. H1: β_k =! 0, utilizamos un t-estadístico estándar: t_obs=β ̂_k/s_(β_k ) , donde el p-valor está dado por p-value=2*P(t_(n-p)>|t_obs |) y el intervalo de confianza al (1-alfa)*100% para β_k es:
CI_(1-α) (β_k )=[β ̂_k-t_(1-α/2;n-p)*s_(β_k ),β ̂_k+t_(1-α/2;n-p)*s_(β_k ) ].

Ejemplo en R

var.cov.b<-vcov(model)                                 
tobs<-summary(model)$coef[2,3:4]  #sea Ho:β_k=0 vs. H1: β_k=!0                     
# entonces tobs=15.912, Pr(>|t|)=2.92-37, |tobs|>t.95; 198=1.6526, refutamos la hipótesis nula                                                         
confint(model, level=.9) #sea alfa=.1

6. Aproximación ANOVA a la regresión
El mismo modelo de la regresión con error normal puede ser considerado en la herramienta de análisis de varianza. El análisis de varianza se basa en la partición de la suma de cuadrados y los grados de libertad asociados con la variable respuesta Y. La desviación total Y_i-Y ̅ puede descomponerse en dos componentes: 1) la desviación de los valores ajustados alrededor de la media y 2) la desviación de las observaciones alrededor de la línea de regresión.
Esto es, la suma total de cuadrados (SST) es la suma de la suma de cuadrados de la regresión (SSR) y la suma de cuadrados del error (residual, SSE):
SST=SSR+SSE
∑(Y_i-Y ̅ )^2 =∑(Y ̂_i-Y ̅ )^2 +∑(Y_i-Y ̂_i )^2
ANOVA con regresión lineal simple

Tabla 1. Tabla de ANOVA completa para la regresión lineal con un único factor.

Fuente g.l. SS MS E(MS) F P-valor
Regresió n 1 SSreg SSreg/1 σ^2+β^2_1+sum(X^2) E(MSreg)/E(MSresid) F con 1,n-2
Residual n-2 RSS RRS/(n-2) σ^2
Total n-1 SSY SSY/(n-1) σ^2*Y

Cuando dividimos la suma de cuadrados por sus grados de libertad asociados, obtenemos los cuadrados medios (MS). En el caso de la regresión lineal simple:
SSR/1=MSR y SSE/(n-1)=MSE. Además E[MSR]=σ^2+β_1^2 ∑(x_i-x ̅ )^2 y E[MSE]=σ^2.
Sea Ho:β_k=0 vs. H1: β_k =! 0, la prueba estadística es: F_obs=MSR/MSE y MSR/MSR~F_(1,n-2). Valores del estadístico cercanos a 1 tienden a afirmar la hipótesis nula, mientras que valores grandes del estadístico apoyan la hipótesis alternativa. Específicamente, refutamos la hipótesis nula si F_obs>f_(1-α;1,n-2).

Ejemplo en R

anova(model)                                                                       
qf(.95, 1, 198) #f_(.95;1,198)
#se observa que debemos refutar la hipótesis nula, y existe fuerte evidencia para sugerir la relación lineal entre gpa y sat.

ANOVA con regresión lineal múltiple

Tabla 2. ANOVA para la regresión lineal múltiple

Fuente g.l. SS MS F P-valor
Regresión k SSR=(β´) ̂X´Y-(1/n)Y´JY MSR=SSR/(k+1) MSR/MSE F alfa,k+1,n-k-1.
Residual n-(k+1) SSE=Y´Y-(β ̂´) ̂X´Y MSE=SSE/(n-k-1)
Total n-1 SST=Y´Y-(1/n)Y´JY

La tabla de ANOVA para el análisis de regresión múltiple puede generalizarse a siguiente test de hipótesis: Ho:β_1=⋯=β_(p-1) vs. H1: al menos un β_i =! 0,i=1,…,1-p

Ejemplo en R
Analizamos 9 variables para un grupo de 78 luchadores de estudiantes de secundaria.

attach(HSwrestler)                                                             
Y<-HWFAT; n<-78                                                             
x1<-ABS, x2<-TRICEPS; x3<-SUBSCAP; X<-cbind(rep(1,n),x1,x2,x3) 
mod123<-lm(Y~x1+x2+x3) 
anova(mod123)                                                      
anova(mod123)[3,4] #Fobs                                           
summary(mod123)$coeff[4,3]^2 #indica el tobs al cuadrado, el cual es equivalente al Fobs

7. Validación y selección del modelo adecuado
Procedimientos basados en pruebas
Selección del mejor subconjunto de predictores mediante: 1) una estrategia stepwise que compara modelos sucesivos y 2) una aproximación de criterio que intenta maximizar alguna medida de bondad de ajuste.
Stepwise: La regresión stepwise es una combinación de la eliminación backward y la selección forward:
Selección backward
Comienza con un modelo que contiene todas las variables xs potenciales e identifica aquellas con mayor p-valor. El alfa-crítico se llama también el “p a remover” y típicamente es un valor entre 15 o 20%.
Selección forward
Comienza sin variables en el modelo y luego agrega la variable x que produce el menor p-valor debajo del alfa-crítico cuando se incluye en el modelo.

Criterio: por ejemplo R2-ajustado, Cp de Mallows, Criterio de información de Bayes (BIC), criterio de información de Akaike (AIC), et.

Ejemplo en R

library(leaps)                                                                    #criterio R2-ajustado  
a<-regsubsets(as.matrix(HSwrestler[,-c(7,8,9)]),HSwrestler[,7])                
summary(a); summary(a)$adjr2                                                #criterio Cp de Mallows summary(a)$cp                                                               #criterio AIC                                            library(MASS)                          
reg.all<-lm(HWFAT~AGE+HT+ABS+TRICEPS+SUBSCAP) 
mod.aic<-stepAIC(reg.all,direction=”both”,scope=(~ AGE+HT+ABS+TRICEPS+SUBSCAP), k=2) #criterio BIC                      
mod.bic<-stepAIC(reg.all,direction=”both”,scope=(~ AGE+HT+ABS+TRICEPS+SUBSCAP), k=log(length(HWFAT)))

8. Diagnóstico
A pesar de que ajustar un modelo mediante el principio de mínimos cuadrados no requiere supuestos de distribución, utilizar éste modelo para realizar inferencias depende de supuestos específicos. Si asumimos que ε~N(0,σ^2 I) debemos verificar el supuesto de errores independientes y con distribución normal de media cero y varianza constante.

Ejemplo en R: Chequear los supuestos del error

mod.3<-lm(HWFAT~AGE+ABS+TRICEPS)                                   
library(lmtest)
bptest(mod.3); 
qqnorm(mod.3) #evalua normalidad y var cte. 
durbin.watson(mod.3) #evalúa autocorrelación             
shapiro.test(resid(mod.3)) #evalúa normalidad                             
plot(mod.3)

Ejemplo en R: Identificar observaciones inusuales

library(MASS)                                
par(mfrow=c(2,2))                                  
plot(fitted(mod.3),resid(mod.3),title=”Residuals vs Fitted”); abline(h=0) plot(fitted(mod.3),stdres(mod.3),title=”Standardized Residuals vs Fitted”); abline(h=0)                                 plot(fitted(mod.3),studres(mod.3),title=”Studentized Residuals vs Fitted”); abline(h=0)                                                    plot(mod.3,main=”Default Graph 1”)                                  
par(mfrow=c(1,1))

Tabla 3. Resumen de las medidas de observaciones influyentes

Medidas de influencia Fórmula Caso i puede ser influyente si:
Cook Di (r_i^2)/p(h_ii/(1-h_ii )) D_i>f_(.5;p,n-p)
DFFITS r_i^* (h_ii/(1-h_ii )) |DFFITS|>2√(p/n)
DFBETAS ((β_k ) ̂-(β_k(i) ) ̂)/√((σ_((i))^2 ) ̂v_(k+1,k+1) ) |DFBETAS|>2/√n

Ejemplo en R: Chequear los supuestos del error

X<-model.matrix(mod.3) #definir p y n                                    
n<-nrow(X); p<-ncol(X)                                  
hii<-lm.influence(mod.3)$hat #valores hat                                        
plot(hii, ylab=”Leverage”)                                               
cv<-2*p/n; abline(h=cv,tly=2)

plot(lm.influence(mod.3)$coeff[,2],ylab=”Difference in Coefficients”)                                                         
plot(dfbetas(mod.3)[,2],ylab=”DFBETAS”)                               
cv<-2/sqrt(n); abline(h=c(-cv,cv))                      
plot(studres(mod.3),ylab=”Studentized Residuals”)                  
cv<-qt(1-.1/(2*n),n-p-1); abline(h=c(-cv,cv))                                   

DFFITS<-studres(mod.3)*(hii/(1-hii))^.5 
plot(DFFITS,ylab=”DFFITS”)                                              
cv<-2*sqrt(p/n); abline(h=c(-cv,cv))     
cd<-cookd(mod.3)                                                
plot(cd,ylab=”Cook´s Distance”)                                            
CF<-qf(.5,p,n-p); abline(h=CF)                                                   par(mfrow=c(1,1))

9. Transformaciones
Cuando el análisis de residues revela serios problemas, o cuando las relaciones entre la respuesta y los predictors no es lineal, la regresión puede aún producir un modelo razonable si transformamos la variable respuesta, los predictores o ambos. Luego de la transformación del predictor se debe volver a analizar los supuestos de normalidad.

Ejemplo en R

#transformación raíz cuadrada 
par(mfrow=c(2,3))                                                          
plot(x1,Y); lines(x1,x1^.5) #scatterplot                               
plot(lm(Y~x1),which=c(1,2));  #residuals vs. fitted values and quantile-quantile de residues estandarizados                                                 
plot(x1^.5,Y); abline(lm(Y~I(x1^.5))) #scatterplot      
plot(lm(Y~I(x1^.5)),which=c(1,2)) #residuals vs. fitted values and quantile-quantile de residues estandarizados                                          
par(mfrow=c(1,1))

#transformación cuadrática 
par(mfrow=c(2,3))                                                        
plot(x2,Y); lines(x2,x2^2) #scatterplot                                          plot(lm(Y~x2),which=c(1,2))  #residuals vs. fitted values and quantile-quantile de residues estandarizados                                                       

plot(x2^2,Y); abline(lm(Y~I(x2^2))) #scatterplot                     
plot(lm(Y~I(x2^2)),which=c(1,2)) #residuals vs. fitted values and quantile-quantile de residues estandarizados                                           
par(mfrow=c(1,1))

#transformación recíproca 
par(mfrow=c(2,3))                                                        
plot(x3,Y); lines(x3,x3^(-1)) #scatterplot                               
plot(lm(Y~x3),which=c(1,2));  #residuals vs. fitted values and quantile-quantile de residues estandarizados                                                  
plot(x3^2,Y); abline(lm(Y~I(x3^(-1)))) #scatterplot
plot(lm(Y~I(x3^(-1))),which=c(1,2)) #residuals vs. fitted values and quantile-quantile de residues estandarizados                                  
par(mfrow=c(1,1))

10. Colinealidad
Esto ocurre cuando algunos de los predictores es combinación lineal de otros predictores. La multicolinealidad causa problemas en la estimación de los betas y sus interpretaciones. Para detectar la colinealidad calculamos el número condición y el factor de inflación de la varianza.
Un valor del número condicional k entre 30-100 indica que existen dependencias de moderadas a fuertes entre los predictores. K>100 indican serios problemas de multicolinealidad.
Cuando existen dependencia entre los predictores R2-ajustado será cercana a 1 y VIF será grande. VIF>10 sugiere una colinealidad seria.

Ejemplo en R

attach(SimDataXT)                                              
modC<-lm(Y~I(x1^.5)+ I(x1^2))                                                  
summary(modC)                                                                   
X<- model.matrix(modC)                                                 
kappa(X,exact=TRUE)   #estima el número condición de una matriz                  
vif(modC) #calcula los factores de inflación de la varianza

Otros: Transformaciones para errores con distribuciones que no son normales y varianza desigual: Transformación box-cox.

BIBLIOGRAFÍA

• Gotelli & Ellison. 2004. A Primer of Ecological Statistics.
• Ugarte et al. 2008. Probability and Statistics with R.

PRINCIPALES ÓRDENES PARA EL ANÁLISIS DE REGRESIÓN EN R

#########################################################
### Funciones de R para el análisis de regresión
#########################################################
 
# 0) Análisis preliminar de los datos
  data()
  names()
  summary()
  dotchart()
  pairs()
  boxplot()
  xyplot()
 
# 1) Ajuste del modelo
    # Regresión lineal simple
    fit <- lm(y ~ x1 , data=mydata)
    # Regresión lineal múltiple
    fit <- lm(y ~ x1 + x2 , data=mydata)
    summary(fit) # resultados
 
    coefficients(fit) # coeficientes del modelo. También coef(). library(stats) 
    coeftest(fit) # test de los coeficientes estimados. library(lmtest)
    confint(fit, level=0.95) # intervalos de confianza para los parámetros del modelo. library(stats)
    fitted(fit) # valores predichos
    residuals(fit) # residuos
    anova(fit) # anova. library(stats)
    vcov(fit) # matriz de varianza-covarianza para los parámetros del modelo
    influence(fit) # diagnóstico   
 
# 2) Gráficos de diagnóstico
    #Revisar heteroscedasticidad, normalidad, y observaciones extremas.
    layout(matrix(c(1,2,3,4),2,2)) # para graficar todo en una misma hoja
    plot(fit)
 
    # Evaluar outliers.
    outlierTest(fit) # p-valor de Bonferonni para las observaciones extremas. library(car)
    qqPlot(fit, main="QQ Plot") #qq plot para residuos estudentizados
    leveragePlots(fit) # gráfico de leverage 
 
    # Observaciones influyentes
    influence(fit)
    # agregamos gráficos
    av.Plots(fit)
    # distancia D de Cook: identificar valores D > 4/(n-k-1)
    cutoff <- 4/((nrow(mtcars)-length(fit$coefficients)-2))
    plot(fit, which=4, cook.levels=cutoff)
    # Gráficos de influencia
    influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportial to Cook's Distance" )
    # Otra forma de realizar gráficos:
    library(MASS)                                
    par(mfrow=c(2,2))                                  
    plot(fitted(fit),resid(fit),title=”Residuals vs Fitted”)
    abline(h=0) 
    plot(fitted(fit),stdres(fit),title=”Standardized Residuals vs Fitted”)
    abline(h=0)                                 
    plot(fitted(fit),studres(fit),title=”Studentized Residuals vs Fitted”) 
    abline(h=0)                                                    
    plot(fit,main=”Default Graph 1”)                                  
    par(mfrow=c(1,1))
 
    # Normalidad de los residuales
    # qq plot para residuos esstudentizados
    qqPlot(fit, main="QQ Plot")
    # distribución de los residuos esstudentizados
    library(MASS)
    sresid <- studres(fit)
    hist(sresid, freq=FALSE, main="Distribution of Studentized Residuals")
    xfit<-seq(min(sresid),max(sresid),length=40)
    yfit<-dnorm(xfit)
    lines(xfit, yfit)
 
    # Homocedasticidad
    # Test de la varianza de los errores no constante
    ncvTest(fit)
    # plot residuos estudentizados vs. valores ajustados
    spreadLevelPlot(fit)
    # Test de Breusch-Pagan
    library(lmtest)
    bptest(fit)
 
    # Colinealidad
    vif(fit) # factores de inflación de la varianza
    sqrt(vif(fit)) > 2 # problem?
    # Otra forma:
    X<- model.matrix(fit)                                                 
    kappa(X,exact=TRUE) #estima el número condición de una matriz                       
    vif(fit) #calcula los factores de inflación de la varianza
 
    # Nolinealidad
    # componente + gráfico de residuos 
    crPlots(fit)
    # Gráficos Ceres 
    ceresPlots(fit)
    # Test de Shapiro
    shapiro.test(resid(fit))
 
    # Errores autocorrelacionados
    durbinWatsonTest(fit)
    durbin.watson(fit)
    #Otra forma de chequear los supuestos del error:
      X<-model.matrix(fit) #definir p y n                                    
      n<-nrow(X); p<-ncol(X)                                  
      hii<-lm.influence(fit)$hat #valores hat                                        
      plot(hii, ylab=”Leverage”)                                               
      cv<-2*p/n; abline(h=cv,tly=2); 
      plot(lm.influence(fit)$coeff[,2],ylab=”Difference in Coefficients”)                                                         
      plot(dfbetas(fit)[,2],ylab=”DFBETAS”)                               
      cv<-2/sqrt(n)
      abline(h=c(-cv,cv))                      
      plot(studres(fit),ylab=”Studentized Residuals”)                  
      cv<-qt(1-.1/(2*n),n-p-1)
      abline(h=c(-cv,cv))                                   
      DFFITS<-studres(fit)*(hii/(1-hii))^.5 
      plot(DFFITS,ylab=”DFFITS”)                                              
      cv<-2*sqrt(p/n)
      abline(h=c(-cv,cv))     
      cd<-cookd(fit)
      plot(cd,ylab=”Cook´s Distance”)                                            
      CF<-qf(.5,p,n-p) >abline(h=CF)                                                       
      par(mfrow=c(1,1))
 
 
    # Test grlobal de los supuestos del modelo
    library(gvlma)
    gvmodel <- gvlma(fit)
    summary(gvmodel)
 
# 3) Comparar modelos
    #You can compare nested models with the anova( ) function. The following code provides a simultaneous test that x3 and x4 add to linear prediction above and beyond x1 and x2.
    # compare models
    fit1 <- lm(y ~ x1 + x2 + x3 + x4, data=mydata)
    fit2 <- lm(y ~ x1 + x2)
    anova(fit1, fit2)
 
# 4) Validación cruzada
    #You can do K-Fold cross-validation using the cv.lm( ) function in the DAAG package.
    # K-fold cross-validation
    library(DAAG)
    cv.lm(df=mydata, fit, m=3) # 3 fold cross-validation
 
    #You can assess R2 shrinkage via K-fold cross-validation. Using the crossval() function from the bootstrap package, do the following:
    # Assessing R2 shrinkage using 10-Fold Cross-Validation
    fit <- lm(y~x1+x2+x3,data=mydata)
    library(bootstrap)
    # define functions
    theta.fit <- function(x,y){lsfit(x,y)}
    theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef}
    # matrix of predictors
    X <- as.matrix(mydata[c("x1","x2","x3")])
    # vector of predicted values
    y <- as.matrix(mydata[c("y")])
    results <- crossval(X,y,theta.fit,theta.predict,ngroup=10)
    cor(y, fit$fitted.values)**2 # raw R2
    cor(y,results$cv.fit)**2 # cross-validated R2
 
# 5) Selección de variables
    # Selecting a subset of predictor variables from a larger set (e.g., stepwise selection) is a controversial topic. You can perform stepwise selection (forward, backward, both) using the stepAIC( ) function from the MASS package. stepAIC( ) performs stepwise model selection by exact AIC.
 
    # Stepwise Regression
    library(MASS)
    fit <- lm(y~x1+x2+x3,data=mydata)
    step <- stepAIC(fit, direction="both")
    step$anova # display results
 
    # Otras formas:
      library(leaps)                                                                    
    #criterio R2-ajustado  
      a<-regsubsets(as.matrix(HSwrestler[,-c(7,8,9)]),HSwrestler[,7]); summary(a); summary(a)$adjr2                                                
    #criterio Cp de Mallows 
      summary(a)$cp                                                               
    #criterio AIC                                             
      library(MASS)                          
      reg.all<-lm(HWFAT~AGE+HT+ABS+TRICEPS+SUBSCAP) 
      mod.aic<-stepAIC(reg.all,direction=”both”,scope=(~ AGE+HT+ABS+TRICEPS+SUBSCAP), k=2) 
    #criterio BIC                      
      mod.bic<-stepAIC(reg.all,direction=”both”,scope=(~ AGE+HT+ABS+TRICEPS+SUBSCAP), k=log(length(HWFAT)))
 
    Alternatively, you can perform all-subsets regression using the leaps( ) function from the leaps package. In the following code nbest indicates the number of subsets of each size to report. Here, the ten best models will be reported for each subset size (1 predictor, 2 predictors, etc.).
 
    # All Subsets Regression
    library(leaps)
    attach(mydata)
    leaps<-regsubsets(y~x1+x2+x3+x4,data=mydata,nbest=10)
    # view results
    summary(leaps)
    # plot a table of models showing variables in each model.
    # models are ordered by the selection statistic.
    plot(leaps,scale="r2")
    # plot statistic by subset size
    library(car)
    subsets(leaps, statistic="rsq")
 
# 6) Importancia relativa
    #The relaimpo package provides measures of relative importance for each of the predictors in the model. See help(calc.relimp) for details on the four measures of relative importance provided.
 
    # Calculate Relative Importance for Each Predictor
    library(relaimpo)
    calc.relimp(fit,type=c("lmg","last","first","pratt"),rela=TRUE)
 
    # Bootstrap Measures of Relative Importance (1000 samples)
    boot <- boot.relimp(fit, b = 1000, type = c("lmg","last", "first", "pratt"), rank = TRUE,diff = TRUE, rela = TRUE)
    booteval.relimp(boot) # print result
    plot(booteval.relimp(boot,sort=TRUE)) # plot result
 
# 7) Transformaciones     
    #transformación raíz cuadrada 
    par(mfrow=c(2,3))                                                          
    plot(x1,Y); lines(x1,x1^.5) #scatterplot                               
    plot(lm(Y~x1),which=c(1,2));  #residuals vs. fitted values and quantile-quantile de residues estandarizados 
    plot(x1^.5,Y); abline(lm(Y~I(x1^.5))) #scatterplot      
    plot(lm(Y~I(x1^.5)),which=c(1,2)) #residuals vs. fitted values and quantile-quantile de residues estandarizados                                          
    par(mfrow=c(1,1))
 
    #transformación cuadrática 
    par(mfrow=c(2,3))                                                        
    plot(x2,Y); lines(x2,x2^2) #scatterplot                                                                
    plot(lm(Y~x2),which=c(1,2));  #residuals vs. fitted values and quantile-quantile de residues estandarizados 
    plot(x2^2,Y); abline(lm(Y~I(x2^2))) #scatterplot                     
    plot(lm(Y~I(x2^2)),which=c(1,2)) #residuals vs. fitted values and quantile-quantile de residues estandarizados                                           
    par(mfrow=c(1,1))
 
    #transformación recíproca 
    par(mfrow=c(2,3))                                                        
    plot(x3,Y); lines(x3,x3^(-1)) #scatterplot                               
    plot(lm(Y~x3),which=c(1,2));  #residuals vs. fitted values and quantile-quantile de residues estandarizados                                                  
    plot(x3^2,Y); abline(lm(Y~I(x3^(-1)))) #scatterplot 
    plot(lm(Y~I(x3^(-1))),which=c(1,2)) #residuals vs. fitted values and quantile-quantile de residues estandarizados                                  
    par(mfrow=c(1,1))
 
    #Otros: Transformaciones para errores con distribuciones que no son normales y varianza desigual: Transformación box-cox.
 
# 8) Centrado y estandarizado de las variables (para simplificar la interpretación del modelo) .     #La transformación lineal de los predictores no afecta el ajuste del modelo de regresión.  
    # Centrado mediante la resta de la media
    cX<-X-mean(X) # no cambia la desviación estándar de los residuos ni el R^2, y el coeficiente y error estándar de la interacción no cambia, pero los efectos principales y el intercepto sí se modifican y ahora pueden ser interpretados comparando con la media de los datos.
    # Centrado utilizando un valor estándar de referencia
    c2X<-X-.5 # por ejemplo .5 es la diferencia media predicha entre un indivuo con 1 y 0 meses.
    # Estandarizado mediante la resta de la media y la división por 2 desviaciones estándar
    zX<-(X-mean(X))/(2*sd(X)) # permite interpretar los coeficientes en una escala común, excepto para el intercepto que ahora corresponde a la media predicha cuando todas las variables predictoras tienen valores medios.