Métodos de Regresión.

3o Estadı́stica y Empresa

Práctica 8. Diagnosis y Validación de Modelos

Continuamos con el ejemplo de la práctica 1 cuyos datos están en el fichero practica1.RData.
En él se encuentran 180 observaciones de un paı́s europeo con información sobre la demanda
de vino de Rioja (miles de litros), Y , el precio medio en euros del vino, X1 , el precio medio en
euros de la cerveza, X2 , los gastos en publicidad en cientos de euros realizados por las bodegas,
X3 y un indicador del poder adquisitivo en el paı́s, X4 . Estima los parámetros del modelo:

Y = β0 + β1 X1 + β2 X2 + β3 X3 + β4 X4 + u

1. Comenzamos con los gráficos de los residuos frente a las variables explicativas y a los
valores ajustados yb, utilizando la función

residualPlots(RegModel.1)

2. Se puede realizar una primera diagnosis del modelo con la opción del menú:

Modelos → Gráficas → Gráficas básicas de diagnóstico

3. Linealidad. En cuanto a los gráficos especı́ficos de linealidad, tenemos:

a) El gráfico de regresión parcial o variables añadida se obtiene con la opción del menú:
Modelos → Gráficas → Gráficas de variables agregadas
b) El gráfico residual más componente o residual parcial se obtiene con la opción del
menú:
Modelos → Gráficas → Gráficas de componentes+residuos
c) Los gráficos de modelo marginal se obtienen con la función mmp() si es para una única
variable o con la función mmps() para obtener los de todas las variables. El código
para esta opción, incluyendo un gráfico sobre los valores ajustados, yb, es
mmps(RegModel.1)

Estas funciones se encuentran en el paquete car.

d ) Y el gráfico de respuesta inversa se obtiene con la función inverseResponsePlot()
del paquete car. Ası́, hay que escribir el siguiente código
inverseResponsePlot(RegModel.1)

e) Para realizar el test RESET, se utiliza la opción del menú:

Modelos → Diagnósticos numéricos → Test RESET de no linealidad...
En la ventana que aparece hay que indicar las potencias a incluir y el tipo de prueba
(sobre las variables explicativas, sobre los valores ajustados, que es la que selecciona-
mos, o sobre la primera componente principal). También se puede hacer directamente
con el código (cargando previamente el paquete lmtest):
resettest(RegModel.1)

1 Antonio Conde Sánchez

 f ) Si hubiera que tratar la no linealidad, habrı́a que comprobar si es necesario transfor-
mar las variables o no. Un análisis previo se puede llevar a cabo analizando la simetrı́a
y normalidad de las variables. Se puede hacer para todas las variables con la matriz
de diagramas de dispersión:
Gráficas → Matriz de diagramas de dispersión...
o bien con la siguente opción del menú, para cada una de las variables explicativas:
Gráficas → Estimar densidad...
Si sospechamos que hay que transformar alguna de las variables podemos utilizar
las transformaciones de Box-Cox. Para ello se utiliza la función powerTransform del
paquete car, la cual se puede llevar a cabo con la opción del menú:
Estadı́sticos → Resúmenes → Transformación para normalizar...
Mientras que para la transformación univariante de Y se utiliza, o bien el gráfico de
respuesta inversa (como en apartado 3d ), o la opción del menú :
Modelos → Diagnósticos numéricos → Transformación de la respuesta...
Otra opción para determinar la transformación Box-Cox de Y es a través de la fun-
ción boxCox() del paquete car, que proporciona un gráfico de la log-verosimilitud
frente al valor de λ, incluyendo un intervalo de confianza al 95 %. Si queremos que se
muestren los resultados numéricos en vez del gráfico se puede utilizar el argumento
plotit=FALSE. Tambiém se puede indicar el intervalo donde se mueve λ, que por
defecto va de -2 a 2. Por ejemplo, para que λ estuviera entre -1 y 4, el código es
boxCox(RegModel.1,lambda=seq(-1,4,0.1))

4. Normalidad. Se puede comprobar mediante gráficos (de densidad, QQ-plot,...) o tests de

normalidad.

a) Un QQ-plot aparece entre los gráficos básicos de diagnosis. También se pueden guar-
dar los residuos estandarizados y obtener con la opción del menú:
Gráficas → Gráfica de comparación de cuantiles...
b) Se puede realizar un gráfico de densidad de los mismos con la opción del menú:
Gráficas → Estimar densidad...
c) Los tests de normalidad (test de Shapiro-Wilk, Lilliefors, . . . ) sobre dichos residuos
se realizan con la opción del menú
Estadı́sticos → Resúmenes → Test de normalidad...
O directamente con el código
shapiro.test(rstandard(RegModel.1))

También es posible utilizar la función ols test normality() del paquete olsrr. Si la
aplicamos directamente sobre el modelo se realizarı́an sobre los residuos de mı́nimos
cuadrados, pero podemos hacerlo sobre los residuos estandarizados con el siguiente
código
ols_test_normality(rstandard(RegModel.1))

2 Antonio Conde Sánchez

5. Heterocedasticidad. Un gráfico adecuado se muestra en los gráficos básicos de diagnosis.

a) El contraste de Goldfeld-Quandt se puede realizar con la función gqtest() del pa-

quete lmtest. Por defecto se divide la muestra en dos subconjuntos y no se eliminan
observaciones. Se puede determinar el número de observaciones a eliminar con el
argumento fraction. También hay que indicar cuál es la variable que provoca la
heterocedasticidad con el argumento order.by. En concreto podemos escribir el si-
guiente código considerando la variable Renta:
gqtest(RegModel.1, fraction=0.25, order.by = ~ Renta, data = Datos)

b) El contraste de Breusch-Pagan se lleva a cabo con la opción del menú:

Modelos → Diagnósticos numéricos → Test de Breusch-Pagan para
heterocedasticidad...
Hay que indicar en la ventana que aparece si se hace el test estudentizado y qué
variables se utilizan (los regresores, los valores ajustados o funciones de las variables).
Se está utilizando la función bptest() del paquete lmtest, por lo que directamente
lo podemos obtener con el código (harı́a el test estudentizado):
bptest(RegModel.1)

c) El contraste de White se realiza de la misma forma que el test de Breusch-Pagan

estudentizado introduciendo las variables explicativas, sus cuadrados y sus productos
cuadrados. Para ello en el campo Otro (especificar) introducimos la fórmula:
(Gasto+Precio.cer+Precio.vin+Renta)^2+I(Gasto^2)+I(Precio.cer^2)+
I(Precio.vin^2)+I(Renta^2)
o bien
fitted(RegModel.1)+I(fitted(RegModel.1)^2)

d ) Para estimar por mı́nimos cuadrados ponderados se utiliza el argumento weights en

la función lm().
e) Otra opción es utilizar la función gls del paquete nlme. En concreto se puede utili-
zar el siguiente código, para estimar por MCGF si hay heterocedasticidad (hay más
opciones, como VarPower o VarConstPower):
library(nlme)
RegModel.1h <- gls(Demanda~Gasto+Precio.cer+Precio.vin+Renta, data=Datos,
weights = varExp(form=~fitted(.)))
summary(RegModel.1h)

6. Autocorrelación.

a) Los gráficos de autocorrelación y autocorrelación parcial se obtienen con las funciones

acf() y pacf(), respectivamente. En concreto hay que escribir el código
acf(RegModel.1$residuals)
pacf(RegModel.1$residuals)

3 Antonio Conde Sánchez

 b) El test de Durbin-Watson se realiza con la opción del menú:
Modelos → Diagnósticos numéricos → Test de Durbin-Watson para
autocorrelación...
c) Otro test que se puede realizar es el test de Breusch-Godfrey que permite comprobar
si hay autocorrelación de cualquier orden. Para ello se utiliza la función bgtest() del
paquete lmtest. Ası́, se puede realizar el test para orden 4 con el código
(bg<-bgtest(RegModel.1,order=4))
coeftest(bg)

Por defecto, si no se indica el argumento order, realiza el contraste para autocorre-

lación de orden 1.
d ) En caso de que exista autocorrelación de orden 1 se pueden estimar los parámetros
utilizando el procedimiento de Cochrane-Orcutt, el cual se puede llevar a cabo con la
función cochrane-orcutt del paquete orcutt. El código a utilizar es:
library(desk, pos=18)
cochorc(RegModel.1)

e) Otra opción es utilizar la función gls del paquete nlme. En concreto hay que utilizar
el siguiente código para estimar por MCGF si hay autocorrelación de orden 1:
library(nlme)
RegModel.1a <- gls(Demanda~Gasto+Precio.cer+Precio.vin+Renta, data=Datos,
correlation = corAR1(), method="ML")
summary(RegModel.1a)

7. Errores estándar robustos. Se pueden obtener con la opción del menú

Modelos → Resumir el modelo

y marcando la casilla para Usar el estimador sandwich.... Se pueden seleccionar diferentes

tipos de estimadores.

8. Validación. Como no tenemos dos subconjuntos se utiliza el método jacknife. Para calcu-
lar el error cuadrático medio y el R2 se escribe el código, utilizando las funciones ols press
y ols pred rsq del paquete olsrr:

(press<-ols_press(RegModel.1))
df<-RegModel.1$df.residual
(mse_jack<-press/df)
(mse<-summary(RegModel.1)$sigma^2)
(r2_pred<-ols_pred_rsq(RegModel.1))
(r2<-summary(RegModel.1)$r.squared)

4 Antonio Conde Sánchez