Métodos Estadísticos II

(Grado en Ciencia e Ingeniería de Datos. Primer curso)

Curso 23/24 Práctica 7: Modelos lineales: Regresión
[Ejercicios propuestos (y soluciones)]

Contenidos
Modelos lineales: Regresión lineal simple 1
Ejercicio 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Ejercicio 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Ejercicio 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Ejericio 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Ejercicio 5 (caso de variable explicativa dicotómica). . . . . . . . . . . . . . . . . 21
Ejercicio 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Regresión lineal múltiple 27

Ejercicio 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Regresión logística 31
Ejercicio 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Ejercicio 9 (propuesto). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

Modelos lineales: Regresión lineal simple

Ejercicio 1.

Un investigador agrícola está analizando la relación entre la producción de papas ( y) y el nivel

de fertilizantes utilizado ( x). Para ello, ha dividido la extensión agrícola donde se desarrolla
el estudio en 8 partes iguales y ha aplicado diferentes niveles de fertilizante en cada una. El
nivel de fertilizante y la producción de cada zona viene dada en la siguiente tabla

1
 Nivel de fertilizante Producción
x y
1 25
1.5 31
2 27
2.5 28
3 36
3.5 35
4 32
4.5 34
Se pide:
1. Calcular los coeficientes de la recta de regresión por mínimos cuadrados y representar
gráficamente los puntos y dicha recta.
2. Calcular la varianza estimada para la recta mínimo cuadrática.
3. Supongamos que conocemos que la producción y dado un nivel de fertilizante x es normal
con media α + β × x y varianza σ 2 = 3.02 conocida. Utilizar una a priori normal con
media 2 y varianza 22 para β y obtener su distribución a posteriori.
4. Calcular un intervalo de credibilidad al 95% para β.
5. Desarrollar el test

H0 : β ≤ 0 vs H1 : β > 0

Solución.
1 y 2. Calculemos las cantidades necesarias con R.
x=c(1,1.5,2,2.5,3,3.5,4,4.5)
y=c(25,31,27,28,36,35,32,34)
n=length(x)
round(mean(x),4)

## [1] 2.75
round(mean(y),4)

## [1] 31
# Calculamos los SS
SSxy=sum((x-mean(x))*(y-mean(y)))
SSx=sum((x-mean(x))ˆ2)
SSy=sum((y-mean(y))ˆ2)
# Cálculo de los parámetros MCO
betahat=SSxy/SSx
alfahat=mean(y)-mean(x)*betahat
coeficientes=c(alfahat,betahat)
round(coeficientes,4)

2
## [1] 24.4524 2.3810
#Varianza estimada del ajuste
sigma2=sum((y-(alfahat+betahat*x))ˆ2)/(n-2)
round(sigma2,4)

## [1] 8.746
#Error residual estándar
round(sqrt(sigma2),4)

## [1] 2.9574
# Gráfica del ajuste
plot(x,y,pch=20, col="Blue",ylab="Producción",xlab="Nivel fertilizantes")
abline(alfahat,betahat, col="Red")
36
34
32
Producción

30
28
26

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5

Nivel fertilizantes
# Coeficiente de determinación (R2)
r2=SSxyˆ2/(SSx*SSy)
round(r2,4)

## [1] 0.5315
Es decir, la recta de regresión mínimo cuadrática es y = 24.45 + 2.38 × x. La interpretación
de β̃ = 2.38 es la siguiente: por cada unidad de incremento en fertilizantes, la producción
aumenta 2 − 38 unidades. La producción constante para cualquier cantidad de fertilizante es
de 24.45 unidades.
La regresión también puede obtenerse utilizando el comando lm de R.

3
x=c(1,1.5,2,2.5,3,3.5,4,4.5)
y=c(25,31,27,28,36,35,32,34)
recta=lm(y~x)
summary(recta)

##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.405 -2.036 -1.500 2.405 4.405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.4524 2.7189 8.993 0.000106 ***
## x 2.3810 0.9127 2.609 0.040185 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.957 on 6 degrees of freedom
## Multiple R-squared: 0.5315, Adjusted R-squared: 0.4534
## F-statistic: 6.806 on 1 and 6 DF, p-value: 0.04019
# Coeficiente de regresión
recta$coefficients

## (Intercept) x
## 24.452381 2.380952
# Residual standard error
res=summary(recta)
sigma=res$sigma
round(sigma,4)

## [1] 2.9574
# Varianza estimada para la recta
sigma2=res$sigmaˆ2
round(sigma2,4)

## [1] 8.746
# Coeficiente de deteminación R2
r2=res$[Link]
round(r2,4)

## [1] 0.5315

4
# Gráfica
plot(x,y,pch=20, col="blue",ylab="Producción",xlab="Nivel fertilizantes")
abline(recta, col="red")
36
34
32
Producción

30
28
26

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5

Nivel fertilizantes

Como hemos podido observar, se han obtenido las mismas cantidades que las obtenidas “a
mano”.
3. Los parámetros y la gráfica de la densidad a posteriori de β serán:
# datos
x=c(1,1.5,2,2.5,3,3.5,4,4.5)
y=c(25,31,27,28,36,35,32,34)
sigma2=3ˆ2
# Parámetros de la a priori para beta
beta_0=2
sigma02_beta=2ˆ2
# Parámetros de la a posteriori para beta
SSxy=sum((x-mean(x))*(y-mean(y)))
SSx=sum((x-mean(x))ˆ2)
betahat=SSxy/SSx
beta_1=beta_0*(sigma2/(sigma2+SSx*sigma02_beta)) +
betahat*(SSx*sigma02_beta/(sigma2+SSx*sigma02_beta))
sigma12_beta=sigma2*sigma02_beta/(sigma2+SSx*sigma02_beta)
par_post=c(beta_1, sigma12_beta)
round(par_post,4)

## [1] 2.3137 0.7059

5
curve(dnorm(x,beta_1,sqrt(sigma12_beta)), from=-1, to=6, xlab="beta",
ylab="Densidad Post.",lty=1, lwd=2) # posteriori
0.4
Densidad Post.

0.3
0.2
0.1
0.0

−1 0 1 2 3 4 5 6

beta

4. Un intervalo de credibilidad al 95% se obtendrá fácilmente de esta distribución de la

siguiente manera:
intervalo_beta=qnorm(c(0.025,0.975),beta_1,sqrt(sigma12_beta))
round(intervalo_beta,2)

## [1] 0.67 3.96

Es decir, con probabilidad 0.95 la pendiente de la recta de regresión está comprendida entre
0.67 y 3.96, indicando claramente su positividad.

5. El test de hipótesis

H0 : β ≤ 0 vs H1 : β > 0

lo realizamos como habitualmente, calculando las probabilidades a priori de cada hipótesis y

a continuación sus probabilidades a posteriori.
pi_0= pnorm(0,beta_0, sqrt(sigma02_beta))
pi_1=1-pi_0
cociente_priori=pi_0/pi_1
tabla_test_priori=cbind(pi_0, pi_1, cociente_priori)
colnames(tabla_test_priori)=c("Prob. priori H0", "Prob. priori H1",
"Cociente priori")
round(tabla_test_priori,3)

6
## Prob. priori H0 Prob. priori H1 Cociente priori
## [1,] 0.159 0.841 0.189

El cociente a priori nos indica que claramente, antes de incorporar los datos, preferimos H1 .
El cálculo de las probabilidades a posteriori será:
p_0= pnorm(0,beta_1,sqrt(sigma12_beta))
p_1=1-p_0
cociente_posteriori=p_0/p_1
factorBayes=cociente_posteriori/cociente_priori
tabla_test_post=cbind(p_0, p_1, cociente_posteriori,factorBayes)
colnames(tabla_test_post)=c("Prob. post. H0", "Prob. post. H1",
"Cociente post.", "Factor Bayes")
round(tabla_test_post,3)

## Prob. post. H0 Prob. post. H1 Cociente post. Factor Bayes

## [1,] 0.003 0.997 0.003 0.016

La probabilidad de ser cierta de H1 es 0.997. Elegiremos H1 , es decir, la pendiente de la recta

de regresión, β, es positiva (con probabilidad 0.997). 

Ejercicio 2.

Deseamos establecer la relación existente entre la reacción a un estímulo (variable dependiente)

y el tiempo de exposición a ese estímulo (variable predictora o explicativa). Disponemos de
14 observaciones dadas por:

Estímulo (y): 4.50 4.00 3.82 2.50 2.00 1.40 0.55 0.67 0.36 0.41 0.24 0.12 0.29 0.11
Tiempo (t): 0.5 1.00 2.00 3.00 3.50 4.00 5.00 6.00 7.00 8.00 9.00 10.00 11.00 12.00
Comparar un ajuste lineal con uno exponencial para explicar y como función de t.

Solución.
1. Realizaremos el ajuste lineal y el exponencial y compararemos los valores de r2 en cada
uno de ellos.
y=c(4.50,4.00,3.82,2.50,2.00,1.40,0.55,0.67,0.36,0.41,0.24,0.12
,0.29,0.11)
t=c(0.5,1.00,2.00,3.00,3.50,4.00,5.00,6.00,7.00,8.00,9.00,10.00,
11.00,12.00)
reg_lin=lm(y~t) # regresión lineal
reg_lin

##
## Call:
## lm(formula = y ~ t)

7
##
## Coefficients:
## (Intercept) t
## 3.6963 -0.3754
summary(reg_lin)

##
## Call:
## lm(formula = y ~ t)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.26959 -0.62732 -0.07424 0.71164 0.99133
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.69634 0.39441 9.372 7.18e-07 ***
## t -0.37535 0.05733 -6.547 2.74e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7739 on 12 degrees of freedom
## Multiple R-squared: 0.7813, Adjusted R-squared: 0.763
## F-statistic: 42.86 on 1 and 12 DF, p-value: 2.743e-05
y2=log(y) # tomamos logaritmos de y
reg_exp=lm(y2~t)
summary(reg_exp)

##
## Call:
## lm(formula = y2 ~ t)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6169 -0.1077 0.0114 0.1570 0.7403
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6834 0.1816 9.27 8.07e-07 ***
## t -0.3329 0.0264 -12.61 2.78e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3564 on 12 degrees of freedom

8
## Multiple R-squared: 0.9298, Adjusted R-squared: 0.924
## F-statistic: 159 on 1 and 12 DF, p-value: 2.779e-08
Como vemos las rectas de regresión para y y su logaritmos log(y) son:

y = 3.6963 − 0.3754 × t,

log(y) = 1.6834 − 0.3329 × t,

2 2
con coeficientes de deteminación, rlineal = 0.7813 y rexp = 0.9298, respectivamente. Elegiremos
por tanto el modelo exponencial. Gráficamente tenemos:
par(mfrow=c(1,2))
plot(t,y,pch=20, col="Blue",ylab="Estímulo",
xlab="Tiempo de exposición", [Link]=1, main="Lineal")
abline(reg_lin, col="Red")

y2=log(y)
reg_exp=lm(y2~t)
plot(t,y,pch=20, col="Blue",ylab="Estímulo",
xlab="Tiempo de exposición", [Link]=1, main="Exponencial")
lines(t,exp(reg_exp$coefficients[1]+reg_exp$coefficients[2]*t),type="l",
col="Red")
Lineal Exponencial
4

4
3

3
Estímulo

Estímulo
2

2
1

1
0

2 4 6 8 10 2 4 6 8 10

Tiempo de exposición Tiempo de exposición 

9
Ejercicio 3.

Un investigador desea conocer la posible relación lineal entre el consumo de carburante y la

velocidad de conducción. Se han tomado seis medidas para diferentes velocidades y distancias
recorridas con un litro de carburante. La velocidad (en kilómetros por hora) y distancias (en
kilómetros) recorridas vienen dadas en la siguiente tabla:
Velocidad Distancia
x y
80 55.7
90 55.4
100 52.5
110 52.1
120 50.5
130 49.2
Se pide:
1. Realizar el ajuste lineal, obtener las cantidades de interés para dicho ajuste y representar
gráficamente.
2. Supongamos que y dado x sigue una distribución normal de media α + β × x y varianza
σ 2 = 0.572 conocida. Usar una previa para β, normal con media (a priori), 0 y varianza
(a proiri), 1 y obtener la a posteriori de β.
3. Realizar el test

H0 : β ≥ 0 vs H1 : β < 0

4. Realizar el test

H0 : β = 0 vs H1 : β 6= 0

al 95%.
5. Predecir puntualmente la distancia recorrida y para una futura observación que haya
alcanzado una velocidad de x = 115 kms por hora. Dar un intervalo de credibilidad para
la predicción al 95%. Nota: usar una a priori para αx̄ normal de media 100 varianza 1.

Solución.

1. Procedemos de forma análoga al Ejercicio 1.

x=c(80,90,100,110,120,130)
y=c(55.7,55.4, 52.5, 52.1, 50.5, 49.2)
n=length(x)
round(mean(x),4)

## [1] 105

10
round(mean(y),4)

## [1] 52.5667
SSxy=sum((x-mean(x))*(y-mean(y)))
SSx=sum((x-mean(x))ˆ2)
SSy=sum((y-mean(y))ˆ2)
betahat=SSxy/SSx
alfahat=mean(y)-mean(x)*betahat
coeficientes=c(alfahat,betahat)
round(coeficientes,4)

## [1] 66.8467 -0.1360

sigma2=sqrt(sum((y-(alfahat+betahat*x))ˆ2)/(n-2))
round(sigma2,4)

## [1] 0.5713
plot(x,y,pch=20, col="Blue",ylab="Distancia",xlab="Velocidad")
abline(alfahat,betahat, col="Red")
55
54
Distancia

53
52
51
50
49

80 90 100 110 120 130

Velocidad
r2=SSxyˆ2/(SSx*SSy)
round(r2,4)

## [1] 0.9612

La recta de regresión ajustada es y = 66.847 − 0.136 × x, o equivalentemente su forma

alternativa y = 52.5667 − 0.136(x − 105). El coeficiente de determinación es alto, r2 = 0.9612,

11
 √
indicando que la relación lineal parece adecuada. El coeficiente de correlación es r = r2 =
0.98, indicando una dependencia directa entre y y x. La pendiente de la recta nos indica
que por cada incremento unitario en la velocidad y, la distancia recorrida con un litro de
combustible, disminuye en 0.316 kilómetros.

Utilizando el comando lm de R obtendremos los mismos resultados.

x=c(80,90,100,110,120,130)
y=c(55.7,55.4, 52.5, 52.1, 50.5, 49.2)
recta=lm(y~x)
summary(recta)

##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## 1 2 3 4 5 6
## -0.26667 0.79333 -0.74667 0.21333 -0.02667 0.03333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.84667 1.45268 46.016 1.33e-06 ***
## x -0.13600 0.01366 -9.959 0.000571 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5713 on 4 degrees of freedom
## Multiple R-squared: 0.9612, Adjusted R-squared: 0.9515
## F-statistic: 99.19 on 1 and 4 DF, p-value: 0.000571
# Coeficientes de regresión
recta$coefficients

## (Intercept) x
## 66.84667 -0.13600
# Residual standard error
res=summary(recta)
sigma=res$sigma
round(sigma,4)

## [1] 0.5713
#Varianza estimada
sigma2=res$sigmaˆ2
round(sigma2,4)

12
## [1] 0.3263
# Coeficiente de determinación R2
r2=res$[Link]
r2

## [1] 0.9612354
# Coeficiente de correlación
r=sqrt(r2)
r

## [1] 0.9804261
# Gráfica
plot(x,y,pch=20, col="blue",ylab="Distancia",xlab="Velocidad")
abline(recta, col="red")
55
54
Distancia

53
52
51
50
49

80 90 100 110 120 130

Velocidad
2. Las cantidades a posteriori y la densidad a posteriori se obtienen con el código siguiente:
sigma2=0.57ˆ2
beta_0=0
sigma02_beta=1ˆ2
beta_1=beta_0*(sigma2/(sigma2+SSx*sigma02_beta))+
betahat*(SSx*sigma02_beta/(sigma2+SSx*sigma02_beta))
sigma12_beta=sigma2*sigma02_beta/(sigma2+SSx*sigma02_beta)
par_post=c(beta_1, sigma12_beta)
round(par_post,4)

## [1] -0.1360 0.0002

13
curve(dnorm(x,beta_1,sqrt(sigma12_beta)), from=-0.5, to=0.5, xlab="beta",
ylab="Densidad Post.",lty=1, lwd=2) # posteriori
25
20
Densidad Post.

15
10
5
0

−0.4 −0.2 0.0 0.2 0.4

beta

Como vemos la distribución a posteriori de β está muy concentrada en torno a −0.1.

3. Para el test de hipótesis tenemos:

pi_0= 1-pnorm(0,beta_0, sqrt(sigma02_beta))
pi_1=1-pi_0
cociente_priori=pi_0/pi_1
tabla_test_priori=cbind(pi_0, pi_1, cociente_priori)
colnames(tabla_test_priori)=c("Prob. priori H0", "Prob. priori H1",
"Cociente priori")
round(tabla_test_priori,3)

## Prob. priori H0 Prob. priori H1 Cociente priori

## [1,] 0.5 0.5 1
A priori no tenemos una hipótesis preferida.
p_0= 1-pnorm(0,beta_1,sqrt(sigma12_beta))
p_1=1-p_0
cociente_posteriori=p_0/p_1
factorBayes=cociente_posteriori/cociente_priori
tabla_test_post=cbind(p_0, p_1, cociente_posteriori,factorBayes)
colnames(tabla_test_post)=c("Prob. post. H0", "Prob. post. H1",
"Cociente post.", "Factor Bayes")
round(tabla_test_post,3)

14
## Prob. post. H0 Prob. post. H1 Cociente post. Factor Bayes
## [1,] 0 1 0 0

Es obvio, que aceptaremos H1 (es decir, β < 0).

4. Para el test

H0 : β = 0 vs H1 : β 6= 0,

calcularemos el intervalo de credibilidad de β al 95%.

intervalo_beta=qnorm(c(0.025,0.975),beta_1,sqrt(sigma12_beta))
round(intervalo_beta,2)

## [1] -0.16 -0.11

Dado que dicho intervalo no contiene al 0, no tenemos evidencia para aceptar que β = 0 al
95%.

5. Finalmente si queremos predecir un valor futuro de yn+1 dado xn+1 , sabemos que la
distribución predictiva es normal de media,

µn+1 = α1 + β1 · (xn+1 − x̄),

y varianza,
2 2
σn+1 = σ1,α + (xn+1 − x̄)2 σ1,β
2
+ σ2.

alfa_0=100
sigma02_alfa=1ˆ2
alfa_1=alfa_0*(sigma2/(sigma2+n*sigma02_alfa))+
mean(y)*(n*sigma02_alfa/(sigma2+n*sigma02_alfa))
sigma12_alfa=sigma2*sigma02_alfa/(sigma2+n*sigma02_alfa)
x_fut=115
mu_pred=alfa_1+beta_1*(x_fut-mean(x))
sigma2_pred=sigma12_alfa+sigma12_beta*(x_fut-mean(x))ˆ2+sigma2
round(mu_pred,2)

## [1] 53.64
intervalo_pred=qnorm(c(0.025,0.975),mu_pred,sqrt(sigma2_pred))
round(intervalo_pred,2)

## [1] 52.41 54.88

Puntualmente, la predicción para la distancia será de 53.64 kilómetros. Con probabilidad

0.95, la próxima observación que haya llevado una velocidad de 115 kms. por hora, recorrerá
entre 52.41 y 54.88 kilómetros. 

15
Ejericio 4.

Se desea establecer la relación que pueda existir entre el coste del transporte de un material
altamente inflamable y su tamaño. Estamos interesados en construir un modelo lineal que
prediga el coste del transporte en función del tamaño de dicho material. Disponemos de datos
para los últimos 10 viajes realizados.
Coste transporte 11.8 16.1 14.0 20.5 20.7 16.1 18.0 7.9 9.5 13.7
Tamaños (cms. cuadrados) 2.9 5.0 3.8 6.5 8.0 4.8 5.6 1.1 2.4 4.2
Se trata de dar un modelo lineal para predecir el precio del transporte en función del tamaño
del material y representar gráficamente la relación encontrada.

Solución. Buscaremos un ajuste lineal para estos datos.

y=c(11.8,16.1,14.0, 20.5, 20.7, 16.1, 18.0, 7.9, 9.5, 13.7)
x=c(2.9,5.0,3.8, 6.5,8.0, 4.8, 5.6, 1.1, 2.4, 4.2)
n=length(x)
xbar=mean(x)
ybar=mean(y)
medias=c(xbar,ybar)
round(medias,4)

## [1] 4.43 14.83

reg_lin=lm(y~x)
reg_lin

##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept) x
## 5.587 2.087
summary(reg_lin)

##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5789 -0.4830 0.1215 0.4946 1.3509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)

16
## (Intercept) 5.5867 0.7431 7.518 6.81e-05 ***
## x 2.0865 0.1539 13.557 8.42e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9344 on 8 degrees of freedom
## Multiple R-squared: 0.9583, Adjusted R-squared: 0.9531
## F-statistic: 183.8 on 1 and 8 DF, p-value: 8.419e-07
plot(x,y,pch=20, col="Blue",ylab="Coste",xlab="Tamaño", ylim=c(5,22))
abline(reg_lin, col="Red")
20
15
Coste

10
5

1 2 3 4 5 6 7 8

Tamaño
La recta de regresión obtenida es

yi = 5.5867 + 2.0865xi ,

El término α = 5.5867, nos indica que, en media, transportar cualquier cantidad de material
inflamable costará como mínimo 5.59 euros. La pendiente estimada nos indica que un aumento
en 1 cm. cuadrado del material supondrá, en media, un incremento del precio de 2.1 euros.
El coeficiente de determinación r2 = 0.9583 próximo a 1 indica un buen ajuste.

El modelo en su forma reparametrizada alternativa es:

yi = α̃x̄ + β̃(x − x̄) = 14.83 + 2.0865 · (x − 4.43)

La salida proporcionada por summary contiene la siguiente información:

17
 1. call presenta la expresión de la regresión que hemos introducido. En nuestro caso, y
como función lineal de x.
2. Residuals contiene valores descriptivos de los residuos que nos pueden resultar de
utilidad para inspeccionar la hipótesis de normalidad de los errores. En este caso,
tenemos que la mediana está proxima a 0 y los cuartiles primero y tercero tienen
prácticamente el mismo valor pero de signo contrario. El mismo razonamiento vale para
para los valores extremos, máximo y mínimo, indicando en general un comportamiento
simétrico en torno a la mediana.
3. El apartado Coefficients contiene la estimación puntual de los coeficientes de la
recta de regresión. En la columna Estimate tenemos la de α en la fila Intercept y la
de β en la de x que como sabemos son las estimaciones dadas por los estimadores de
máxima verosimilitud de los parámetros cuya desviación estándar (Std. Error) nos
informa sobre la variabilidad de éstas. En nuestro caso valen 0.7431 y 0.1539. Con estas
cantidades podemos realizar los test clásicos de la regresión lineal. Básicamente son los
referidos a la nulidad de los coeficientes, H0 : α = 0 y H0 : β = 0, utilizando para ello
los valores del estadístico t de Student proporcionados en la columna t value (7.518
y 13.557, en nuestro caso) y sus correspondientes p−valores en la columna Pr(>|t|)
(prácticamente 0 en este ejemplo, rechazando por tanto en ambos casos la hipótesis de
que vale 0). Suele decirse entonces que la variable es significativa pues tiene influencia
sobre la predicción de la variable respuesta.
4. Residual standard error, proporciona una estimación de la desviación típica de los
residuos.
5. Multiple R-squared corresponde al coeficiente de determinación (r2 ). Además de
nuestra interpretación sobre la bondad del ajuste lineal (0.9583 en este ejercicio indicando
un buen ajuste), suele interpretarse como la proporción de la variable dependiente (y)
explicada por el modelo. En nuestro caso, el modelo recoge el 95.83% de la variabilidad
de los datos. Este estadístico tiene el inconveniente de que si añadimos nuevas variables
al modelo, independientemente de su poder predictivo, su valor aumenta. Por esta
razón, el valor de Adjusted R-squared tiene en cuenta el número de variables del
modelo (k = 1 en este ejercicio) y el tamaño de la muestra n (10 en este caso) y tiene
una interpretación análoga a la de r2 (en este ejercicio su valor es también alto, 0.9531).
6. F-statistic es el valor del estadístico para el contraste de que todos los coeficientes
conjuntamente valen 0 frente a la alternativa de que alguno de ellos no es nulo. En
nuestro ejercicio, dado que el p−valor es prácticamente 0 rechazaremos la hipótesis
nula y por tanto el modelo es globalmente significativo.

2. Para el caso no informativo tenemos que las a posteriori de cada uno de los parámetros
tienen distribuciones normales:
!
σ2
π(β|y, x) = N β | β,
b .
SSx

18
 !
σ2
π(αx̄ |y, x) = N αx̄ | α
b x̄ , .
n
n 
SSxy 1 X 2
donde βb = b x̄ = ȳ, y σ 2 viene dado por
,α yi − (α̃x̄ + β̃(xi − x̄)) .
SSx n − 2 i=1
SSx=sum((x-xbar)ˆ2)
SSxy=sum((x-xbar)*(y-ybar))
betahat=SSxy/SSx
sigma2=sum((y-(ybar+betahat*(x-xbar)))ˆ2)/(n-2)
beta_1=betahat
sigma12_beta=sigma2/SSx
alfa_1=mean(y)
sigma12_alfa=sigma2/n
# Intervalos de credibilidad
intervalo_alfa=qnorm(c(0.025,0.975),alfa_1,sqrt(sigma12_alfa))
round(intervalo_alfa,3)

## [1] 14.251 15.409

intervalo_beta=qnorm(c(0.025,0.975),beta_1,sqrt(sigma12_beta))
round(intervalo_beta,3)

## [1] 1.785 2.388

Es interesante observar que con la instrucción residuals podemos obtener el valor del error
cuadrático σc2 . Análogamente con deviance.
sum(residuals(reg_lin)ˆ2)/(n-2)

## [1] 0.8731383
sigma2

## [1] 0.8731383
deviance(reg_lin)/(n-2)

## [1] 0.8731383
Si nos situamos en el caso más realista de varianza desconocida, sabemos que entonces
debemos usar la distribución t de Student.
ext_inf_alfa=alfa_1+qt(0.025,n-2)*sqrt(sigma12_alfa)
ext_sup_alfa=alfa_1-qt(0.025,n-2)*sqrt(sigma12_alfa)
int_alfa=c(ext_inf_alfa, ext_sup_alfa)
round(int_alfa,4)

## [1] 14.1486 15.5114

19
ext_inf_beta=beta_1+qt(0.025,n-2)*sqrt(sigma12_beta)
ext_sup_beta=beta_1-qt(0.025,n-2)*sqrt(sigma12_beta)
int_beta=c(ext_inf_beta,ext_sup_beta)
round(int_beta,4)

## [1] 1.7316 2.4414

Dado que sabemos que el caso no informativo coincide con el caso de máxima verosimilitud
podemos utilizar directamente las salidas de R. Por ejemplo, el comando confint nos
proporciona los intervalos de confianza de los parámetros que coincidirán con los de credibilidad
(bayesianos).
confint(reg_lin, level=0.95)

## 2.5 % 97.5 %
## (Intercept) 3.873185 7.300299
## x 1.731604 2.441424
Que como vemos coincide con el de credibilidad de probabilidad 0.95 obtenido anteriormente.
En el caso del intercept se presenta el intervalo para α y nosotros hemos obtenido anterior-
mente el intervalo para el modelo centrado en x̄, es decir, αx̄ . Para ver que coinciden, basta
hacer la regresión para la recta centrada en x̄.
x1=x-xbar
recta=lm(y~x1)
confint(recta,level=0.95)

## 2.5 % 97.5 %
## (Intercept) 14.148601 15.511399
## x1 1.731604 2.441424

3. Para los test de hipótesis tenemos lo siguiente. En el caso, del test

H0 : β = 0 vs H1 : β 6= 0,
ó
H0 : α = 0 vs H1 : α 6= 0,
tenemos que los intervalos de credibilidad al 95% para ambos parámetros eran (3.873, 7.300)
para α y (1.731, 2.441) para β y dado que ninguno de los dos contiene al 0 se rechazan las
hipótesis de que son 0. Tal y como ocurría con los p−valores en summary(recta) que daban
los dos valores muy pequeños (altamente significativos) y por tanto se rechazaba la hipótesis
de nulidad de cada parámetro.
Finalmente, el test de hipótesis
H0 : β > 0 vs H1 : β ≤ 0,
lo realizamos como habitualmente con el cálculo de las probabilidades (a posteriori) de cada
hipótesis:

20
p_0=1-pnorm(0,beta_1,sqrt(sigma12_beta))
p_1=1-p_0
cociente_post=p_0/p_1
test_post=cbind(p_0,p_1,cociente_post)
colnames(test_post)=c("Prob. post. H0","Prob. post. H1", "Cociente post.")
round(test_post,3)

## Prob. post. H0 Prob. post. H1 Cociente post.

## [1,] 1 0 Inf
Y por tanto, aceptaremos H0 .

4. La predicción del coste para un nuevo transporte de material de tamaño 6.3 centímetros
2
cuadrados tendrá una distribución Normal de media µn+1 y varianza σn+1 con:
x_fut=6.3
mu_pred=alfa_1+beta_1*(x_fut-mean(x))
sigma2_pred=sigma12_alfa+sigma12_beta*(x_fut-mean(x))ˆ2 + sigma2
round(mu_pred,3)

## [1] 18.732
# Intervalo para la predicción. Usamos el caso varianza desconocida y por tanto
# la distribución t de Student
lb_pred=mu_pred+qt(0.025,n-2)*sqrt(sigma2_pred)
ub_pred=mu_pred-qt(0.025,n-2)*sqrt(sigma2_pred)
int_pred=c(lb_pred,ub_pred)
round(int_pred,3)

## [1] 16.376 21.087

Podríamos obtenerlo directamente desde R con la instrucción predict (como sabemos ambos
procedimientos coinciden cuando realizamos un análisis basado en datos, sin incorporación
de información a priori).
prediccion=[Link](x=c(6.3))
ajuste=predict(reg_lin,newdata=prediccion, interval="prediction")
round(ajuste,3)

## fit lwr upr

## 1 18.732 16.376 21.087


Ejercicio 5 (caso de variable explicativa dicotómica).

Los datos smoke, disponibles en el paquete wooldridge, son referidos al consumo diario de
tabaco (variable cigs) de 807 individuos y una serie de características socioeconómicas así

21
como el lugar de residencia. Una de esta variables corresponde a la variable dicotómica (0 ó
1) que informa sobre si el estado en el que reside el individuo tiene o no restricciones legales
para fumar en los restaurantes (variable restaurn). Construir un modelo de regresión lineal
para analizar el efecto de imponer dichas restricciones legales en los restaurantes.

Solución. Debemos cargar el paquete wooldridge con la instrucción [Link](wooldridge)

y a continuación los datos smoke.
library(wooldridge)
data(smoke)

El modelo de regresión lineal expresará la variable respuesta (cigs) como función lineal de
la variable explicativa, restricción para fumar en los restaurantes (restaurn), que es una
variable dicotómica que toma el valor 1 si el estado no permite fumar en los restaurantes, y 0
en caso contrario.
reg_lin=lm(cigs~restaurn, data=smoke)
summary(reg_lin)

##
## Call:
## lm(formula = cigs ~ restaurn, data = smoke)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.370 -9.370 -6.598 10.630 70.630
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.3701 0.5547 16.892 <2e-16 ***
## restaurn -2.7721 1.1171 -2.482 0.0133 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.68 on 805 degrees of freedom
## Multiple R-squared: 0.007592, Adjusted R-squared: 0.006359
## F-statistic: 6.158 on 1 and 805 DF, p-value: 0.01328
Lo interesante de este ejercicio no es el modelo de regresión obtenido, que como vemos es
pobre, dado que la variable explicativa sólo toma dos posibles valores. Lo importante ahora
es la interpretación de los coeficientes. El término constante α̂ = 9.37, nos indica que la
media de cigarrillos diarios consumidos por los individuos de estados donde se permite fumar
es de 9.37 y la pendiente β̂ = −2.77, nos indica la diferencia del consumo diario medio de
cigarrillos entre los estados con restricción y los que no la tienen. En efecto,
mean(smoke$cigs[smoke$restaurn==0])

## [1] 9.370066

22
mean(smoke$cigs[smoke$restaurn==1])-mean(smoke$cigs[smoke$restaurn==0])

## [1] -2.772076
Para ver si estas diferencias son signficativas estadísticamente lo podemos hacer de dos formas:
a) Dado el carácter dicotómico de x, las diferencias entre grupos (tener o no restricciones)
equivale ver si H0 : β = 0 ó H1 : β 6= 0. El p−valor de 0.0133 nos indica que al 95%
rechazamos la hipótesis nula de que β = 0 y por tanto el consumo de cigarrillos es diferente
entre estados con o sin restricción. Análogamente, si calculamos los intervalos de credibilidad
para β (o confianza, ahora coinciden como sabemos) tenemos:
confint(reg_lin, level=0.95)

## 2.5 % 97.5 %
## (Intercept) 8.281219 10.458912
## restaurn -4.964764 -0.579388
Es decir, (−4.96, −0.58) que no contiene al 0 y por tanto se rechaza H0 al 95%.
b) Otra forma equivalente, consiste en realizar un test de hipótesis [Link] como el realizado
en el Tema 1 para diferencias de medias (caso no informativo) con varianzas iguales.
[Link](cigs~restaurn, data=smoke, [Link]=TRUE)

##
## Two Sample t-test
##
## data: cigs by restaurn
## t = 2.4816, df = 805, p-value = 0.01328
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.579388 4.964764
## sample estimates:
## mean in group 0 mean in group 1
## 9.370066 6.597990
Vemos que el p−valor y el valor de estadístico t son iguales (cambiado de signo en t pero esto
se debe a que contrasta el grupo de mayor valor con el de menos, es decir, sin restricciones
frente a restricciones). La decisión es la misma, se rechaza la igualdad de medias. 

Ejercicio 6.

Este ejercicio trata sobre la importancia de representar nuestro conjunto de datos antes de
realizar una regresión lineal y la precaución que debe tomarse. Los datos anscombe del paquete
MASS, contienen pares de variables (Xi , Yi ), i = 1, . . . , 4 para construir las correspondientes
rectas de regresión. Realizarlas y representarlas gráficamente.

23
Solución. Cargamos el paquete MASS con [Link](MASS) y posteriormente los
datos anscombe.
library(MASS)

##
## Attaching package: 'MASS'
## The following object is masked from 'package:wooldridge':
##
## cement
data(anscombe)

Ahora hacemos las regresiones lineales.

reg_1=lm(y1~x1, data=anscombe)
reg_2=lm(y2~x2, data=anscombe)
reg_3=lm(y3~x3, data=anscombe)
reg_4=lm(y4~x4, data=anscombe)
summary(reg_1)

##
## Call:
## lm(formula = y1 ~ x1, data = anscombe)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.92127 -0.45577 -0.04136 0.70941 1.83882
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0001 1.1247 2.667 0.02573 *
## x1 0.5001 0.1179 4.241 0.00217 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared: 0.6665, Adjusted R-squared: 0.6295
## F-statistic: 17.99 on 1 and 9 DF, p-value: 0.00217
summary(reg_2)

##
## Call:
## lm(formula = y2 ~ x2, data = anscombe)
##
## Residuals:

24
## Min 1Q Median 3Q Max
## -1.9009 -0.7609 0.1291 0.9491 1.2691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.001 1.125 2.667 0.02576 *
## x2 0.500 0.118 4.239 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared: 0.6662, Adjusted R-squared: 0.6292
## F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002179
summary(reg_3)

##
## Call:
## lm(formula = y3 ~ x3, data = anscombe)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1586 -0.6146 -0.2303 0.1540 3.2411
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0025 1.1245 2.670 0.02562 *
## x3 0.4997 0.1179 4.239 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared: 0.6663, Adjusted R-squared: 0.6292
## F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002176
summary(reg_4)

##
## Call:
## lm(formula = y4 ~ x4, data = anscombe)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.751 -0.831 0.000 0.809 1.839
##
## Coefficients:

25
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0017 1.1239 2.671 0.02559 *
## x4 0.4999 0.1178 4.243 0.00216 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared: 0.6667, Adjusted R-squared: 0.6297
## F-statistic: 18 on 1 and 9 DF, p-value: 0.002165
Los cuatro modelos tienen coeficientes y p−valores muy similares con sus respectivas variables
explicativas significativas. La bondad de los ajustes también es muy similar, con valores de
los coeficientes de determinación alrededor de 0.66. En general, parece que todas las variables
independientes tienen el mismo poder explicativo y que todos los modelos presentan el mismo
nivel de precisión. Veamos sus gráficas:
par(mfrow=c(2,2))
plot(anscombe$x1,anscombe$y1, pch=20, col="Blue", xlab="X_1",ylab="Y_1")
abline(reg_1, col="Red")
plot(anscombe$x2,anscombe$y2, pch=20, col="Blue", xlab="X_2",ylab="Y_2")
abline(reg_2, col="Red")
plot(anscombe$x3,anscombe$y3, pch=20, col="Blue", xlab="X_3",ylab="Y_3")
abline(reg_3, col="Red")
plot(anscombe$x4,anscombe$y4, pch=20, col="Blue", xlab="X_4",ylab="Y_4")
abline(reg_4, col="Red")
11

9
7
Y_1

Y_2
4 6 8

5
3

4 6 8 10 12 14 4 6 8 10 12 14

X_1 X_2
10
10
Y_3

Y_4

6
6

4 6 8 10 12 14 8 10 12 14 16 18

X_3 X_4

Como vemos, sólo las relaciones (X1 , Y1 ) y (X3 , Y3 ) pueden ajustarse linealmente, ya que el

26
resto no cumple el supuesto de linealidad. La relación entre Y3 y X3 parece afectada por una
observación atípica. Para estos caso, se recomienda el uso de estimadores robustos que no se
vean tan afectados por este tipo de observaciones. El comando rlm realiza este ajuste.
reg_rob=rlm(y3~x3, data=anscombe)
plot(y3~x3, data=anscombe, pch=20, col="Blue", xlab="X_3", ylab="Y_3")
abline(reg_3,lty=2)
abline(reg_rob,col="Red")
12
10
Y_3

8
6

4 6 8 10 12 14

X_3 

Regresión lineal múltiple

Ejercicio 7.

Utilizando el conjunto de datos hprice1 en el paquete wooldridge vamos a realizar una

regresión lineal múltiple que nos permita explicar el precio (en logaritmo) de las viviendas
(variable price) en función del número de habitaciones (bdrms), las dimensiones de la parcela
donde se ubica (lotsize) en logaritmo, el tamaño de la casa en pies cuadrados (sqrft) en
logaritmo y el estilo de la casa (colonial). Esta última variable toma el valor 1 si es colonial
y 0 en caso contrario.

Solución. Tenemos un modelo de regresión lineal múltiple con 4 variables explicativas. La

resolución con R es análoga al caso univariante.
library(wooldridge)
data(hprice1)

27
reg_multiple=lm(log(price)~bdrms+log(lotsize)+log(sqrft)+colonial, data=hprice1)
summary(reg_multiple)

##
## Call:
## lm(formula = log(price) ~ bdrms + log(lotsize) + log(sqrft) +
## colonial, data = hprice1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69479 -0.09750 -0.01619 0.09151 0.70228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.34959 0.65104 -2.073 0.0413 *
## bdrms 0.02683 0.02872 0.934 0.3530
## log(lotsize) 0.16782 0.03818 4.395 3.25e-05 ***
## log(sqrft) 0.70719 0.09280 7.620 3.69e-11 ***
## colonial 0.05380 0.04477 1.202 0.2330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1841 on 83 degrees of freedom
## Multiple R-squared: 0.6491, Adjusted R-squared: 0.6322
## F-statistic: 38.38 on 4 and 83 DF, p-value: < 2.2e-16
Como vemos, tan sólo el término constante junto con las variables sqrft (tamaño de la
casa) y lotsize (dimensiones de la parcela) influyen estadísticamente en el precio (son
significativas), es decir, se rechaza la hipótesis nula de que valen 0. El valor del coeficiente
de determinación es 0.6491 indicando un relativamente buen ajuste lineal y que el modelo
es capaz de explicar un 64.91% de la variabilidad del precio de las viviendas (en logaritmo).
El p−valor en el estadístico F también nos indica que rechazamos que conjuntamente las
variables sean 0 y por tanto, el modelo tiene cierta capacidad predictora conjunta.
Es decir, nuestro modelo lineal estimado responde a la ecuación:

log(pricei ) = −1.35+0.027×bdrmsi +0.17×log(lotsizei )+0.71×log(sqrfti )+0.05×coloniali +εi

La interpretación de los coeficientes es algo diferente al caso univariante. Fijémonos por

ejemplo en la variable sqrft. En este caso, diremos que si mantenemos las demás variables
constantes, el precio de la vivienda incrementa 0.71% dólares si se incrementa el tamaño de la
casa un 1%. La interpretación por tanto debe incorporar el hecho de que las demás variables
deben mantenerse constantes.
Una práctica habitual en muchas aplicaciones consiste en lo siguiente:
1. Comenzar nuestro modelo un número “grande” de variables X1 , . . . , Xk .

28
 2. Obtener la regresión lineal de Y sobre ellas.
3. Eliminar toda variable Xi para la cual el estimador β̂i tenga un p−valor asociado alto
(mayor de 0.1).
4. Volver a realizar la regresión sobre las variables restantes.
reg_mult2=lm(log(price)~log(lotsize)+log(sqrft), data=hprice1)
summary(reg_mult2)

##
## Call:
## lm(formula = log(price) ~ log(lotsize) + log(sqrft), data = hprice1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.65332 -0.11051 -0.00648 0.11822 0.66417
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.64007 0.60188 -2.725 0.00781 **
## log(lotsize) 0.16846 0.03846 4.380 3.37e-05 ***
## log(sqrft) 0.76237 0.08089 9.425 7.44e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1855 on 85 degrees of freedom
## Multiple R-squared: 0.6353, Adjusted R-squared: 0.6267
## F-statistic: 74.04 on 2 and 85 DF, p-value: < 2.2e-16
Presenta por tanto una estimación muy parecida a la anterior, pero en la que hemos reducido
el número de variables explicativas.
Uno de los problemas habituales en la regresión múltiple, se conoce como la multicolineali-
dad. Básicamente, ésta se presenta cuando dos variables explicativas tienen correlaciones
alta y ésto hace que sus coeficientes asociados sean inestables. También, puede presentarse
cuando una de las variables está altamente correlacionada con una combinación lineal del
resto. Para examinar la presencia de multicolinealidad suele usarse el factor de inflación
de la varianza (variance inflation factor, vif). Este indicador, nos informa sobre el grado
en el que la varianza de una coeficiente se incrementa por la presencia de multicolinealidad.
Su expresión es:
1
vif(Xi ) = 2
,
1 − rX −i

2
donde rX −i
es el coeficiente de determinación del modelo lineal que predice Xi en función
del resto de variables explicativas. Su valor es siempre mayor o igual que 1 (cuando el resto
de variables están completamente incorreladas con la variable analizada) y tiende a infinito
cuando la correlación tienda a 1. Valores del vif para una variable Xi mayor de 10 aconsejan

29
eliminar esa variable.
library(car)

## Loading required package: carData

vif(reg_multiple)

## bdrms log(lotsize) log(sqrft) colonial

## 1.499012 1.107418 1.479106 1.106776
vif(reg_mult2)

## log(lotsize) log(sqrft)
## 1.107306 1.107306
Los modelos no presentan ninguna variable con multicolinealidad elevada.
Supongamos ahora que deseamos predecir el valor de una casa. Por ejemplo, la predicción
en el mercado de las siguiente dos casas:
Casa 1: 3 habitaciones, 1845 pies cuadrados, colonial y tamaño de la parcela 8328 pies
cuadrados.
Casa 2: 2 habitaciones, 1972 pies cuadrados en parcela de 7040 pies cuadrados de tipo no
colonial.

Procedemos de la siguiente forma:

prediccion=[Link](bdrms=c(3,2),lotsize=c(1845,1972), sqrft=c(8328,7040),
colonial=c(1,0))
ajuste=predict(reg_multiple, newdata=prediccion, interval="prediction")
ajuste

## fit lwr upr

## 1 6.430836 5.935302 6.926371
## 2 6.242562 5.750379 6.734745
exp(ajuste)

## fit lwr upr

## 1 620.6929 378.1541 1018.7901
## 2 514.1742 314.3098 841.1289
Vemos que el precio estimado para la Casa tipo 1, es de 620692.9 dólares, con un intervalo al
95% de (378154.1 − 1018790.1). Análogamente para la casa tipo 2. 

30
Regresión logística

Ejercicio 8.

Los datos del fichero [Link] contienen información sobre las solicitudes de admisión
de candidatos al grado en la Universidad de California en Los Ángeles (UCLA). La variable
binaria admit toma el valor 0 cuando el alumno no es admitido y 1 cuando si lo es. Vamos a
modelizar la relación entre la probabilidad de admisión y las demás variables: gre (graduate
requirements), gpa (grade point average) y rank. Y predecir la probabilidades de admisión de
un candidato con gre=790, gpa=3.8 y rank=1.

Solución. Vamos en primer lugar a leer los datos y realizar la estimación con la instrucción
glm. Debemos convertir la variable rank en un factor.
binary=[Link]("~/Desktop/[Link]")
df=binary
str(df) # vemos los valores de las variables

## '[Link]': 400 obs. of 4 variables:

## $ admit: int 0 1 1 1 0 1 1 0 1 0 ...
## $ gre : int 380 660 800 640 520 760 560 400 540 700 ...
## $ gpa : num 3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ...
## $ rank : int 3 3 1 4 4 2 1 2 3 2 ...
summary(df) # vemos una descriptiva

## admit gre gpa rank

## Min. :0.0000 Min. :220.0 Min. :2.260 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:520.0 1st Qu.:3.130 1st Qu.:2.000
## Median :0.0000 Median :580.0 Median :3.395 Median :2.000
## Mean :0.3175 Mean :587.7 Mean :3.390 Mean :2.485
## 3rd Qu.:1.0000 3rd Qu.:660.0 3rd Qu.:3.670 3rd Qu.:3.000
## Max. :1.0000 Max. :800.0 Max. :4.000 Max. :4.000
xtabs(~admit+rank, data=df)
# vemos la distribución de admitidos por rank

## rank
## admit 1 2 3 4
## 0 28 97 93 55
## 1 33 54 28 12
# Realizamos la estimación logit
logit=glm(admit~gre+gpa+[Link](rank), data=df, family=binomial)
summary(logit)

##

31
## Call:
## glm(formula = admit ~ gre + gpa + [Link](rank), family = binomial,
## data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6268 -0.8662 -0.6388 1.1490 2.0790
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.989979 1.139951 -3.500 0.000465 ***
## gre 0.002264 0.001094 2.070 0.038465 *
## gpa 0.804038 0.331819 2.423 0.015388 *
## [Link](rank)2 -0.675443 0.316490 -2.134 0.032829 *
## [Link](rank)3 -1.340204 0.345306 -3.881 0.000104 ***
## [Link](rank)4 -1.551464 0.417832 -3.713 0.000205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.52
##
## Number of Fisher Scoring iterations: 4
library(mfx)
logitmfx(admit~gre+gpa+[Link](rank), data=df)

## Call:
## logitmfx(formula = admit ~ gre + gpa + [Link](rank), data = df)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## gre 0.00047428 0.00022816 2.0787 0.03764 *
## gpa 0.16840573 0.06914174 2.4357 0.01486 *
## [Link](rank)2 -0.13608352 0.06128009 -2.2207 0.02637 *
## [Link](rank)3 -0.24649509 0.05466303 -4.5094 6.502e-06 ***
## [Link](rank)4 -0.25096705 0.04914826 -5.1063 3.285e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##

32
## [1] "[Link](rank)2" "[Link](rank)3" "[Link](rank)4"
Para cuantificar el efecto de cada variable sobre la probabilidad de admisión utilizamos la
instrucción logitmfx (cargando previamente el paquete mfx). Por ejemplo, para el caso de la
variable gpa, nos indica que evaluando el resto de las variables en su media, la probabilidad
de ser admitido aumenta en 0.17 si el individuo aumenta un punto su gpa.

Para la predicción procedemos:

x_new=[Link](gre=790, gpa=3.8, rank=[Link](1))
predict(logit,x_new)

## 1
## 0.85426
p
El valor z = 0.85426 corresponde al valor de log( ) = z, luego para conocer la probabilidad
1−p
p debemos despejar p en esta expresión. Es decir:

p p exp(z)
log( ) = z ⇐⇒ = exp(z) ⇐⇒ p = .
1−p 1−p 1 + exp(z)

z=predict(logit,x_new)
exp(z)/(1+exp(z))

## 1
## 0.70146
Como vemos, estimamos la probabilidad de admisión para este candidato alrededor del 70%.
El valor de la probabilidad puede obtenerse directamente introduciendo el comando
type="response" en predict.
predict(logit,x_new,type="response")

## 1
## 0.70146

El mismo análisis con el modelo probit se hace de forma análoga.

probit=glm(admit~gre+gpa+[Link](rank), data=df,
family=binomial(probit))
summary(probit)

##
## Call:
## glm(formula = admit ~ gre + gpa + [Link](rank),
## family = binomial(probit), data = df)
##
## Deviance Residuals:

33
## Min 1Q Median 3Q Max
## -1.6163 -0.8710 -0.6389 1.1560 2.1035
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.386836 0.673946 -3.542 0.000398 ***
## gre 0.001376 0.000650 2.116 0.034329 *
## gpa 0.477730 0.197197 2.423 0.015410 *
## [Link](rank)2 -0.415399 0.194977 -2.131 0.033130 *
## [Link](rank)3 -0.812138 0.208358 -3.898 9.71e-05 ***
## [Link](rank)4 -0.935899 0.245272 -3.816 0.000136 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.41 on 394 degrees of freedom
## AIC: 470.41
##
## Number of Fisher Scoring iterations: 4
z_2=predict(probit,x_new)
exp(z_2)/(1+exp(z_2))

## 1
## 0.6260374
Como vemos los coeficientes de la estimación cambian y la predicción para la probabilidad de
admisión es algo inferior a la anterior, alrededor del 62%. 

Ejercicio 9 (propuesto).

Los datos hapiness disponibles en el paquete wooldridge contienen información sobre el

nivel de felicidad de individuos encuestados entre 1994 y 2006 (Wooldridge, 2013) para una
muestra de n = 17137. Las variables contenidas son vhappy (toma el valor 1 si el encuestado
declara ser muy feliz y 0 en caso contrario), educ (años de educación), income (ingresos
anuales), female (mujer), unem10 (desempleado). Obtener un modelo de regresión logística
que explique la probabilidad de ser muy feliz en función de las otras variables. Predecir la
probabilidad de ser muy feliz de una mujer con 15 años de educación con ingresos superiores
a 25000 dólares anuales (income="25000 or more")

34