• Author / Uploaded
• Manuel Hernandez

Citation preview

Taller Modelos Manuel Hern´andez Mart´ınez Maryoris Perez Sanchez Katy Pacheco Manchego Lili Bautista Arelllano Fernando Lagares Narvaes 8 de junio de 2017

1

2 ´ TALLER MODELOS DE REGRESION

Regresi´ on Lineal Simple 5. Que sucede con los estimadores del modelo de regresi´on lineal simple si en vez de trabajar con la variable independiente X se trabaja con la variable centrada, es decir, Xi − X, suponiendo que tambi´en Y es centrado, es decir, el modelo se trabaja con Yc = Y − Y : para el modelo de regresi´on lineal simple demuestre que: a)

n X

Ybi ei = 0

i=1

b) c)

n X i=1 n X

Xi ei = 0 X i Yi =

n X

i=1

Xi Ybi

i=1

6 Demuestre que para el modelo lineal simple a) r = βb1 b)

n X i=1

q

Yi =

Sxx Syy n X

Ybi

i=1

c) La estad´ıstica de prueba para la hip´otesis H0 : β1 = 0; con CM R Fc = CM es equivalente a E βb2 Fc = 2 σ bβb

1

11 En la siguiente tabla aparecen datos sobre el desempe˜ no de los 26 equipos de las Ligas nacionales de Futbol en 1976. Se cree que la cantidad de yardas ganada por tierra por los contrarios (X8 ) tiene un efecto sobre la cantidad de juegos que gana un equipo (y) a) Ajustar un modelo de regresi´on lineal multiple que relacione los juegos ganados, Y con la cantidad de yardas ganada por tierra por los contrarios (X8 )

3 b) Formar la tabla de an´alisis de varianza y probar el significado de la regresi´on. c) Determinar un intervalo de confianza de 95 % para la pendiente . d ) ¿Qu´e porcentaje de variabilidad total da y, y explica este modelo? e) Determinar un intervalo de confianza de 95 % para la cantidad promedio de juegos ganados, si la distancia ganada por tierra por los contrarios se limita a 2000 yardas Donde: Y Juegos ganados (por temporada de 14 juegos) x1 : Yardas por tierra (temporrada) x2 : Yardas por aire (temporada) x3 : Promedio de pateo (yardas/patadas) x4 : Porcentaje de goles (GC hechos/ GC intentados, temporada) x5 : Diferencia de perdidas de bal´on (perdidas ganadas/ perdidas) x6 : Yardas de castigo (temporadas) x7 : Porcentaje de carreras (jugadas por tierra/ jugadas totales) x8 : Yardas por tierra del contrario (temporada) x9 : Yardas por aire del contrario (temporada)

Regresi´ on M´ ultiple 0 En la siguiente tabla aparecen datos sobre el desempe˜ no de los 26 equipos de las Ligas nacionales de Futbol en 1976. Se cree que la cantidad de yardas ganada por tierra por los contrarios (X8 ), las Yardas por tierra (X1 ) y las Yardas por aire (X2 ) tienen un efecto sobre la cantidad de juegos que gana un equipo (y) a) Ajustar un modelo de regresi´on lineal multiple que relacione los juegos ganados, Y con la cantidad de yardas ganada por tierra por los contrarios (X8 ), las Yardas por tierra (X1 ) y las Yardas por aire (X2 ) b) Formar la tabla de an´alisis de varianza y probar el significado de la regresi´on. c) Determinar un intervalo de confianza de 95 % para los parametros β1 , β2 y β3 .

4

Cuadro 1: Desempe˜ no de equipos en ligas Nacionales de F´ utbol 1976

d ) ¿Qu´e porcentaje de variabilidad total da y, y explica este modelo? e) Determinar un intervalo de confianza de 95 % para la cantidad promedio de juegos ganados, si la distancia ganada por tierra por los contrarios se limita a 2000 yardas, las yardas por tierra son 1900 y las yardas por aire son 1950 5) Sea X ; Np (V µ, σ 2 V ), con V una matriz positiva sim´etrica con −1 ρ(V ) = P . Definamos B = V −1 − V −1 1p (1p 0 V −1 1p ) 1p V −1 a) Demuestre que B es sim´etrica y que BV es idempotente b) Encuentre la distribuci´on de Y = βX c) Encuentre la distribuci´on de Y 0 ΣY d ) Es Y y Y 0 ΣY iid

6) Considere el modelo Yi = β1 Xi1 + β2 Xi2 + εi ; εi ; N (0, σ 2 ), use los siguientes datos

5 Onservaci´on 1 2 3 4 5 6 Y 82 79 74 83 80 81 X1 10 9 9 11 11 10 X2 15 14 13 15 14 14

7 8 84 81 10 12 16 13

a) Estime β1 , β2 , y σ b) Realice el respectivo an´alisis de varianza y los an´alisis parciales y secuenciales c) Obtenga intervalos de confianza del 95 % para β1 y β1 + β2 d ) Encuentre un p − valor aproximado para las hip´otesis H0 : β1 − β2 = 0 e) Lleve a cabo el an´alisis de selecci´on de variables, usando el R2 , 2 Radj , Cp , PRESS, AIC, BIC, y la varianza residual. 7) Considere el modelo Yi = β1 Xi1 + β2 Xi2 + i ; i ; N (0, σ 2 ), use los datos del ejercicio 6) para obtener a) SCR(X1 , X2 |X0 ) b) Corroborar la hip´otesis H0 : β0 = 0 vs H1 : B0 6= 0 c) Lleve a cabo el ANAVA, los an´alisis parciales y secuenciales. 2 d ) Calcule las estad´ısticas de selecci´on de variables, R2 , Radj , Cp , PRESS, AIC, BIC, la varianza residual, etc.

9) Considere los siguientes datos (Ver Graybill, p´ag 227), que se asumen vienen del modelo Yi = β0 + β1 Xi1 + β2 Xi2 + εi Y X1 X2 Y X1 X2

con εi ; N (0, σ 2 )

12.1 5.5 4.6 4.5 10.8 4.9 0.870 0.202 0.203 0.198 0.730 0.510 1.69 1.17 1.17 1.21 1.63 1.59 6.0 4.2 5.3 6.7 4.0 6.1 0.205 0.670 0.205 0.271 0.203 0.264 1.14 1.92 1.22 1.71 1.16 1.37

6 a) Encuentre las estimativas de β0 , β1 y β2 b) Encuentre los errores est´andar de βb0 , βb1 y βb2 e interpretar estas estimativas c) Encuentre la estimaci´on puntual de 2β0 + β1 y β1 − β2 , adem´as los errores est´andar. Interprete los resultados d ) Lleve a cabo el respectivo an´alisis de varianza respectivo y encuentre el grado de ajuste. e) Lleve a cabo los an´alisis parciales y secuenciales. f ) Encuentre las estad´ısticas de selecci´on de variables, Cp , PRESS, AIC, BIC, la varianza residual, etc.

An´ alisis de Diagn´ osticos 9 Un punto se considera influyente si ´este es influyente para la SCE. Entonces para evaluar la influencia de un punto se puede analizar la estad´ıstica P (i) = SCE − SCE(i) con SCE(i) la suma de cuadrado del error despu´es de haber eliminado la i-´esima observaci´on ei a) Demuestre que P (i) = 1 − hii b) Encuentre la distribuci´on de P (i) c) Demuestre que r

(1 − hii )P (i) ∼ tglee CM E

d ) Establezca el criterio para clasificar un punto como influyente. 10 Para los datos en adjunto (Cuadro 1) tomados del libro de Draper and Smith, pagina 111-113. a) Estime el modelo de regrsi´on lineal simple b) Realice el an´alisis de diagn´ostico respectivo c) Realice los gr´aficos de ei ∼ Ybi y ei ∼ Xi

7 d ) Con los datos de X que se repiten, realice el gr´afico X i vs Si2 e) Estime la relaci´on 2

Si2 = α0 + α1 X i + α2 X i

f ) Para i = 1, ..,28 reemplace cada Xi en la ecuaci´on encontrada en cc = 1 v) y obtener W Si2 g) Defina P = diag{W −1/2 } y estime nuevamente el modelo con estos pasos. c 1/2 (Yi − Ybi ) = Yip − Ybip , donde c 1/2 Ybi y ebi = W h) viii) Defina Ybip = W i c 1/2 Yi , realice nuevamente el los valores observados son Yip = W an´alisis de diagnostico y las pruebas de validaci´on. X 1.15 1.90 3.00 3.00 3.00 3.00 3.00 5.34 5.38 5.40 5.40 5.45 7.70 7.80 7.81 7.85 7.87

Y 0.99 0.98 2.60 2.67 2.66 2.78 2.80 5.92 5.35 4.33 4.89 5.21 7.68 9.81 6.52 9.71 9.82

ci W 1.24028 2.18224 7.84930 7.84930 7.84930 7.84930 7.84930 7.43652 6.99309 6.78574 6.78574 6.30514 0.89204 0.84420 0.83963 0.82171 0.81296

X 7.94 9.03 9.07 9.11 9.14 9.16 9.37 10.17 10.18 10.22 10.22 10.22 10.18 10.50 10.23 10.03 10.23

Y 8.50 9.47 11.45 12.14 11.50 10.65 10.64 9.78 12.39 11.03 8.0 11.90 8.68 7.25 13.46 10.19 9.93

ci W 0.78342 0.47385 0.46621 0.45878 0.45327 0.44968 0.41435 0.3118 0.3107 0.3067 0.3067 0.3067 0.3107 0.2803 0.3057 0.3268 0.3057

Cuadro 2: Datos libro de Draper and Smith

13 Demuestre que si centran las variables en el modelo de regresi´on lineal, entonces

8

H=

1 J + Hc n

donde Hc es la matriz hatcon los datos centrados. 14 De acuerdo a los resultados del ejercicio 13) verifique que si hik es un elemento de H y hcik es un elemento de Hc entonces Ybi =

X

hik Yi = Y +

X

hcik Yi

i

i

17 Lleve a cabo el an´alisis de diagn´ostico y las pruebas de validaci´on del ejercicio 7) y 10) del cap´ıtulo de modelos de regresi´on lineal multiple.

Modelos Polin´ omicos 1 Ajuste el modelo de regresi´on polinomica 2 2 Yˆi = βˆ0 + βˆ1 Xi1 + βˆ2 Xi2 + βˆ11 Xi1 + βˆ22 Xi2 + βˆ12 Xi1 X2i

Con los siguientes datos, y 26 24 175 160 163 55 62 100 26 30 70 71

x1 1.0 1.0 1.5 1.5 1.5 0.5 1.5 0.5 1.0 0.5 1.0 0.5

x2 1.0 1.0 4.0 4.0 4.0 2.0 2.0 3.0 1.5 1.5 2.5 2.5

Cuadro 3:

9 3 En un experimento de ingener´ıa qu´ımica que se refiere a la transferencia de calor en una cama fluidizada superficial, se recolectan los datos de las siguientes cuatro variables de regresi´on: tasa del flujo de gas fluidizante, Ib/hr(x1 ), la tasa del flujo del gas flotante Ib/hr(x2 ); abertura de la entrada de gas flotante, mil´ımetros (x3 ); temperatura de entrada del gas flotante, 0 F (x4 ). Las respuestas medidas son: eficiencia de la transferencia de calor (y1 ); la efiencia t´ermica(y2 ). Los datos son los siguientes:

Observaci´on

y1

y2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

41.852 155.329 99.628 49.409 72.958 107.702 97.239 105.856 99.348 111.907 100.008 175.380 117.800 217.409 41.725 151.139 220.630 131.666 80.537 152.966

38.75 51.87 53.79 53.84 49.17 47.61 64.19 52.73 51.00 47.37 43.18 71.23 49.30 50.87 54.44 47.93 42.91 66.60 64.94 43.18

x1

x2

x3

69.69 170.83 45 113.46 230.06 25 113.54 228.19 65 118.75 117.73 65 119.72 117.69 25 168.38 173.46 45 169.85 169.85 45 169.85 170.86 45 170.89 173-92 80 71.31 173.34 25 171.43 171.43 45 171.59 263.49 45 171.63 171.63 45 171.93 170.91 10 173.92 71.73 45 221.44 217.39 65 222.74 221.73 25 228.90 114.40 25 231.19 113.52 65 236.84 167.77 45

x4 219.74 181.22 179.06 281.30 282.20 216.14 223.88 222.80 218.84 218.2 219.20 168.62 217.58 219.92 296.60 189.14 186.08 285.80 286.34 221.72

10 Considere el modelo para pronosticar la respuesta coefiente de transferencia de calor 4 4 X X XX y1i = β0 + βj xji + βjj x2ji + j=1

jF) x8 1 178.09 178.092 31.103 7.381e-06 *** Residuals 26 148.87 5.726 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Con la informaci´on suministrada en la anterior tabla y los valores estimados a partir del modelo, se obtiene que Fc = 31,103; entonces F c > F(0,05;1;26) = 4,225; lo cual ´ındica que el modelo lineal si explica la cantidad de juegos ganados por temporada a partir la distancia ganada por tierra por los contrarios. c) Intervalos de confianza al 95 % para los par´ametros del modelo son

25 > confint(modelo,level=0.95) 2.5 % 97.5 % (Intercept) 16.246064040 27.330437725 x8 -0.009614347 -0.004435854 obteni´endose el intervalo -0.00961 y -0.00444 para la pendiente, es decir, que con un 95 % de confiabilidad, en el rango de -0.00961 y -0.00444 puede disminuir la cantidad de juegos ganados por temporada. d ) Para conocer que porcentaje de variabilidad total de Y , explica este modelo se an´aliza el valor del coeficiente de determinaci´on R2 , Residual standard error: 2.393 on 26 degrees of freedom Multiple R-squared: 0.5447,Adjusted R-squared: 0.5272 F-statistic: 31.1 on 1 and 26 DF, p-value: 7.381e-06 de donde, la cantidad de juegos ganados por temporada es explicada en un 54,47 % la distancia ganada por tierra por los contrarios, lo cual llevaria a pesar intuitiva e inferencialmente que el modelo no tiene un buen ajuste. e) Una estimaci´on y un intervalo de confianza al 95 % para la cantidad promedio de juegos ganados, si la distancia ganada por tierra por los contrarios se limita a 2000 yardas > predict(modelo,newdata=data.frame(x8=2000),interval="confidence") fit lwr upr 7.73805 6.765753 8.710348 Por tanto, con una confianza de 95 % una estimaci´on para el valor real de la cantidad promedio de juegos ganados, si la distancia ganada por tierra por los contrarios se limita a 2000 yardas es de 8 partidos y pudiendo variar entre 7 y 9 partidos.

Regresi´ on M´ ultiple 0.

a) Tenemos el modelo lineal multiple Yi = β0 + β1 Xi1 + β2 Xi2 + β3 Xi8 + εi

26 por tanto ajustando el modelo por M.C.O tendr´ıamos el modelo. Yˆi = βˆ0 + βˆ1 Xi1 + βˆ2 Xi2 + βˆ3 Xi8 el cual lo ajust´andolo en R nos da los siguientes resultados > mod2 mod2 Call: lm(formula = y ~ X1 + X2 + X8) Coefficients: (Intercept) 6.464714

X1 0.002229

X2 0.003232

X8 -0.005251

es decir el modelo ajustado es : Yˆi = 6,4647140 + 0,002229Xi1 + 0,003232Xi2 + −0,00525Xi8 b) Se desea probar la hip´otesis Ho : β = 0 contra β 6= 0 para esto se realizan los respectivos an´alisis de varianza. Note que  0   6,464714 195     ˆ 0 Y =  0,002229  432994 = 1612,266 SCR = βX  0,003232  437445 −0,005251 386127 y SCT = Y0 Y = 1685 de donde SCE = SCT − SCR = 1685 − 1612,266 = 72,734 as´ı obtenemos los resultados de la tabla?? de an´alisis de varianza.

27 F.V. Regresi´on Error Total

gl SC CM Fc F5 % 4 1612.266 403.0665 132.99964 2.77629 24 72.734 3.03058 28 1685

Cuadro 8: Medias de las variables analizadas El valor de la estad´ıstico Fc =132.99964 es mucho mayor que el tabulado a un nivel de significancia α= 0.05, por tanto H0 : β = 0 es rechazada, concluy´endose que al menos una de las variables explicativas es significativas en el modelo. c) Queremos hallar un intervalo de confianza de 95 % para β1 , β2 y β3 entonces utilizamos la siguiente c´odigo en R nos da: > confint(mod2)[2:4,] 2.5 % 97.5 % X1 -0.0001636138 0.004621932 X2 0.0018392027 0.004625779 X8 -0.0077770757 -0.002725183 por tanto tenemos; −0,0001636138 ≤β1 ≤ 0,004621932 0,0018392027 ≤β2 ≤ 0,004625779 −0,0077770757 ≤β3 ≤ −0,002725183 Concluy´endose que a un nivel de confianza del 95 % tenemos que β1 es estadisticamente igual a cero, ya que el intervalo de este par´ametro contiene a cero, mientras que β2 y β3 son estadisticamente distintas de cero. d ) la proporci´on de la variable respuesta Y , que es explicada por las variables X1 , X2 y X3 respectivamente, en conjunto; para esto se calcula el coeficiente de determinaci´on m´ ultiple, cuyo valor para este caso particular es > summary(mod2)$r.squared [1] 0.7775477

28 de donde se concluye que el ajuste del modelo es del 77, 75 %, siendo casi un 78 % lo que X1 , X2 y X3 logran explicar la variable respuesta. Una medida mas realista para el valor del del grado de ajuste es el coeficiente de determinaci´on m´ ultiple ajustado Ra2 , cuyo valor en este caso es > summary(mod2)$adj.r.squared [1] 0.7497411 e) para calcular un intervalo de confianza de 95 % para la cantidad promedio de juegos ganados, si la distancia ganada por tierra por los contrarios se limita a 2000 yardas, las yardas por tierra son 1900 y las yardas por aire son 1950 usamos el siguiente c´odigo en R de donde tenemos > predict(mod2,data.frame(X1=1900,X2=1950,X8=2000) + ,interval="prediction") fit lwr upr 1 6.501215 2.761713 10.24072 por tanto para este caso tenemos 2,761713 ≤ Y ≤ 10,24072 Por tanto, cuando las variables X1 , X2 y X8 toman valores de 1900, 1950 y 2000 se estima que la cantidad de juegos que gana un equipo son aproximadamente 7 y con una confianza del 95 % se concluye que el valor real de dicha cantidad estimada se encuentra entre 3 y 10 aproximadamente, para este caso se redondearon los resultados ya que la cantidad de juegos que gana un equipo deben ser valores enteros. 5 como tenemos X ; Np (V µ, σ 2 V ), con V una matriz positiva sim´etrica −1 con ρ(V ) = P y definiendo B = V −1 − V −1 1p (1p 0 V −1 1p ) 1p 0 V −1 a) Demostraremos que B es sim´etrica y que BV es idempotente Notemos que B, es simetrica, cuando se cuando que B 0 = B. En efecto,

29

 0 −1 B 0 = V −1 − V −1 1p (10p V −1 1p ) 10p V −1 0
 −1 −1 0 −1 0 −1 −1 0 − V 1p (1p V 1p ) 1p V = V 0

−1 0

0

0

= V −1 − (V −1 ) (10p ) ((10p V −1 1p ) ) (10p )(V −1 ) −1

= V −1 − V −1 1p (1p 0 V −1 1p ) 1p 0 V −1 Por lo tanto B es simetrica. Veamos que BV es idempotente. En efecto,   −1 0 −1 −1 0 −1 −1 (BV ) = V V − V 1p (1p V 1p ) 1p V V   −1 V −1 V − V −1 1p (10p V −1 1p ) 10p V −1 V   −1 = I − V −1 1p (10p V −1 1p ) 10p I   −1 0 −1 0 −1 I − V 1p (1p V 1p ) 1p I 2

−1

−1

=I − V −1 1p (10p V −1 1p ) 10p − V −1 1p (10p V −1 1p ) 10p + −1

−1

V −1 1p (10p V −1 1p ) 10p V −1 1p (10p V −1 1p ) 10p −1

−1

=I − V −1 1p (10p V −1 1p ) 10p − V −1 1p (10p V −1 1p ) 10p + −1

V −1 1p (10p V −1 1p ) 10p −1

=I − V −1 1p (10p V −1 1p ) 10p =I − A −1

con A = V −1 1p (10p V −1 1p ) 10p . Pero como A es idempotente, pues    −1 −1 A2 = V −1 1p (10p V −1 1p ) 10p V −1 1p (10p V −1 1p ) 10p
−1 −1 = V −1 1p (10p V −1 1p ) 10p V −1 1p (10p V −1 1p ) 10p −1

= V −1 1p (10p V −1 1p ) 10p =A As´ı A es idempotente e I es idempotente, de modo que, A − I es idempotente. Por lo tanto BV es idempotente.

30 b) Encontremos la distribuci´on de Y = BX Sabemos que E(X) = V µ y V ar(X) = σ 2 V . Luego, E(BX) = BE(X) = BV µ   −1 = V −1 V − V −1 1p (10p V −1 1p ) 10p V −1 V µ   −1 = I − V −1 1p (10p V −1 1p ) 10p µ = Cµ con C = I − A y V ar(BX) =B V ar(X)B 0 =B σ 2 V B =σ 2 BV B   −1 0 2 −1 0 −1 =σ I − V 1p (1p V 1p ) 1p   −1 V −1 − V −1 1p (1p 0 V −1 1p ) 10p V −1  −1 =σ 2 V −1 − V −1 1p (1p 0 V −1 1p ) 10p V −1 −1

− V −1 1p (1p 0 V −1 1p ) 10p V −1 −1

−1

+V −1 1p (1p 0 V −1 1p ) 10p V −1 1p (1p 0 V −1 1p ) 10p V −1   −1 =σ 2 V −1 − V −1 1p (10p V −1 1p ) 10p V −1



=σ 2 CV −1 As´ı como X ∼ Np (V µ, σ 2 V ) y Y es una combinaci´on lineal de X tenemos que Y ∼ Np (Cµ, σ 2 CV −1 ) −1

con C = I − V −1 1p (10p V −1 1p ) 10p 7.

a) Por definici´on sabemos que SCR(X1 , X2 |X0 ) = SCR(X1 , X2 , X0 ) − SCR(X0 )

31 tal que SCR(X1 , X2 , X0 ) corresponde a la suma de cuadrado de regresi´on para el modelo lineal Yi = β0 Xi0 + β1 Xi1 + β2 Xi2 + i mientras que SCR(X0 )es la suma de cuadrado de regresi´on para el modelo lineal Yi = β0 Xi0 +i Entonces, considerando los siguientes vectores de datos Y=c(82, 79, 74, 83, 80, 81, 84, 81) X0=c(rep(1,8)) X1=c(10, 9 ,9 ,11, 11, 10, 10, 12) X2=c(15, 14, 13, 15, 14, 14, 16, 13) de donde se sigue la siguiente matriz de datos x B_012=solve(t(x)%*%x)%*%t(x)%*%Y; B_012 [,1] X0 30.000 X1 1.625 X2 2.375 

ˆ es decir β 012 luego,

 30,000 =  1,625  2,375

32 > SCR_012=t(B_012)%*%t(x)%*%Y;SCR_012 [,1] [1,] 51900.25 as´ı, SCR(X0 , X1 , X2 ) = 51900,25 Ahora, para el modelo Yi = β0 Xi0 + i que el parametro β0 es > B_0=solve(t(X0)%*%X0)%*%t(X0)%*%Y;B_0 [,1] [1,] 80.5 es decir βˆ0 = 80,5 luego, > SCR_0=t(B_0)%*%t(X0)%*%Y;SCR_0 [,1] [1,] 51842 es decir, SCR(X0 ) = 51842 Por tanto, SCR(X1 , X2 |X0 ) = SCR(X1 , X2 , X0 ) − SCR(X0 ) = 51900,25 − 51842 = 58,25 item Para corrobora la hip´otesis H0 : β0 = 0 contra H1 : β0 6= 0 con β0 del modelo lineal Yi = β0 Xi0 + β1 Xi1 + β2 Xi2 + i ; considere la suma de cuadrados de la media y la suma de cuadrados de la regresi´on corregida por la media son, respectivamente 2

SCM = nY = 8(80,5)2 = 51842, y SCRm = SCR − SCM = 51900,25 − 51842 = 58,25 Por otro lado la suma de cuadrado total es

33 b) Para corrobora la hip´otesis H0 : β0 = 0 contra H1 : β0 6= 0 con β0 del modelo lineal Yi = β0 Xi0 + β1 Xi1 + β2 Xi2 + i ; considere la suma de cuadrados de la media y la suma de cuadrados de la regresi´on corregida por la media son, respectivamente 2

SCM = nY = 8(80,5)2 = 51842, y SCRm = SCR − SCM = 51900,25 − 51842 = 58,25 Por otro lado la suma de cuadrado total es > SCT=t(Y)%*%Y;SCT [,1] [1,] 51908 as´ı, SCE = SCT − SCM − SCRm = 51908 − 51842 − 58,25 = 7,75. Los resultados anteriores llevan a la siguiente tabla de an´alisis de varianza

F.V. gl SC CM Fc F5 % Media 1 51842 51842 33446.45 6.608 Regresi´on corr 2 58.25 37.58065 6.988 5.7861 Error 5 7.75 1.55 Total 8 51908 Suponiendo que el modelo no tiene ninguna variable explicativa, es decir β1 y β2 son iguales a cero. Esto es Yi = β0 + i , entonces

34 2 βˆ0 = Y y la suma de cuadrados del modelo ser´ıa nY . En tal caso E(Y ) = β0 . Por tanto, bajo la hip´otesis nula y dado que el valor calculado de la estad´ıstica F (M ) = 33446,45 es mucho mayor que el tabulado a un nivel de significancia de α = 0,05; rechazando H0 , se concluye que el par´ametro β0 es significativo para el modelo.

c)

• An´ alisis de varianza Para probar el significado de la  regresi´  on, es decir verificar la β1 0 hip´otesis nula H0 : β 12 = = considere los siguienβ2 0 tes resultados de la tabla ANAVA de la regresi´on Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value F 5% regresi´ on 2 58.25 29.125 18.79032 5.7861 Residuals 5 7.750 1.550 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 As´ı, con una confianza de 95 % dada la informaci´on suministrada en la anterior tabla y los valores estimados a partir del modelo, se obtiene que el valor Fc = 18,79032, es mayor que el tabulado F(0,05;2;5) = 5,7861; lo que lleva al rechazo de la hip´otesis nula, es decir el modelo lineal logra explicar la variable Y por medio de las variables X1 y X2 . • An´ alisis de varianza parcial Ahora el interes radica en verificar cual de las variables explicativa logra explicar a la variable Y , para esto con la funci´on summary obtenemos los valos de las estadisticas para verificar las siguientes hip´otesis H0 : β0 = 0 contra H1 : β0 6= 0 H0 : β1 = 0 contra H1 : β1 6= 0 H0 : β2 = 0 contra H1 : β2 6= 0 En efecto, > summary(mod) Call: lm(formula = Y ~ X1 + X2)

35 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.0000 8.2583 3.633 0.01502 * X1 1.6250 0.4556 3.567 0.01610 * X2 2.3750 0.4556 5.213 0.00343 ** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Dado que el p-valor de la estad´ıstica de cada respectiva prueba es menor que el nivel de significancia α = 0,05, se rechazan todas las tres hip´otesis nula; concluyendo con una confianza de 95 % que las dos variables explicativas aportan a la explicaci´on de la variable Y . • An´ alisis secuencial A pesar que el an´alisis parcial para todos los par´ametro resultaron significativo, podemos reafirmar este hecho con el siguiente an´alisis, > step(mod,direction="backward") Start: AIC=5.75 Y ~ X1 + X2 Df Sum of Sq

- X1 - X2

1 1

RSS AIC 7.750 5.746 19.717 27.467 13.868 42.117 49.867 18.639

Call: lm(formula = Y ~ X1 + X2) Coefficients: (Intercept) X1 X2 30.000 1.625 2.375 La regresi´on por pasos nos recomienda que el mejor modelo es Y en funci´on de X1 y X2 , ya que con este se obtiene el menor valor para el AIC. d ) A continuaci´on se calculan estad´ısticas de selecci´on de variables, R2

36 $which 1 2 1 FALSE TRUE 1 TRUE FALSE 2 TRUE TRUE

solo con X_2 solo con X_1 con X_1 y X_2

$label [1] "(Intercept)" "1"

"2"

$size [1] 2 2 3 $r2 [1] 0.5838384 0.2444444 0.8825758 R2ajt $which 1 2 1 FALSE TRUE 1 TRUE FALSE 2 TRUE TRUE

solo con X_2 solo con X_1 con X_1 y X_2

$label [1] "(Intercept)" "1"

"2"

$size [1] 2 2 3 $adjr2 [1] 0.5144781 0.1185185 0.8356061 Cp $which 1 2 1 FALSE TRUE 1 TRUE FALSE 2 TRUE TRUE

solo con X_2 solo con X_1 con X_1 y X_2

37 $label [1] "(Intercept)" "1"

"2"

$size [1] 2 2 3 $Cp [1] 13.72043 28.17204

3.00000

Adem´as, Modelo solo con X1 solo con X2 con X1 y X2

AIC 43.34231 38.57127 30.44903

PRESS 101.4691 59.99297 25.60019

BIC 38.80959 43.58063 30.76679

En conclusi´on, cada uno de los criterios de selecci´on de variables coincide en que el modelo que mejor logra explicar a la variable dependiente Y , es el modelo con las dos variables independientes X1 y X2 , es decir este es el modelo que mejor se ajusta a los datos explicando la mayor variabilidad total. 6. Haciendo los calculos con R tenemos reg|t|) (Intercept) 30.0000 8.2583 3.633 0.01502 x1 1.6250 0.4556 3.567 0.01610 x2 2.3750 0.4556 5.213 0.00343 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05

8 0.625

* * ** ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.245 on 5 degrees of freedom Multiple R-squared: 0.8826,Adjusted R-squared: 0.8356 F-statistic: 18.79 on 2 and 5 DF, p-value: 0.004725 El p-valor de la estad´ıstica F c = 18,79 es menor al nivel de significancia α, tambi´en se tiene (p − valor = 18,79 < 0,05) por tanto H0 : β = 0

40 no es rechazado, concluy´endose que las variables explicativas no logran explicar la variable respuesta. El an´alisis parcial arroja que la variables X1 y X2 es significativa en el modelo. Los siguientes errores est´andar 8.2583, 0.4556, 0.4556 para βb0 , βb1 y βb2 resrespectivamente, para obtener S usamos summary(reg)$sigma summary(mod1)$sigma 0.1781172 Obteniendo σ b = 0,1781172 Sumas de cuadrados parciales Para obtener las sumas de cuadrados parciales de la regresi´on, inicialmente se ajusta el modelo completo, esto es b i = 30 + 1,625Xi1 + 2,375Xi2 Y

(1)

Implementando en R el an´alisis secuencial tenemos step(reg,direction="backward") Start: AIC=5.75 y ~ x1 + x2 Df Sum of Sq

- x1 - x2

1 1

RSS AIC 7.750 5.746 19.717 27.467 13.868 42.117 49.867 18.639

Call: lm(formula = y ~ x1 + x2, x = T) Coefficients: (Intercept) 30 1.625

x1

x2 2.375

note que el mejor modelo es con las dos variable ya que tiene un AIC menor, observe que el modelo sin la variable X1 tiene un AIC=18.639

41 y sin la variable X2 tiene un AIC=13.868. Calculamos ahora un intervalo del 95 % de confianza para β1 y β1 + β2 en el caso de datos As´ı, un intervalo al 95 % de confianza para β1 , es √ √ βb1 − t(8−3,0,05/2) CM Ea11 ≤ β1 ≤ βb1 + t(8−3,0,05/2) CM Ea11 ,entonces p p 1,625 − 2,5706 (10381,6)(0,134) ≤ β1 ≤ 1,625 + 2,5706 (10381,6)(0,134) −94,22679 ≤ β1 ≤ 97,47679 De lo anterior se concluye que a un nivel confianza del 95 % tenemos que β1 es estad´ısticamente igual a cero, ya que el intervalo de este par´ametro contiene a cero. ahora un intervalo del 95 % de confianza para β1 + β2 q q −1 0 0ˆ 0ˆ 0 0 α α K β−t(n−p, 2 ) CM E K (X X) K ≤ K β ≤ K β−t(n−p, 2 ) CM E K0 (X0 X)−1 K 

 30
b = 1,625 el intervalo obtenido es el donde K 0 = 0 1 1 y β 2,375 siquiente, −136,0005 ≤ β1 + β2 ≤ 144,0005 9

a) para hallar las estimativas de β0 , β1 y β2 ajustamos el modelo en R y mediante el siguiente c´odigo se tiene: > mod mod$coefficients (Intercept) X1 7.888132 11.073717

X2 -4.130306

por tanto tenemos que βˆ0 = 7,888, βˆ1 = 11,074 y βˆ2 = −4,130 b) Los errores est´andar para βˆ0 , βˆ1 y βˆ2 los calculamos con el siguiente c´odigo en R: > summary(mod)$coefficients[,2] (Intercept) X1 X2 3.821288 3.798024 3.416633

42 es decir tenemos que σ ˆβˆ0 = 3,821, σ ˆβˆ1 = 3,798 y σ ˆβˆ2 = 3,417 Para la interpretaci´on de los par´ametros estimados, se puede concluir que con un error de 3.798 por cada unidad que se incremente X1 manteniendo fijo X2 la variable Y incrementa en 11.074 mientras que con un error de 3.417 si se mantiene fijo X1 a variable Y disminuye -4.130, ademas con un error del 3.821 se estima que la cantidad inicial bajo ausencia de X1 y X2 es 7.888 c) La estimaci´on para V ar(β) mediante R es: > X=cbind(1,X1,X2) > XtX solve(XtX)*summary(mod)$sigma^2 X1 X2 14.602240 9.275843 -12.58924 X1 9.275843 14.424983 -10.40458 X2 -12.589240 -10.404577 11.67338 es decir tenemos que: 

 14,602 9,276 −12,589 ˆ =  9,276 14,425 −10,405 V\ ar(β) −12,589 −10,405 11,673 Cabe notar, que cuando se incrementa una unidad de la cantidad X1 y teniendo en cuenta 2β0 , tenemos que: 2βˆ0 + βˆ1 = 2(7,888) + 11,074 = 15,776 + 11,074 = 26,85 y σ ˆ22βˆ0 +βˆ1 = 22 · σ ˆβ2ˆ0 + σ ˆβ2ˆ1 + σ ˆβ2ˆ0 βˆ1 = 4(3, 821)2 + (3, 417)2 + 2(9,276) = 4(14, 602) + 14, 425 + 2(9,276) = 29,204 + 14, 425 + 18, 552 = 62,181 entonces σ ˆ2βˆ0 +βˆ1 =

q p σ ˆ22βˆ +βˆ = 62,181 = 7,8855 0

1

43 Luego, con un error de 7.8855 se estima que un aumento de una unidad X1 y teniendo en cuenta el intercepto 2β0 se tiene una estimaci´on de 26.85 Ahora si se incrementa una unidad de X1 y se disminuyen dos unidades de X2 tenemos: βˆ1 + βˆ2 = 11,074 − (−4,130) = 11,074 + 4,130 = 15,204 y σ ˆβ2ˆ1 −βˆ2 = σ ˆβ2ˆ1 + σ ˆβ2ˆ2 − σ ˆβ2ˆ0 βˆ1 = (3, 798)2 + (3, 417)2 − 2(−12,589) = 14,425 + 11,676 + 25,178 = 51,279 entonces σ ˆ2βˆ0 +βˆ1 =

q p σ ˆ22βˆ +βˆ = 51,279 = 7,1609 0

1

Es decir con un error del 7.1609 se estima que un aumento de una unidad de X1 y si se disminuye una unidad de X2 , entonces se estima un aumento en 15,204. d ) Se desea probar la hip´otesis H0 : β = 0 contra β 6= 0 para esto se realizan los respectivos an´alisis de varianza. Note que  0   7,888 74,700 ˆ 0 Y =  11,074   33,214  = 507,136 SCR = βX −4,130 108,930 y SCT = Y0 Y = 538,55 de donde SCE = SCT − SCR = 538,55 − 507,136 = 31,414

44 F.V. Regresi´on Error Total

gl SC CM Fc F5 % 3 507.136 169.045 48.437 3.862548 9 31.414 3.49 12 538.55

Cuadro 10: Medias de las variables analizadas as´ı obtenemos los resultados de la tabla?? de an´alisis de varianza. El valor de la estad´ıstico Fc =48.437 es mucho mayor que el tabulado a un nivel de significancia α= 0.05, por tanto H0 : β = 0 es rechazada, concluy´endose que al menos una de las variables explicativas es significativas en el modelo. la proporci´on de la variable respuesta Y , que es explicada por las variables X1 y X2 respectivamente, en conjunto; para esto se calcula el coeficiente de determinaci´on m´ ultiple, cuyo valor para este caso particular es > summary(mod)$r.squared [1] 0.5728469 de donde se concluye que el ajuste del modelo es del 57,28 %, siendo casi un 57 % lo que X1 y X2 logran explicar la variable respuesta. Una medida mas realista para el valor del del grado de ajuste es el coeficiente de determinaci´on m´ ultiple ajustado Ra2 , cuyo valor en este caso es > summary(mod)$adj.r.squared [1] 0.477924 e) Dado que la hip´otesis nula H0 : β = 0 fue rechazada, entonces al menos un subconjunto de las 2 variables explicativas logran explicar la variable respuesta, por lo tanto es de inter´es conocer que variables son las que logran explicar la variable respuesta y cuales no. En estos casos se lleva a cabo un an´alisis parcial, es decir, una prueba de hip´otesis para cada variable mediante el test H0 : βj = 0 contra H1 : βj 6= 0 con j = 0, 1, 2; luego, bajo la

45 hip´otesis nula se utilizar´ıa el estad´ıstico de prueba tj = √

βˆj CM Eajj

con ajj el j −esimo elemento de la diagonal de la matriz (X0 X)−1 . En efecto, Para j = 0 se tiene la hip´otesis H0 : β0 = 0 contra H1 : β0 6= 0 y adem´as que a00 = 4,184 y βˆ0 = 7,888. Luego, el estad´ıstico de prueba es 7,888 βˆ0 t0 = √ =p = 2,064 (3,49)(4, 184) CM Ea00 Para j = 1 se tiene la hip´otesis H0 : β1 = 0 contra H1 : β1 6= 0 y adem´as que a11 = 4,133 y βˆ1 = 11,074. Luego, el estad´ıstico de prueba es: t1 = √

βˆ1 CM Ea11

11,074 = 2,916 =p (3,49)(4,133)

Para j = 2 se tiene la hip´otesis H0 : β1 = 0 contra H1 : β2 6= 0 y adem´as que a22 = 3,344 y βˆ2 = −4,130. Luego, el estad´ıstico de prueba es: βˆ1 −4,130 t1 = √ =p = −1,209 (3,49)(3,344) CM Ea22 Al comparar los estad´ısticos de prueba antes hallados con t(1 − α/2, n − p) = t(0,975, 10) = 2,228, se rechazan las hip´otesis nulas H0 : β1 = 0 dado que los valores absolutos de las estad´ıstica de prueba, |t1 |, es mayor que el valor tabulado; mientras que no se rechaza las hip´otesis nula H0 : β0 = 0 y H0 : β2 = 0 ya que |t0 | y |t2 | son menor que el valor tabulado. En conclusi´on, con un nivel de confianza del 95 % se tiene que la variable que logra explicar significativamente a la variable respuesta Y es la variable X1 , es decir el aporte de X2 al modelo

46 estimado es muy pobre. Para el an´alisis secuencial implementaremos R: > step(mod,direction="backward") Start: AIC=17.55 Y ~ X1 + X2 Df Sum of Sq RSS AIC - X2 1 5.1009 36.515 17.354

31.414 17.548 - X1 1 29.6723 61.086 23.529 Step: AIC=17.35 Y ~ X1 Df Sum of Sq

- X1

1

RSS AIC 36.515 17.354 37.028 73.543 23.756

Call: lm(formula = Y ~ X1) Coefficients: (Intercept) 3.434

X1 7.392

El cual nos dice que el mejor modelo es Yˆ = β0 + βˆ1 X1 Pero dado que AIC disminuye muy poco no se recomienda eliminar la variable X2 f ) Ahora encontramos mediante R las estad´ısticas de selecci´on de variables, Cp, PRESS, AIC, BIC, la varianza residual. > require(leaps) > leaps(x=Xd, y=Y, method = "r2") $which 1 2 1 TRUE FALSE

47 Modelo solo con X1 solo con X2 con X1 y X2

PRESS 64.10259 105.9383 109.819

AIC BIC 53.40826 54.86298 59.58307 61.03779 53.60265 55.54228

Cuadro 11: Resumen de Selecci´on de variables 1 FALSE 2 TRUE

TRUE TRUE

$label [1] "(Intercept)" "1"

"2"

$size [1] 2 2 3 $r2 [1] 0.5034870 0.1693757 0.5728469 > leaps(x=Xd, y=Y,method = "Cp" ) $which 1 2 1 TRUE FALSE 1 FALSE TRUE 2 TRUE TRUE $label [1] "(Intercept)" "1"

"2"

$size [1] 2 2 3 $Cp [1] 2.461395 9.501029 3.000000 En conclusi´on, cada uno de los criterios de selecci´on de variables coincide en que el modelo que mejor logra explicar a la variable dependiente Y , es el modelo con la variables independiente X1 , sin

48 embargo los resultados con respecto al modelo con las variables X1 y X2 no var´ıan mucho por lo que se recomienda tomar este modelo es decir este es el modelo que mejor se ajusta a los datos explicando la mayor variabilidad total.

An´ alisis de Diagn´ osticos 9

a) En efecto se tiene que Pi = SCE − SCE(i) X e2k(i) = SCE − k6=i

"



ei 1 − hii # " k=1 n X e2i 2 = SCE − ek − 1 − hii k=1   e2i = SCE − SCE − 1 − hii n X e2i 2 = SCE − ek + 1 − hii k=1 = SCE −

=

10.

X

e2k −

e2i 1 − hii

a) tenemos el modelo de regresi´on lineal simple Y = β0 + β1 X 1 + εi El modelo ajustado por MCO es: Yb = −0,579 + 1,135X1 En R,

#

2 (1 − hii )

49 > mod influence.measures(mod) Influence measures of lm(formula = Dep.Y ~ Indep.X, data = dia) :

1 2 3 4 5 6 7 8 9 10 11 12 13

dfb.1_ dfb.In.X dffit cov.r cook.d hat inf 0.085691 -0.077960 0.085852 1.267 3.80e-03 0.1629 * -0.169005 0.150599 -0.170077 1.211 1.48e-02 0.1323 * -0.051062 0.043627 -0.052218 1.173 1.40e-03 0.0946 -0.035326 0.030183 -0.036126 1.174 6.73e-04 0.0946 -0.037574 0.032103 -0.038425 1.174 7.61e-04 0.0946 -0.010614 0.009069 -0.010854 1.175 6.07e-05 0.0946 -0.006122 0.005231 -0.006261 1.175 2.02e-05 0.0946 0.052925 -0.036521 0.063627 1.105 2.08e-03 0.0426 -0.021478 0.014714 -0.025987 1.109 3.48e-04 0.0421 -0.146849 0.100230 -0.178255 1.061 1.60e-02 0.0418 -0.078933 0.053875 -0.095815 1.095 4.70e-03 0.0418 -0.046648 0.031536 -0.057104 1.104 1.68e-03 0.0411 -0.014003 -0.007154 -0.057466 1.088 1.70e-03 0.0290

50 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

0.039106 -0.044697 0.032310 0.033284 0.029793 0.001258 0.004566 -0.042436 -0.063364 -0.046448 -0.022897 -0.020798 0.081544 -0.097834 -0.000359 0.227978 -0.061571 0.163760 0.386768 -0.178673 0.038613 0.078564

0.029472 -0.034956 0.029248 0.032391 0.033489 0.001576 -0.014321 0.127436 0.182841 0.130356 0.063139 0.049332 -0.144081 0.172470 0.000627 -0.398340 0.107581 -0.288690 -0.640126 0.311513 -0.070603 -0.136975

0.185998 -0.216016 0.166907 0.178062 0.171584 0.007684 -0.028140 0.246124 0.347303 0.244630 0.117557 0.085193 -0.204750 0.244658 0.000883 -0.561130 0.151547 -0.409522 -0.862996 0.438066 -0.103003 -0.192622

1.021 0.997 1.036 1.028 1.033 1.096 1.105 1.010 0.929 1.014 1.085 1.100 1.080 1.061 1.128 0.837 1.102 0.954 0.629 0.937 1.111 1.087

1.72e-02 2.30e-02 1.40e-02 1.58e-02 1.47e-02 3.04e-05 4.08e-04 2.98e-02 5.70e-02 2.95e-02 7.05e-03 3.72e-03 2.12e-02 2.99e-02 4.02e-07 1.40e-01 1.17e-02 7.96e-02 2.86e-01 9.01e-02 5.44e-03 1.88e-02

0.0293 0.0293 0.0295 0.0295 0.0297 0.0298 0.0386 0.0390 0.0395 0.0399 0.0402 0.0430 0.0566 0.0568 0.0576 0.0576 0.0576 0.0568 0.0635 0.0578 0.0539 0.0578

*

51

Representaci´on gr´afica Notemos que hay algunas observaciones que perturban los datos, como los son las observaciones 1, 2 y 32. Tomando estas observaciones respectivamente. Realizemos un an´alisis para cada una de estas observaciones, para notar que sucede en las estimaciones de los par´ametros del modelo y notar que modelo se ajusta mejor. El estadistico que se usar´a es: Z=

θˆ − θˆ(Ai ) θˆ

Donde el mayor valor de Z es la observaci´on que evidencia mas perturbaci´on y debe ser eliminado. Excluyendo la observaci´ on 1. Notemos que el nuevo modelo estimado viene dado por:

52

Yb = −0,638 + 1,142X1 > datos1 mod datos2 mod datos3 mod datos4 mod datos5 mod datos6 mod datos7 mod anova(mod) Analysis of Variance Table Response: Dep.Y Df Sum Sq Mean Sq F value Pr(>F) Indep.X 1 385.57 385.57 236.69 2.433e-16 *** Residuals 32 52.13 1.63 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 En donde se puede concluir que β1 es altamente significativa para nuevo modelo. Ahora obtengamos algunas estad´ısticas muy importantes para el vector de par´ametros β lo son los respectivos errores est´andar para cada par´ametro del modelo y los valores t para probar la hip´otesis H0 : βi = 0,adem´as obtenemos la estimaci´on de la desviaci´on est´andar del error σ > summary(mod) Call: lm(formula = Dep.Y ~ Indep.X, data = datos1) Residuals: Min 1Q Median -3.2889 -0.5844 -0.0772

3Q 1.0153

Max 2.1650

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.80913 0.59917 -1.35 0.186 Indep.X 1.18376 0.07694 15.38 2.43e-16 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.276 on 32 degrees of freedom

57 Multiple R-squared: 0.8809,Adjusted R-squared: 0.8772 F-statistic: 236.7 on 1 and 32 DF, p-value: 2.433e-16 Note que estas variables explicativas logran explicar al modelo en un 88.09 %, notemos que el nuevo modelos sin la observaci`on 32 logra ser un mejor el modelo ya que aumenta el ajuste del modelo, por ende es mejor quedarse con los datos iniciales eliminando la observaci`on 32 . Lo cual se verifica en R, por medio de > outlierTest(mod) No Studentized residuals with Bonferonni p < 0.05 Largest |rstudent|: rstudent unadjusted p-value Bonferonni p 32 -3.313649 0.0022947 0.080315 c) Realizamos los gr´aficos de i ∼ Yˆi y i ∼ Xi

(a)i ∼ Yˆi .

(b) i ∼ Xi

El gr´afico residual de residuos contra X independiente muestra que la varianza de error tiende a aumentar con X. Los residuos

58 tienden a situarse en una banda que diverge a medida que se mueve a lo largo del eje x. En general, si la banda dentro de la cual los residuos residen diverge (se hace m´as ancha) a medida que X aumenta, la varianza de error tambi´en aumenta con X. Por otra parte, si la banda converge (se hace m´as estrecha), la varianza del error disminuye con x. La gr´afica del residuo con respecto a la variable independiente {plot (X, resid (mod))} que acaba de trazar aqu´ı - se˜ nala la presencia de errores heterosced´asticos ya que los residuos aumentan con x. 13. Demostraremos que si centramos las variables en el modelo de regresi´on lineal entonces, 1 H = Jn + H c n 0 −1 0 en efecto, note que H = (X X) X es una matriz sim´etrica e ignipotente. Ahora definamos el modelo centrado, donde yi puede escribirse en terminos de variables x centradas como yi = β0 + β1 Xi1 + β2 Xi2 + · · · + βk Xik + i = α + β1 (Xi1 − X 1 ) + β2 (Xi2 − X 2 ) + · · · + βk (Xik − X k ) + i con i = 1, 2, · · · , n donde α = β0 + β1 X 1 + · · · + βk X k en forma matricial el modelo centrado para y1 , y2 , · · · , yk se convierte en   α y = (J, Xx ) + β1 donde   1 0 β1 = (β1 , β2 , · · · , βk ) y Xc = I − J X1 n y ademas   X11 X12 · · · X1k  X21 X22 · · · X2k    X1 =  .. .. .. ..   . . . .  Xn1 Xn2 · · · Xnk

59 la matriz I − ( n1 )J la llamaremos la matriz centrada . Ahora bien para nuestro modelos tomemos 0 0 α ˆ = y y βˆ1 = (Xc Xc )−1 Xc y

note que

0

0

yˆ = Xβˆ = X(X X)−1 X y = Hy

(1)

Ahora si definimos, 0

0

Hc = Xc (Xc Xc )−1 Xc )−1 Xc Ahora nuestro modelo estimado viene dado por 0

0

yˆ = yJ + Xc (X Xc )−1 X y 1 0 = ( J y) + Hc y n 1 = ( J + Hc )y n asi

1 yˆ = ( J + Hc )y n Ahora de ?? y ?? tenemos que H=

(2)

1 1 0 0 J + Hc = J + Xc (Xc Xc )−1 Xx n n

lo que finaliza la prueba 17. primero se realiza el an´alisis para Ejercicio 7 de regresi´on lineal multiple An´ alisis de Diagnostico El modelo que se esta estudiando es: Y = β0 + β1 X1 + β2 X2 + εi El modelo ajustado por MCO es:

60

Yb = 30,000 + 1,625X1 + 2,375X2 Ahora se hace un an´alisis de diagnostico para estudiar la presencia de observaciones o conjuntos de observaciones sobre la estimaci´on de los par´ametros y supuestos del modelo. En los modelos de regresi´on existe algunas medidas estad´ısticas propuestas para identificar y medir conjuntos de observaciones influyentes, las cuales son: Distancia de Cook, DFFITS, DFBETAS y COVRATIO. Para obtener las medidas de influencias antes mencionadas (Usando R). Para el ejemplo que llevamos, tenemos que: > influence.measures(mod) Influence measures of lm(formula = Y ~ X1 + X2) : dfb.1_ dfb.X1 dfb.X2 dffit cov.r cook.d hat inf 1 -1.71e-02 -0.00828 0.0304 0.0513 2.439 0.00109 0.205 2 5.40e-01 -0.66398 -0.1747 0.8452 1.261 0.22309 0.348 3 -3.00e+00 2.16025 2.1602 -3.3466 0.153 1.50538 0.571 * 4 1.94e-01 -0.15067 -0.1507 -0.2751 2.380 0.03011 0.286 5 8.60e-02 -0.30623 0.0835 -0.5171 1.233 0.08852 0.205 6 2.33e-01 -0.15088 -0.1509 0.5843 0.659 0.09409 0.143 7 1.76e-01 0.01904 -0.2475 -0.2856 3.991 0.03340 0.536 * 8 4.48e-16 1.02091 -0.6942 1.4057 3.776 0.68253 0.705 * Notemos que hay muchas observaciones que perturban los datos, como los son las observaciones 3, 7 y 8. Tomando estas observaciones respectivamente. Realizamos un an´alisis para cada una de estas observaciones, para notar que sucede en las estimaciones de los par´ametros del modelo y notar que modelo se ajusta mejor. El estad´ıstico que se usar´a es:

61

Z=

θˆ − θˆ(Ai ) θˆ

Donde el mayor valor de Z es la observaci´on que evidencia mas perturbaci´on y debe ser eliminado. Excluyendo la observaci´ on 3. Notemos que el nuevo modelo estimado viene dado por: Yb = 45,75 + 1,00X1 + 1,75X2 > datos1 modelo1 datos2 modelo2 datos3 modelo3 datos4 modelo4 datos5 modelo5 datos6 modelo6 datos7 modelo7 anova(modelo1) Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) X1 1 0.5143 0.5143 0.8229 0.415652 X2 1 14.7000 14.7000 23.5200 0.008341 ** Residuals 4 2.5000 0.6250 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

66 En donde se puede concluir que β2 es altamente significativa para nuevo modelo. Ahora obtengamos algunas estad´ısticas muy importantes para el vector de par´ametros β lo son los respectivos errores est´andar para cada par´ametro del modelo y los valores t para probar la hip´otesis H0 : βi = 0,adem´as obtenemos la estimaci´on de la desviaci´on est´andar del error σ > summary(modelo1) Call: lm(formula = Y ~ ., data = datos1) Residuals: 1 2 -1.021e-14 -2.500e-01

4 5 2.680e-15 -1.250e+00

6 7.500e-01

7 2.500e-01

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 45.7500 7.5519 6.058 0.00375 ** X1 1.0000 0.3608 2.771 0.05026 . X2 1.7500 0.3608 4.850 0.00834 ** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7906 on 4 degrees of freedom Multiple R-squared: 0.8589,Adjusted R-squared: 0.7883 F-statistic: 12.17 on 2 and 4 DF, p-value: 0.01992 Note que estas variables explicativas logran explicar al modelo en un 85.89 %, notemos que el nuevo modelos sin la observacion 3 logra ser un buen modelo a pesar de que disminuye el ajuste de este. No obstante, es mejor quedarse con el modelo inicial sin eliminar alguna observaci´on Lo cual se verifica en R, por medio de

8 5.000e-01

67 > library(car) > outlier.test(mod) No Studentized residuals with Bonferonni p < 0.05 Largest |rstudent|: rstudent unadjusted p-value Bonferonni p 3 -2.898275 0.044194 0.35355 Warning message: ’outlier.test’ is deprecated. Use ’outlierTest’ instead. See help("Deprecated") and help("car-deprecated"). Validaci´ on de Supuestos del Modelo • Linealidad Para el ajuste obtenido con los datos

Yb = 30,000 + 1,625X1 + 2,375X2 construiremos los gr´aficos de regresi´on parcial y extraeremos conclusiones sobre la relaci´on lineal entre regresores y respuesta. As´ı:

Efectivamente, la tendencia que se aprecia que tanto en la variable X1 como en X2 es lineal, por lo que no tenemos indicios de problemas con la hip´otesis de linealidad.

68 • Normalidad Este supuesto es muy importante puesto que es u ´til para probar hip´otesis, calcular bandas de confianza para los pron´osticos de nuestro modelo ajustado debido a que es el objetivo principal de los modelos de regresi´on. Las hip´otesis a contrastar son las siguientes: H0 = Los datos provienen de una distribuci´on normal H1 = Los datos no provienen de una distribuci´on normal veamos graficamente como se comporta los residuales para el modelo

Prueba de Shapiro- Wilk > shapiro.test(resid(mod)) Shapiro-Wilk normality test data: resid(mod) W = 0.97229, p-value = 0.9154 Como tenemos un p-valor = 0.9154 > 0.05 no se rechaza la hip´otesis de normalidad, esto es, los residuales del modelo siguen una distribuci´on normal. • Independencia

69 Para realizar contrastes de autocorrelaci´on entre los residuales hay que especificar la hip´otesis alternativa que defina un esquema de autocorrelaci´on. La hip´otesis a probar es: H0 : Cov(i ; 0i ) = 0 para todo i 6= i0 H1 : Cov(i ; 0i ) 6= 0 para alg´ un i 6= i0 Realizando el test de independencia para los residuales del modelo ajustado, se tiene: Durbin Watson > durbinWatsonTest(mod) lag Autocorrelation 1 -0.3165323 D-W Statistic p-value 2.580645 0.426 Alternative hypothesis: rho != 0 Como el p-valor asociado a la prueba de independencia es mayor que 0.05 no rechazamos la hip´otesis H0 : Cov(i ; 0i ) = 0, por tanto se concluye que existe independencia entre los residuales del modelo. Rachas > X=as.factor(sign(resid(mod))) > runs.test(X) Runs Test data: X Standard Normal = 0, p-value = 1 alternative hypothesis: two.sided Como el p-valor asociado al test de independencia mediante la pruebas de rachas resulta ser superior a 0.05 concluimos que existe independencia en los errores.

70 • Homogeneidad de Varianza Primero miremos el gr´afico de dispersi´on de los residuales del modelo

En el gr´afico hay indicios de que la varianza de los errores tiene un patr´on a crecer y decrecer, para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de Breusch-Pagan, teniendo como resultado. La hip´otesis nula en ´esta prueba es:

H0 : Los residuales tienen varianza constante H1 : Los residuales no tienen varianza constante > bptest(modelo1) studentized Breusch-Pagan test data: modelo1 BP = 1.0587, df = 2, p-value = 0.589 De acuerdo con los resultados de la prueba, el p-valor es 0.589 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´otesis nula, es

71 decir que el modelo no presenta problemas de homocedasticidad. Por otro lado, el analisis diagnostico encontro que las observaciones 3, 7 y 8 ejercian infuencia en las estimativas de los parametros, entonces surge la inquietud si el problema de falta de homogeneidad entre las varianzas que se presentan en los residuos pueden ser ocasionados por la alta influencia que ejercen estas observaciones Realizamos ahora el analisis del modelo eliminando la observacio´on 3. Se ajustara un nuevo modelo eliminando dichas observaciones. A continuaci´on se muestran las estimativas de los nuevos par´ametros, R2 y pruebas de validacion entre otras. Validaci´ on de los Supuestos sin la obs. 3 • Linealidad Para el ajuste obtenido con los datos, sin incluir la observaci´on 17, se tiene: Yb = 45,75 + 1,00X1 + 1,75X2 construiremos los gr´aficos de regresi´on parcial y extraeremos conclusiones sobre la relaci´on lineal entre regresores y respuesta. As´ı:

72 Not´ese que, la tendencia que se aprecia para la variable X1 es evidentemente lineal, igualmente que para X2 , por lo que no se ve afectada la linealidad de los datos al momento de eliminar la observaci´on 17. • Normalidad veamos graficamente como se comporta los residuales para el nuevo modelo:

Para verificar dicha normalidad usamos el test de ShapiroWilk (n shapiro.test(resid(modelo1)) Shapiro-Wilk normality test data: resid(modelo1) W = 0.91054, p-value = 0.3996 como se tiene que el p-valor= 0.3996> 0.05 se concluye que no se rechaza la hipotesis de normalidad,es decir; los errores del modelo siguen una distribucio normal.

73 • Independencia • Prueba de rachas > X=as.factor(sign(resid(modelo1))) > runs.test(X) Runs Test data: X Standard Normal = -0.3638, p-value = 0.716 alternative hypothesis: two.sided Como el p-valor= 0.716> 0.05, se concluye que existe independencia entre los errores. • Durbin-Watson > dwtest(modelo1) Durbin-Watson test data: modelo1 DW = 2.4, p-value = 0.8025 alternative hypothesis: true autocorrelation is greater than 0 Como p-valor= 0.8025> 0.05 , se rechaza la hipotesis nula. Por lo tanto se concluye que hay independencia entre los residuales del modelo. • Homogeneidad de varianza Grafico de dispersion de los residuales del modelo.

74

Para probar si existe homogeneidad de varianzas para los errores del modelo de regresion ajustado se usa test de BreuschPagan. > bptest(modelo1,varformula= NULL,studentize=F) Breusch-Pagan test data: modelo1 BP = 1.1473, df = 2, p-value = 0.5635 Para este caso la hipotesis nula es la varianza constante, el p-valor= 0.5635 de la prueba comparado con el nivel de significancia α=0,05 es mayor, por lo que no se permite rechazar homocedasticidad. Luego el problema de heterocedasticidad en el ajuste inicial queda corregido. • ahora se realiza el an´alisis para el ejerci´o 10 de regresi´on lineal m´ ultiple Y X1 X2 Y X1 X2

12.1 5.5 4.6 4.5 10.8 4.9 0.870 0.202 0.203 0.198 0.730 0.510 1.69 1.17 1.17 1.21 1.63 1.59 6.0 4.2 5.3 6.7 4.0 6.1 0.205 0.670 0.205 0.271 0.203 0.264 1.14 1.92 1.22 1.71 1.16 1.37

75 ´ ANALISIS DE DIAGNOSTICO El modelo que se esta estudiando es: Y = β0 + β1 X1 + β2 X2 + εi El modelo ajustado por MCO es: Yb = 7,888 + 11,074X1 − 4,130X2 Ahora se hace un an´alisis de diagnostico para estudiar la presencia de observaciones o conjuntos de observaciones sobre la estimaci´on de los par´ametros y supuestos del modelo. En los modelos de regresi´on existe algunas medidas estad´ısticas propuestas para identificar y medir conjuntos de observaciones influyentes, las cuales son: Distancia de Cook, DFFITS, DFBETAS y COVRATIO. Para obtener las medidas de influencias antes mencionadas (Usando R). Para el ejemplo que llevamos, tenemos que:

> influence.measures(mod) Influence measures of lm(formula = Y ~ X1 + X2) : dfb.1_ 1 0.3588 2 0.0303 3 -0.1043 4 -0.0580 5 0.1424 6 0.1236 7 0.1022 8 1.8884 9 0.0169 10 -3.9993 11 -0.2210 12 -0.0150

dfb.X1 1.096563 0.000285 -0.001839 0.022142 0.478202 -0.016241 0.017770 0.586945 -0.006832 -4.201851 -0.015932 -0.090309

dfb.X2 -0.54979 -0.02007 0.06935 0.02822 -0.21605 -0.13445 -0.07516 -1.79618 -0.00793 4.61593 0.15243 0.05632

dffit 1.3367 0.0488 -0.1670 -0.1274 0.6523 -0.4471 0.1425 -2.5712 0.0384 5.0497 -0.3369 0.1840

cov.r 1.769 1.677 1.593 1.590 1.430 0.977 1.669 0.229 1.636 0.128 1.372 1.455

cook.d 0.559642 0.000892 0.010268 0.006011 0.141915 0.063366 0.007521 1.130521 0.000552 3.136869 0.039631 0.012292

hat inf 0.531 * 0.155 0.155 0.138 0.299 0.120 0.176 0.409 * 0.133 0.609 * 0.161 0.113

76

Notemos que las observaciones 1, 8 y 10, perturban los datos, de tal manera que vamos a realizar un an´alisis a fondo, para verificar si realmente son influyentes o son efectos naturales de los datos. Realizemos un an´alisis para cada una de estas observaciones, para notar que sucede en las estimaciones de los par´ametros del modelo y notar que modelo se ajusta mejor. El estadistico que se usar´a es: Z=

θˆ − θˆ(Ai ) θˆ

Donde el mayor valor de Z es la observaci´on que evidencia mas perturbaci´on y debe ser eliminado. Excluyendo la observaci´ on 1.

77 Notemos que el nuevo modelo estimado viene dado por: Yb = 6,559 + 7,037X1 − 2,309X2 > modelo1 datos2 modelo2 datos3 modelo3 datos4 modelo4 datos5 modelo5 datos6 modelo6 datos7 modelo7 anova(modelo1) Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) X1 1 8.3123 8.3123 2.5344 0.1501 X2 1 1.3384 1.3384 0.4081 0.5408 Residuals 8 26.2384 3.2798 En donde se puede concluir que ni β1 , ni β2 son significativos para el nuevo modelo. Ahora obtengamos algunas estad´ısticas muy importantes para el vector de par´ametros β lo son los respectivos errores est´andar para cada par´ametro del modelo y los valores t para probar la hip´otesis H0 : βi = 0,adem´as obtenemos la estimaci´on de la desviaci´on est´andar del error σ > summary(modelo1) Call: lm(formula = Y ~ ., data = datos1)

82 Residuals: Min 1Q -2.6395 -0.9970

Median 0.1159

3Q 0.7391

Max 2.8686

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.559 3.852 1.703 0.127 X1 7.037 4.887 1.440 0.188 X2 -2.309 3.615 -0.639 0.541 Residual standard error: 1.811 on 8 degrees of freedom Multiple R-squared: 0.2689,Adjusted R-squared: 0.08613 F-statistic: 1.471 on 2 and 8 DF, p-value: 0.2857 Note que estas variables explicativas logran explicar al modelo en un 26.89 %, notemos que el nuevo modelos sin la observacion 1 logra ser un mal modelo pues el ajuste del mismo disminuye notoriamente. Por lo tanto es mejor quedarse con el modelo inicial sin eliminar la observaci´on Validaci´ on De Supuestos Del Modelo • Linealidad Para el ajuste obtenido con los datos

Yb = 7,888 + 11,074X1 − 4,130X2 construiremos los gr´aficos de regresi´on parcial y extraeremos conclusiones sobre la relaci´on lineal entre regresores y respuesta. As´ı:

83

Efectivamente, la tendencia que se aprecia que tanto en la variable X1 como en X2 es lineal, a pesar de que en X2 esta linealidad se da de manera inversa, por lo que no tenemos indicios de problemas de linealidad. • Normalidad Este supuesto es muy importante puesto que es u ´til para probar hip´otesis, calcular bandas de confianza para los pron´osticos de nuestro modelo ajustado debido a que es el objetivo principal de los modelos de regresi´on. Las hip´otesis a contrastar son las siguientes: H0 = Los datos provienen de una distribuci´on normal H1 = Los datos no provienen de una distribuci´on normal veamos gr´aficamente como se comporta los residuales para el modelo

84

Prueba de Shapiro- Wilk > shapiro.test(resid(mod)) Shapiro-Wilk normality test data: resid(mod) W = 0.98694, p-value = 0.9985 Como tenemos un p-valor = 0.9985 > 0.05 no se rechaza la hip´otesis de normalidad, esto es, los residuales del modelo siguen una distribuci´on normal. • Independencia Para realizar contrastes de autocorrelaci´on entre los residuales hay que especificar la hip´otesis alternativa que defina un esquema de autocorrelaci´on. La hip´otesis a probar es: H0 : Cov(i ; 0i ) = 0 para todo i 6= i0 H1 : Cov(i ; 0i ) 6= 0 para alg´ un i 6= i0 Realizando el test de independencia para los residuales del modelo ajustado, se tiene: Durbin Watson

85 > durbinWatsonTest(mod) lag Autocorrelation D-W Statistic p-value 1 -0.3702329 2.634657 0.182 Alternative hypothesis: rho != 0 Como el p-valor asociado a la prueba de independencia es mayor que 0.05 no rechazamos la hip´otesis H0 : Cov(i ; 0i ) = 0, por tanto se concluye que existe independencia entre los residuales del modelo. Rachas > X=as.factor(sign(resid(mod))) > runs.test(X) Runs Test data: X Standard Normal = 1.3533, p-value = 0.1759 alternative hypothesis: two.sided Como el p-valor asociado al test de independencia mediante la pruebas de rachas resulta ser superior a 0.05 concluimos que existe independencia en los errores. • Homogeneidad De Varianza Primero miremos el gr´afico de dispersi´on de los residuales del modelo En el gr´afico hay indicios de que la varianza de los errores tiene un patr´on a crecer y decrecer, para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de Breusch-Pagan, teniendo como resultado. La hip´otesis nula en ´esta prueba es:

H0 : Los residuales tienen varianza constante H1 : Los residuales no tienen varianza constante

86

Figura 1: dispersi´on de los residuales > bptest(mod,varformula= NULL,studentize=F) Breusch-Pagan test data: mod BP = 8.0274, df = 2, p-value = 0.01807 De acuerdo con los resultados de la prueba, el p-valor es 0.589 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´otesis nula, es decir que el modelo no presenta problemas de homocedasticidad. Apoyandonos en el criterio, se decide que no es necesario eliminar ninguna observaci´on dado que los datos ajustan un buen modelo y al realizar otro, el ajuste en vez de aumentar, disminuye significativamente.

87

0.8

1.0

1.2

1.4

100

150

0.6

1.2

1.4

50

Y

3.0

3.5

4.0

0.6

0.8

1.0

X1

1.0

1.5

2.0

2.5

X2

50

100

150

1.0

1.5

2.0

2.5

3.0

3.5

4.0

Figura 2: Dispersi´on entre las variables

Modelos Polin´ omicos 1. Notamos un comportamiento lineal entre X2 e Y mientras que es posible que entre las variables X1 e Y no exista una relaci´on de tipo lineal ya que en su dispersi´on no se ve ese comportamiento figura ??. Ajustamos el modelo de regresi´on polinomica: 2 2 Yˆi = βˆ0 + βˆ1 Xi1 + βˆ2 Xi2 + βˆ12 Xi1 Xi2 + βˆ11 Xi1 + βˆ22 Xi2

ahora implementando R tenemos > mod=lm(Y~X1+X2+I(X1*X2)+I(X1^2)+I(X2^2));mod Call: lm(formula =Y~ X1+ X2 + I(X1*X2) + I(X1^2) + I(X2^2)) Coefficients: (Intercept) X1 24.410 -38.033

X2 0.72

I(X1 * X2) -9.986

I(X1^2) 34.975

I(X2^2) 11.066

88 2 2 Yˆi = 24,410−38,033Xi1 +0,72Xi2 −9,986Xi1 Xi2 +34,975Xi1 +11,066Xi2

luego un resumen de este modelo ajustado es: > res=summary(mod);res Call: lm(formula =Y~ X1+ X2 + I(X1*X2) + I(X1^2) + I(X2^2)) Residuals: Min 1Q Median -6.3504 -2.7354 -0.3533

3Q 2.7022

Max 8.9329

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 24.410 26.592 0.918 0.3941 X1 -38.033 40.452 -0.940 0.3834 X2 0.720 11.687 0.062 0.9529 I(X1 * X2) -9.986 8.742 -1.142 0.2969 I(X1^2) 34.975 21.556 1.623 0.1558 I(X2^2) 11.066 3.158 3.504 0.0128 * --Signif.codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 Residual standard error: 6.042 on 6 degrees of freedom Multiple R-squared: 0.9938,Adjusted R-squared: 0.9886 F-statistic: 192.2 on 5 and 6 DF, p-value: 1.556e-06 Para el ajuste a pesar de que los par´ametros no son significativos, ya que el p − valor asociado a la prueba H0 : βj = 0 son mayor que 0.05 para casi todos los par´ametros del modelo. Muy a pesar que todos los coeficientes del modelo no resultan significativos, el coeficiente de determinaci´on es muy bueno (R2 = 0,9938 %), y un (Ra2 = 0,9886 %) dado que la variable independientes explican un (98,86 % la variabilidad total de la variable dependiente Y ; asimismo tambi´en se obtiene un AIC = 82,90815 Ahora se realiza un an´alisis diagn´osticos para analizar si existen observaciones at´ıpicas o muy influyentes sobre las estimativas de los par´ame-

89 tros del modelo, las cuales tambi´en pueden influenciar en las pruebas de validaci´on de los errores del modelo. An´ alisis diagnostico para el modelo: Un resumen del an´alisis de diagnostico graficamente y en R es: > summary(influence.measures(mod)) Potentially influential observations of lm(formula =Y~ X1+ X2 + I(X1*X2) + I(X1^2) + I(X2^2)):

2 7 8 10 11 2 7 8 10 11

dfb.1_ dfb.X1 dfb.X2 dfb.I(*X dfb.I(X1 dfb.I(X2 0.03 0.03 -0.09 -0.05 -0.01 0.09 -0.01 -1.28_* 1.49_* -0.86 2.08_* -0.94 0.10 -0.08 -0.04 0.41 -0.08 -0.27 -0.66 0.58 0.29 -0.53 -0.26 0.18 -1.35_* 1.23_* 1.24_* 0.57 -1.41_* -1.14_* dffit cov.r cook.d hat 0.16 5.23_* 0.00 0.45 3.90_* 329.27_* 2.97_* 0.99 -0.56 12.30_* 0.06 0.78 -0.89 4.74_* 0.15 0.63 1.75 5.26_* 0.52 0.78

Eliminaci´ on de observaciones De los gr´aficos de las figuras ?? y ?? se concluye que las observaciones que est´an disminuyendo la precisi´on en la estimaci´on de los par´ametros son la #7 y la #11, por lo tanto ajustaremos nuevamente el modelo cuadr´atico eliminando estas observaciones y se realizaran las respectivas pruebas de validaci´on de los residuos. Los resultados arrojados en el paquete R son: > > > >

Y1=Y[-c(7,11)]; X11=X1[-c(7,11)]; X12=X2[-c(7,11)]; mod1=lm(Y1~X11+X12+I(X11*X12)+I(X11^2)+I(X12^2));mod1

Call: lm(formula=Y1~X11+X12+I(X11*X12)+I(X11^2)+I(X12^2))

DFBETAS para X1^2

0

2

2 4

4 6

Index

6 8

8

1.5

10

10

10 12

12

12

DFBETAS para X1X2

1.0

8

0.5

DFBETAS para X2

1.0

6

hatvalues(mod)

−1

1

DFFITS 2

3

Index 2

2

2 4

Index 4

4 6

6

6 8

8

8 10

10

10

−0.5

0.0

4

6

−1.0

DFBETAS para el intercepto

0

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

4

−0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6

0.0

0.5

4

Çovratio

0.0

0.0

DFBETAS para X1 −0.5

2

−0.5

DFBETAS para X2^2

−1.0

2

50 100 150 200 250 300

−1.0

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0

0.0

0.5

1.0

1.5

Cook's distance 2.0

2.5

3.0

90

Cook's distance

7

3 11

Obs. number lm(Y ~ X1 + X2 + I(X1 * X2) + I(X1^2) + I(X2^2)) 8 10

12

12

Index

12

Figura 4: Gr´aficos de diagn´ostico 12

Figura 3: Gr´afico de la distancia de Cook

Index 2

2

2 4 Index

Index 4

4 6

Index 6

Index

6 8 10 12

8 10 12

8

10

12

91 Coefficients: (Intercept) X11 131.7 -498.9

X12 I(X11*X12 62.4 -70.4

I(X11^2) I(X12^2) 396.2 4.0

> res=summary(mod1);res Call: lm(formula = Y1~X11+X12+I(X11*X12)+I(X11^2)+I(X12^2)) Residuals: 1 2 3 4 5 1.000e+00 -1.000e+00 9.000e+00 -6.000e+00 -3.000e+00 6 7 8 9 10 3.300e+00 1.100e+00 -6.531e-14 -1.100e+00 -3.300e+00 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 131.70 88.40 1.490 0.211 X11 -498.90 435.39 -1.146 0.316 X12 62.40 82.34 0.758 0.491 I(X11 * X12) -70.40 58.91 -1.195 0.298 I(X11^2) 396.20 348.35 1.137 0.319 I(X12^2) 4.00 12.34 0.324 0.762 Residual standard error: 6.168 on 4 degrees of freedom Multiple R-squared: 0.9956,Adjusted R-squared: 0.9902 F-statistic: 182.1 on 5 and 4 DF, p-value: 8.333e-05 El ajuste del modelo aumenta y a pesar que todos los coeficientes del modelo no resultan significativos, el coeficiente de determinaci´on es muy bueno (R2 = 0,9956 %), y un (Ra2 = 0,9902 %) dado que la variable independientes explican un (99,02 % la variabilidad total de la variable dependiente Y ; asimismo tambi´en se obtiene un AIC = 69,60487 Ahora realizamos los supuestos del modelo en R

92 Normalidad: > library(car) > e1=resid(mod1) > shapiro.test(e1) Shapiro-Wilk normality test data: e1 W = 0.93884, p-value = 0.5401 Homogeneidad de Varianza : > library(lmtest) > bptest(mod1) studentized Breusch-Pagan test data: mod1 BP = 5.4874, df = 5, p-value = 0.3593 Independencia : > durbin.watson(mod1) lag Autocorrelation D-W Statistic p-value 1 -0.3195795 2.561038 0.616 Alternative hypothesis: rho != 0 Notamos que despu´es de eliminadas las observaciones #7 y #11, el modelo cumple con todos los supuestos y adem´as el coeficiente de determinaci´on ajustado aumento al 99,02 % y hubo una disminuci´on considerable en el AIC, el cual paso de 82.91 a 69.61, por lo cual se concluye que el modelo ajusta mucho mejor que el modelo con todas las observaciones. 3.

a) Calcule PRESS para el ajuste de la regresi´on de m´ınimos cuadrados para el modelo anterior Tenemos que P RESS =

n X i=1

ri 1 − hii

(modelo general)

93 Primero ajustando el modelo para la eficiencia de la transferencia de calor (Y1 ), donde obtenemos > mod library(DAAG) > press(mod) [1] 941921.7

ri = 941921 1 − hii

94 b) Ajuste un modelo de segundo orden con X4 eliminado en forma adsoluta (es decir, eliminando todos los t´erminos con X4 ) calcule el criterio de de predicci´on para el modelo reducido.Comente lo apropiado de X4 para pronosticar el coeficiente de transferencia de calor. Un modelo ajustado, sin la variable X4 , viene dado por: > mod2 mod library(DAAG) > press(mod3) [1] 111465.8 Ahora, un modelo ajustado, sin la variable X4 , viene dado por: > mod4F) X 1 42.588 42.588 477103 < 2.2e-16 *** I(X^2) 1 4.637 4.637 51946 8.264e-15 *** Residuals 7 0.001 0.000 --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 scrm=42.588+4.637 scrm [1] 47.225 Ahora SCRm = 42,588 + 4,637 = 47,225 SCE = 0,001, luego m = 23,6125, CM E = SCE = 0,001 = 0,0001428571, CM R = SCR 2 gle 7 FC = 165287,5, F5 % = 4,737414 dado que el valor del estad´ıstico FC > F5 % rechazamos la hip´otesis nula. Concluyendo que por lo menos uno entre β1 y β2 es distinto de cero, esto es el modelo explica una porci´on significante de la variable respuesta Y. c) Para probar la Hip´otesis H0 : β1 = 0, usaremos la funci´on summary, as´ı summary(mod) Call: lm(formula = Y ~ X + I(X^2)) Residuals: Min 1Q Median 3Q Max -0.0164545 -0.0020606 0.0001364 0.0066136 0.0094545 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.641000 0.011112 147.7 1.72e-13 *** X -1.249394 0.018564 -67.3 4.20e-11 *** I(X^2) 1.499394 0.006579 227.9 8.26e-15 *** --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

98

Residual standard error: 0.009448 on 7 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 2.645e+05 on 2 and 7 DF, p-value: < 2.2e-16 Ahora bien bajo la hip´otesis H0 : β11 = 0, y dado que el valor calculado de la estad´ıstica Fc = 227,9 > 8,26e − 15 es mucho mayor que el tabulado a un nivel de significancia de α = 0 : 05; rechazando (H0 ), se concluye que el par´aametro β11 es significativo para el modelo. d) Analizar el ajuste del modelo, calcule el PRESS; AIC;BIC summary(mod) Call: lm(formula = Y ~ X + I(X^2)) Residuals: Min 1Q Median 3Q Max -0.0164545 -0.0020606 0.0001364 0.0066136 0.0094545 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.641000 0.011112 147.7 1.72e-13 *** X -1.249394 0.018564 -67.3 4.20e-11 *** I(X^2) 1.499394 0.006579 227.9 8.26e-15 *** --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Residual standard error: 0.009448 on 7 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 2.645e+05 on 2 and 7 DF, p-value: < 2.2e-16 AIC(mod) [1] -60.42709 BIC(mod) [1] -59.21675 require(DAAG) press(mod) [1] 0.001051006

99 Notemos que para el ajuste del modelo polinomico de segundo grado se tienen par´ametros altamente significativos, ya que el p-valor asociado a la prueba H0 : βj = 0 es menor que 0.05 para todos los par´ametros del modelo. Ahora bien el coeficiente de determinaci´on es muy bueno (R2 = 100 %), esto es que la variable independiente X explica en un 100 % la variabilidad total de la variable respuesta Y. Ademas note que AIC = −60,42709, BIC = −59,21675, P RESS = 0,001051006 AIC(mod) [1] -60.42709 BIC(mod) [1] -59.21675 require(DAAG) press(mod) [1] 0.001051006 e) Ajustaremos el modelo cuadr´atico Yˆ = β0 + β1 (Xi − X) + β11 (Xi − X)2 de la siguiente manera xb=mean(X) Xi=X-xb mod1=lm(Y~Xi+I(Xi^2)) mod1 Call: lm(formula = Y ~ Xi + I(Xi^2)) Coefficients: (Intercept) Xi I(Xi^2) 2.758 2.874 1.499 f) calculamos en R en AIC,BIC,PRESS del modelo AIC(mod1) [1] -60.42709 BIC(mod1) [1] -59.21675

100 press(mod1) [1] 0.001051006 Notamos que los AIC.BIC,PRESS son iguales en ambos modelos

Modelo con Variables Cualitativas 1.

a) Definamos Y = crecimiento de baterias, X2 = tiempo,  1 si la informaci´on proviene de la bateria B X1 = 0 si la informaci´on proviene de la bateria A Se ajusta el modelo lineal Yi = β0 + βi1 Xi1 + β2 Xi2 + β3 Xi1 Xi2 + εi en R, definiendo previamente las siguientes variables tiempo=c(-2, -1, 0, 1, 2) A=c(8.0, 9.0, 9.1, 10.2, 10.4) B=c(10.0, 10.3, 12.2, 12.6, 13.9) Y=c(A,B) x2=c(tiempo,tiempo) x1=c(rep(1,5),rep(0,5)) tal que Call: lm(formula = Y ~ x1 + x2 + I(x1 * x2)) Coefficients: (Intercept) 11.80

x1 -2.46

x2 1.01

I(x1 * x2) -0.41

as´ı, el modelo ajustado es Yˆi = 11,80 − 2,46Xi1 + 1,01Xi2 − 0,41Xi1 Xi2 b) β0 es el intercepto del modelo para las bater´ıas A, β1 es la diferencia entre la bacteria B y la bacteria A, β2 es el intercepto para el modelo cuando las bater´ıas provienen de B y β3 es la diferencia para el intercepto entre las bacterias B y las bacterias A

101 c) A continuaci´on para calcular el modelo para cada tipo de bacteria, consideremos los siguientes casos si Xi1 = 0 tenemos Yˆi = 11,80 − 2,46Xi1 + 1,01Xi2 − 0,41Xi1 Xi2 = 11,80 − 2,46(0) + 1,01Xi2 − 0,41(0)Xi2 = 11,80 + 1,01Xi2

Xi1 = 1 Yˆi = = = =

11,80 − 2,46Xi1 + 1,01Xi2 − 0,41Xi1 Xi2 11,80 − 2,46(1) + 1,01Xi2 − 0,41(1)Xi2 11,80 − 2,46 + 1,01Xi2 − 0,41Xi2 9,34 + 0, 6Xi2

es decir,  11,80 + 1,01Xi2 si la informaci´on proviene de la bateria A ˆ Yi = 9,34 + 0, 6Xi2 si la informaci´on proviene de la bateria B Luego si se toma Xi1 = 0 es decir la informaci´on proviene de la bateria A, tenemos que las tasa de crecimiento de dicha bateria varia respecto al tiempo X2 a raz´on de 1.01; esperando que en un tiempo cero la tasa de crecimiento sea de 11.8; mientras se toma Xi1 = 1 es decir que la informaci´on proviene de la bateria B, tenemos que las tasa de crecimiento de esta bateria varia respecto al tiempo X2 a raz´on de 0,6; en este caso esperando que en un tiempo cero la tasa de crecimiento sea de 9.34. d ) De los anteriores resultados tenemos que existen diferencia en las tasas de crecimiento entre las bacterias de tipo A y B, esto dado que la pendiente como el intersepto en el modelo que describe el crecimiento de la bacteria A es mayor a la del modelo que describe el crecimiento de la bacteria B. e) A continuaci´on se verifica el cumplimiento de los supuestos para el modelo estimado. Prueba de normalidad de los errores

102 Shapiro-Wilk normality test data: modelo$residuals W = 0.93982, p-value = 0.551 Prueba de homogeneidad de varianza de los errores studentized Breusch-Pagan test data: modelo BP = 3.1603, df = 3, p-value = 0.3676 Prueba de aleatoriedad o independencia de los errores Runs Test data: X Standard Normal = 2.6833, p-value = 0.00729 alternative hypothesis: two.sided Como se puede observar, el u ´nico supuesto que presenta problemas (al 5 % de significancia) es el deindependencia de los errores, dado que el valor del p-value esta alejado del nivel de significancia se considera severo dicho problema. Ahora se realiza un an´alisis de diagnostico, para verificar la existencia de datos influyentes. En efecto, > library(faraway) > rs iflSR identical(rs, rstandard(modelo, infl = iflSR)) [1] TRUE Dado el resultado en el paso anterior, se procede a calcular las medidas de influencia para los datos, tal que > inflm.SR summary(inflm.SR) Potentially influential observations of lm(formula = Y ~ x1 + x2 + I(x1 * x2)) : dfb.1_ 1 0.00 5 0.00 7 -1.12_* 10 0.24

dfb.x1 dfb.x2 dfb.I(*x -0.30 0.00 0.42 -0.30 0.00 -0.42 0.79 0.79 -0.56 -0.17 0.33 -0.24

dffit cov.r cook.d hat -0.73 3.93_* 0.15 0.60 -0.73 3.93_* 0.15 0.60 -1.38 0.24 0.30 0.30 0.41 4.74_* 0.05 0.60

De donde los datos infleyentes son el 1, 5, 7, y 10, tal que segun la medida de DFBETA la observaci´on 7 influyen el la estimaci´on del β0 ; mientras que la medida de COVRATIO, las observaciones 1, 5 y 10 son las que resultan influyentes, es decir estas influyen en la estimaci´on de la matr´ız de covarianza de los estimadores de los coeficientes. 2. Definamos Y = log(T ) X2 = log10 (W BC)),  1 para el grupoAg + X1 = 0 para el grupoAg − Primero ajustaremos el modelo Yt = β0 + β1 Xt1 + β2 Xt2 + β3 Xt1 Xt2 + t En R tenemos

104 tiempo=c(Agmastiempo,agmentiempo) y=log(tiempo) wbc=c(agmasWBC,admenWBC) x2=log10(wbc) x1=c(replicate(16,1),replicate(16,0)) data=data.frame(y,x1,x2)

mod3=lm(y~x1+x2+I(x1*x2)) mod3 Call: lm(formula = y ~ x1 + x2 + I(x1 * x2)) Coefficients: (Intercept) 5.2529

x1 4.7445

x2 -0.7046

I(x1 * x2) -0.8866

veamos unas estadisticas importantes para el modelo summary(mod3) Call: lm(formula = y ~ x1 + x2 + I(x1 * x2)) Residuals: Min 1Q -2.04120 -0.85979

Median 0.05471

3Q 0.60911

Max 2.13318

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.2529 2.1719 2.419 0.0223 * x1 4.7445 2.9718 1.596 0.1216 x2 -0.7046 0.5076 -1.388 0.1760 I(x1 * x2) -0.8866 0.7105 -1.248 0.2224 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

105 Residual standard error: 1.154 on 28 degrees of freedom Multiple R-squared: 0.4448,Adjusted R-squared: 0.3853 F-statistic: 7.477 on 3 and 28 DF, p-value: 0.000795 Vemos que ninguno de los par´ametros resulto significativo para el modelo solo β0 y que el modelo esta ajustando un 38.53 %. Luego Los intervalos de confianza del 95 % para los par´ametros estimados son: confint(mod3) 2.5 % 97.5 % (Intercept) 0.8038308 9.7019019 x1 -1.3430294 10.8319903 x2 -1.7442678 0.3350969 I(x1 * x2) -2.3420603 0.5687737 veamos un resumen del analisis de varianza anova(mod3) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x1 1 13.672 13.6720 10.2629 0.003374 ** x2 1 14.137 14.1370 10.6120 0.002942 ** I(x1 * x2) 1 2.075 2.0745 1.5572 0.222405 Residuals 28 37.301 1.3322 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 verifiquemos los supuestos del modelo Normalidad Miremos el gr´afico de probabilidad normal para los residuales del modelo > shapiro.test(mod3$res) Shapiro-Wilk normality test data: mod3$res W = 0.97577, p-value = 0.6707

106

0 −2

−1

resid(mod3)

1

2

Normalidad resi

−2

−1

0

1

2

norm quantiles

Figura 5: Grafico de probabilidad normal para los residuales Para el juego de hip´otesis H0 : los residuales se distribuyen normal vs H1 : los residuales no se distribuyen normal a un nivel de significancia del 0.05 no se rechaza la hip´otesis nula ya que el p − valor = 0,6707 > 0,05 Por tanto los residuales del modelo siguen una distribuci´on normal Independencia Para verificar la independencia para los errores del modelo ajustado usamos el test de Durbin-Watson durbinWatsonTest(mod3) lag Autocorrelation D-W Statistic p-value 1 0.08239369 1.718582 0.214 Alternative hypothesis: rho != 0 La hip´otesis nula en ´esta prueba es son H0 : Los residuales son independientes vs H1 : Los residuales son independientes

107 El p-valor de la prueba es de 0,188 > 0,05. Se acepta la hip´otesis nula, es decir que se cumple el supuesto de que los errores son independientes.Luego para reafirmar la independencia de los residuales usaremos la prueba de rachas X=as.factor(sign(resid(mod3))) > runs.test(X) Runs Test data: X Standard Normal = -1.0782, p-value = 0.2809 alternative hypothesis: two.sided Como el p-valor asociado al test de independencia mediante la pruebas de rachas resulta ser superior a 0.05 concluimos que existe independencia en los errores. Homogeneidad de Varianzas Primero miremos el gr´afico de dispersi´on de los residuales del modelo En el gr´afico parece que la varianza de los errores no tiene ning´ un patr´on a crecer o decrecer, para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de Breusch- Pagan de la libreria lmtest , teniendo como resultado require(het.test) > data=data.frame(x1;x2; y ) Error: unexpected ’;’ in "data=data.frame(x1;" > modelo=VAR(data,p=2) > whites.htest(modelo) White’s Test for Heteroskedasticity: ==================================== No Cross Terms H0: Homoskedasticity H1: Heteroskedasticity

108

0 −2

−1

Residuales

1

2

Grafico de residuales

2

3

4

5

Ajustados

Figura 6: Gr´afico dispersi´on de los residuales

Test Statistic: 62.1805 Degrees of Freedom: 72 P-value: 0.7889 De acuerdo con los resultados de la prueba, el p-valor es 0.7889 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´otesis nula, es decir que el modelo no presenta problemas de homocedasticidad. Ahora hagamos un an´alisis de diagnostico para estudiar la presencia de observaciones o conjuntos de observaciones sobre la estimaci´on de los par´ametros y supuestos del modelo An´ alisis de Diagnostico

109 Realizemos un resumen del an´alisis, en R tenemos la siguiente salida summary(influence.measures(mod3)) Potentially influential observations of lm(formula = y ~ x1 + x2 + I(x1 * x2)) :

2 14 16

dfb.1_ dfb.x1 dfb.x2 dfb.I(*x dffit cov.r cook.d hat 0.00 -0.17 0.00 0.16 -0.26 1.65_* 0.02 0.31 0.00 0.63 0.00 -0.71 -1.19_* 0.80 0.31 0.23 0.00 -0.66 0.00 0.75 1.25_* 0.75 0.34 0.23

Donde vemos que hay algunas datos que tienen influencia sobre el modelo y los supuestos. Ahora para verificar hagamos las siguientes gr´aficas Notamos que la inclusi´on de las observaciones 14 y 16 disminuyen

0.30

Cook's distance 16

0.20 0.15

hatvalues(mod3)

0.2

32

0.0

0.10

0.1

Cook's distance

0.25

0.3

14

0

5

10

15

20

25

30

0

5

10

15

25

30

20

25

30

Covratio

1.0

1.2

1.4

1.0 0.5 0.0 −0.5

0.8

−1.0

DFFITS

20

Index

1.6

Obs. number

0

5

10

15

20

25

30

0

5

10

15

Index

Figura 7: An´alisis de Diagnostico

Index

0.0 −0.6

−0.6

−0.4

−0.2

DFBETAS para X1

0.0 −0.2 −0.4

DFBETAS para el intercepto

0.2

0.2

0.4

0.4

0.6

110

0

5

10

15

20

25

30

0

5

10

15

20

25

30

20

25

30

Index

0.0

DFBETAS para interacción

0.2 0.0 −0.4

−0.5

−0.2

DFBETAS para X2

0.4

0.5

0.6

Index

0

5

10

15

20

25

30

0

5

10

15

Index

Index

Figura 8: An´alisis de Diagnostico la precisi´on en las estimaciones de los par´ametros del modelo, por tanto se ajustar´a un nuevo modelo eliminando dichas observaciones. y el ajuste del nuevo modelo seria datos1=as.data.frame(data[-c(14,16),]);datos1 reg1=lm(y~x1+x2+I(x1*x2),datos1) ;reg1 Call: lm(formula = y ~ x1 + x2 + I(x1 * x2), data = datos1) Coefficients: (Intercept) 5.2529

x1 4.8584

x2 -0.7046

I(x1 * x2) -0.9175

Veamos unas estad´ısticas importantes para el modelo summary(reg1)

111

Call: lm(formula = y ~ x1 + x2 + I(x1 * x2), data = datos1) Residuals: Min 1Q -1.86269 -0.82783

Median 0.05471

3Q 0.46036

Max 2.03126

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.2529 1.9729 2.662 0.0131 * x1 4.8584 3.0132 1.612 0.1190 x2 -0.7046 0.4611 -1.528 0.1385 I(x1 * x2) -0.9175 0.7405 -1.239 0.2264 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.048 on 26 degrees of freedom Multiple R-squared: 0.4986,Adjusted R-squared: 0.4407 F-statistic: 8.618 on 3 and 26 DF, p-value: 0.0003872 Notamos que nuestro modelo mejoro el ajuste al eliminar las observaciones ya que ahora esta ajustando un 49.86 % y antes su ajuste era de 44.48 % Ahora calculemos el anova del modelo anova(reg1) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x1 1 17.2368 17.2368 15.6807 0.000519 *** x2 1 9.4940 9.4940 8.6369 0.006824 ** I(x1 * x2) 1 1.6877 1.6877 1.5353 0.226383 Residuals 26 28.5800 1.0992 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Ahora hagamos la validaci´on de supuestos del modelo ajustado para  Normalidad Miremos el gr´afico de probabilidad normal para los re-

112 siduales del modelo para verificar la normalidad para los errores del

0 −2

−1

resid(reg1)

1

2

Norm alidad resi

−2

−1

0

1

2

norm quantiles

Figura 9: Gr´afico de probabilidad para los residuales modelo ajustado usamos el test de Shapiro-Wilk teniendo como resultado shapiro.test(reg1$res) Shapiro-Wilk normality test data: reg1$res W = 0.97326, p-value = 0.6315 a un nivel de significancia del 0.05 no se rechaza la hip´otesis nula ya que el p-valor = 0.6315¿0.05 Por tanto los residuales del modelo siguen una distribuci´on normal

113 Independencia para verificar la independencia para los errores del modelo ajustado usamos el test de Durbin-Watson durbinWatsonTest(reg1) lag Autocorrelation D-W Statistic p-value 1 -0.01462423 1.876693 0.408 Alternative hypothesis: rho != 0 Homogeneidad de Varianzas para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de BreuschPagan de la librer´ıa lmtest , teniendo como resultado bptest(reg1) studentized Breusch-Pagan test data: reg1 BP = 2.4806, df = 3, p-value = 0.4788 De acuerdo con los resultados de la prueba, el p-valor es 0.4788 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´otesis nula, es decir que el modelo no presenta problemas de homocedasticidad. De lo anterior observamos que eliminando las observaciones influyentes la estimaci´on de los par´ametros y los resultados de la validaci´on de supuestos no cambia demasiado, solo mejora un poco la estimacion de los β Ahora calculemos el modelo para cada grupo, si X1 = 0 tenemos Yˆ = 5,2529 + 4,7445X1 − 0,7046X2 − 0,8866X1 X2 = 5,2529 + 4,7445(0) − 0,7046X2 − 0,8866(0)X2 = 5,2529 − 0,7046X2

114 X1 = 1 Yˆ = = = =

5,2529 + 4,7445X1 − 0,7046X2 − 0,8866X1 X2 5,2529 + 4,7445(1) − 0,7046X2 − 0,8866(1)X2 5,2529 + 4,7445 − (0,7046 + 0,8866)X2 9,9974 − 1,5912X2

As´ı cuando tomamos X1 = 0 es decir antigeno calla en la superficie de los blastos (Ag + ) tenemos que el intercepto para el modelo es de 5.2529 , y el par´ametro β2 nos indica que si no hay presencia de antigeno Ag − tenemos una variabilidad de -1.5912 en el modelo; mientras que si X1 = 1 es decir hay presencia del antigenos es decir (Ag − ) El intercepto aumenta a 9.9974, y la variabilidad del modelo logra disminuir a -1.5912 5. a) Para los datos datos obtenidos en un programa de entrenamiento en • Gram´atica (G) • habilidades lectoras (R) • ortograf´ıa (S) los puntajes obtenidos para G experimental es data1=data.frame(geG,geR,geS);data1 geG geR geS 1 31 12 24 2 52 64 32 3 57 42 21 4 63 19 54 5 42 12 41 6 71 79 64 7 65 38 52 8 60 14 58 9 54 75 69 10 67 22 69 realizaremos un an´alisis de regresi´on para el grupo G. experimental, donde definimos a Y=G, X1 = R y X2 = S Ahora bien ajustando el modelo en R tenemos la siguiente salida modg=lm(y~x1+x2) > modg

115

Call: lm(formula = y ~ x1 + x2) Coefficients: (Intercept) x1 x2 32.7327 0.1008 0.4063 realicemos un summary para ver un resumen de el modelo y estad´ısticas importantes summary(modg) Call: lm(formula = y ~ x1 + x2) Residuals: Min 1Q -14.330 -6.499

Median 3.150

3Q 5.881

Max 11.501

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 32.7327 9.9476 3.291 0.0133 * x1 0.1008 0.1306 0.772 0.4655 x2 0.4063 0.1923 2.113 0.0725 . --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.05 on 7 degrees of freedom Multiple R-squared: 0.4698,Adjusted R-squared: 0.3183 F-statistic: 3.101 on 2 and 7 DF, p-value: 0.1085 Notamos que el modelo para G.Experimental ajusta un 46.98 % ademas notamos que solo nos dio significativo el intercepto ahora bien veamos los intervalos de confianza del 95 % para los par´ametros estimados estos son: confint(modg) 2.5 % 97.5 % (Intercept) 9.21042598 56.2549982 x1 -0.20801330 0.4096021

116 x2 -0.04842817 0.8611281 claramente vemos que X1 , X2 contienen al cero esto nos confirma el hecho que no est´en siendo significativos realicemos el anova del modelo anova(modg) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x1 1 175.59 175.59 1.7383 0.22886 x2 1 450.92 450.92 4.4640 0.07249 . Residuals 7 707.09 101.01 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 An´ alisis de Diagnostico para el primer Grupo veamos que observaciones estan siendo imfluyentes en el modelo influence.measures(modg) Influence measures of lm(formula = y ~ x1 + x2) : dfb.1_ dfb.x1 dfb.x2 dffit cov.r cook.d 1 -1.28973 0.49937 0.8709 -1.3413 0.666 0.452595 2 -0.00720 -0.01134 0.0108 -0.0165 2.507 0.000105 3 1.07705 0.38107 -1.1264 1.3116 0.867 0.463163 4 0.05932 -0.20152 0.1242 0.3208 1.558 0.037178 5 -0.32994 0.32247 0.0659 -0.4914 1.311 0.081474 6 -0.19412 0.31371 0.1144 0.4268 2.316 0.067765 7 0.03434 -0.00995 0.0534 0.2541 1.362 0.022995 8 0.00971 -0.10021 0.0734 0.1439 2.071 0.007977 9 1.07338 -1.18969 -0.8473 -1.9260 0.365 0.746119 10 -0.09389 -0.16912 0.2509 0.3280 2.153 0.040394 realicemos los gr´aficos para verificar los datos influyentes notamos que lasobservaciones 1,3,9 disminuyen la precisi´on en las estimaciones de los par´ametros del modelo, realizaremos la validaci´on de los supuestos y verificaremos si es necesario eliminar las observaciones .hagamos la validaci´on de supuestos del modelo ajustado para 

hat inf 0.355 * 0.367 * 0.392 * 0.183 0.209 0.400 * 0.105 0.256 0.398 * 0.336

117

0.8

0.40

Cook's distance

0.30 0.10

0.0

0.15

0.2

0.20

0.25

hatvalues(modg)

3

1

0.4

Cook's distance

0.6

0.35

9

2

4

6

8

10

2

4

6

8

10

8

10

Index

1.0

1.5

Covratio

−0.5 −1.0

0.5

−1.5 −2.0

DFFITS

0.0

0.5

2.0

1.0

2.5

Obs. number

2

4

6

8

10

2

Index

4

6 Index

Figura 10: An´alisis de Diagnostico Grupo 1 Normalidad Miremos el gr´afico de probabilidad normal para los residuales del modelo En el gr´afico vemos que que no hay datos at´ıpicos y ninguno se sale de las bandas de ajuste por lo que podemos deducir que los residuos se distribuyen normal para verificar la normalidad para los errores del modelo ajustado usamos el test de Shapiro-Wilk teniendo como resultado: shapiro.test(modg$res) Shapiro-Wilk normality test data: modg$res W = 0.90075, p-value = 0.2233

0.0 −0.5

DFBETAS para X1

0.0

−1.0

−1.0

−0.5

DFBETAS para el intercepto

0.5

1.0

0.5

118

2

4

6

8

10

2

4

6

8

10

Index

0.0 −0.5 −1.0

DFBETAS para X2

0.5

Index

2

4

6

8

10

Index

Figura 11: An´alisis de Diagnostico Grupo 1 a un nivel de significancia del 0.05 no se rechaza la hip´otesis nula ya que el p − valor = 0,2233 > 0,05 Por tanto los residuales del modelo siguen una distribuci´on normal Independencia realizaremos la prueba de Durbin-Watson para realizar esta prueba durbinWatsonTest(modg) lag Autocorrelation D-W Statistic p-value 1 -0.08532137 1.919971 0.77 Alternative hypothesis: rho != 0 El p-valor de la prueba es de 0,77 > 0,05. lo que indica que s e acepta la hip´otesis nula, es decir que se cumple el supuesto de que los errores son independientes. Luego para reafirmar la independencia de los residuales usaremos la prueba de rachas, los resultados son los siguientes

119

Figura 12: gr´afico de probabilidad normal para los residuales del modelo del Grupo 1 library(tseries) X=as.factor(sign(resid(modg))) runs.test(X) Runs Test data: X Standard Normal = 0.14049, p-value = 0.8883 alternative hypothesis: two.sided lo que nos confirma que los residuales son independientes Homogenidad de varianza Miremos el gr´afico de dispersi´on de los residuales del modelo para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de Breusch-Pagan de la libreria lmtest , teniendo como resultado

120

−15

−10

−5

Residuales

0

5

10

Grafico de residuales

45

50

55

60

65

Ajustados

Figura 13: gr´afico de residuales del modelo para el Grupo 1 studentized Breusch-Pagan test data: modg BP = 0.86387, df = 2, p-value = 0.6493

De acuerdo con los resultados de la prueba, el p-valor es 0.6493 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´otesis nula, es decir que el modelo no presenta problemas de homocedasticidad.No vamos a eliminar las observaciones dado que el ajuste del modelo no cambia y ademas los supuestos se siguen cumpliendo igualmente y dado que la idea no es eliminar observaciones trabajaremos con nuestro modelo original

121 Ahora bien realicemos el an´alisis de regresi´on para el grupo G.control Asi, data1=data.frame(gcG,gcR,gcS);data1 gcG gcR gcS 1 70 34 24 2 31 50 20 3 60 40 15 4 65 36 12 5 70 29 18 6 78 48 24 7 90 47 26 8 98 18 40 9 95 10 10 donde tenemos que y=gcG x1=gcR x2=gcS modg=lm(y~x1+x2) modg Call: lm(formula = y ~ x1 + x2) Coefficients: (Intercept) x1 x2 87.7411 -0.9233 0.8222 Veamos unas estad´ısticas importantes para el modelo summary(modg) Call: lm(formula = y ~ x1 + x2) Residuals: Min 1Q -27.021 -6.010

Median -3.143

3Q 8.270

Max 24.276

122 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 87.7411 21.2023 4.138 0.00609 ** x1 -0.9233 0.4358 -2.119 0.07842 . x2 0.8222 0.6641 1.238 0.26190 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 16.95 on 6 degrees of freedom Multiple R-squared: 0.4983,Adjusted R-squared: 0.3311 F-statistic: 2.98 on 2 and 6 DF, p-value: 0.1263 notamos que el ajuste del modelo explica un 49.83 %, la variabilidad, el valor del coeficiente de determinaci´on es de 33.11 % nos es un buen ajuste. Luego Los intervalos de confianza del 95 % para los par´ametros estimados son: confint(modg) 2.5 % 97.5 % (Intercept) 35.8610107 139.6212815 x1 -1.9895999 0.1430128 x2 -0.8026715 2.4470873 ahora calculemos un an´alisis de varianza del modelo anova(modg) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x1 1 1272.5 1272.49 4.4265 0.08004 . x2 1 440.7 440.70 1.5331 0.26190 Residuals 6 1724.8 287.47 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 An´ alisis de Diagnostico El modelo que se esta estudiando es yˆ = 87,7411 − 0,9233X1 + 0,8222X2

123 Realicemos un an´alisis de diagnostico para estudiar la presencia de observaciones o conjunto de observaciones sobre la estimaci´on de los par´ametros y supuestos del modelo, para ello usaremos las medidas estad´ısticas propuestas para identificar y medir el conjunto de observaciones influyentes. Usaremos R-studio para ello as´ı. Influence measures of lm(formula = y ~ x1 + x2) : dfb.1_ dfb.x1 dfb.x2 dffit cov.r cook.d hat inf 1 -0.00935 0.00703 -0.0446 -0.1341 1.834 0.007022 0.125 2 0.50045 -1.20850 0.1345 -1.5866 0.177 0.424350 0.268 * 3 -0.03020 -0.02918 0.0494 -0.0899 2.078 0.003212 0.186 4 0.01307 0.00172 -0.0157 0.0217 2.262 0.000188 0.237 5 -0.09840 0.05301 0.0425 -0.1403 1.889 0.007696 0.146 6 -0.26441 0.39427 0.1310 0.5676 1.310 0.107165 0.241 7 -0.61608 0.74095 0.4522 1.1738 0.373 0.300581 0.248 8 0.22466 1.06375 -1.8225 -2.2504 7.491 1.730960 0.856 * 9 1.26279 -0.97402 -0.6555 1.2887 3.711 0.578739 0.692 * Ahora para verificar los valores de los estad´ısticos, los podemos representar gr´aficamente para apreciar mejor los puntos influyentes: Un nuevo modelo eliminando estas observaciones viene dado por yˆ = 32,14980,3767X1 + 1,2769X2 Validaci´ on de Supuestos Ahora se realizara la validaci´on de supuestos del modelo ajustado para Linealidad Para el ajuste para  Linealidad Para el ajuste obtenido con los datos yˆ = 87,7411 − 0,9233X1 + 0,8222X2 Concideremos los gr´aficos de regresi´on parcial para ver la relaci´on lineal entre regresores y respuesta Observamos que la tendencia de residuos parciales no es de tipo lineal, realizaremos la prueba de linealidad para corroborar este indicio

124

0.30

Cook's distance 16

0.20

hatvalues(mod3)

0.15

0.2

32

0.0

0.10

0.1

Cook's distance

0.25

0.3

14

0

5

10

15

20

25

30

0

5

10

20

25

30

20

25

30

Index

Covratio

1.0

1.2

1.4

1.0 0.5 0.0 −0.5

0.8

−1.0

DFFITS

15

1.6

Obs. number

0

5

10

15 Index

20

25

30

0

5

10

15 Index

Figura 14: Analisis de Diagnostico del Grupo 2 Normalidad Observemos el gr´afico de probabilidad normal para los residuales del modelo Donde observamos que la mayor´ıa de residuos del modelo siguen una distribuci´on normal. Prueba de Shapiro Wilk Para verificar la normalidad para los errores del modelo ajustado usamos el test de Shapiro-Wilk teniendo como resultado Shapiro-Wilk normality test data: resid(mod) W = 0.94725, p-value = 0.66 Note que nuestro p-value es de (p − value = 0,66 > 0,05) entonces no se rechaza la hip´otesis de normalidad, esto es, los errores del modelo siguen una distribuci´on normal Independencia Para verificar la independencia para los errores del modelo ajustado usamos el test de Durbin -Watson de la librer´ıa CAR

−0.5

−1.0

−0.5

0.0

DFBETAS para X1

0.5 0.0

DFBETAS para el intercepto

0.5

1.0

1.0

125

2

4

6

8

2

4

6

8

Index

−0.5 −1.0 −1.5

DFBETAS para X2

0.0

0.5

Index

2

4

6

8

Index

Figura 15: An´alisis de Diagnostico del Grupo 2 , teniendo como resultado durbinWatsonTest(mod) lag Autocorrelation D-W Statistic p-value 1 0.1871534 1.564596 0.342 Alternative hypothesis: rho != 0 Notamos que el p-value es de 0,342 > 0,05 por tanto se rechaza la hip´otesis nula, es decir que no se cumple el supuesto de que los errores son independientes. Para corroborar la independencia de los residuales usaremos la prueba de rachas, esta se implementa con la funci´on runs.test() de la librer´ıa tseries Runs Test data: X Standard Normal = 0.40161, p-value = 0.688

30 20

beta*x2+res

−10

0

−60

10

−40

beta*x1+res

−20

40

0

126

10

20

30

40

50

10

15

20

x1

25

30

35

40

x2

Figura 16: Linealidad

0 −10 −20

resid(mod)

10

20

Normalidad resi

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

norm quantiles

Figura 17: Normalidad Grupo 2 alternative hypothesis: two.sided Notamos que en esta prueba el p-value resulta superior a 0.05 entonces concluimos que existe independencia en los errores

127

0 −10 −20

residuals(mod)

10

20

Homogenidad de varianza Miremos el gr´afico de dispersi´on de los residuales del modelo En el grafico parece que la

2

4

6

8

Index

Figura 18: Dispersi´on de los residuos varianza de los errores no tiene ning´ un patr´on a crecer o decrecer, para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de Breusch-Pagan de la libreria lmtest , teniendo como resultado Breusch-Pagan test data: mod BP = 3.3721, df = 2, p-value = 0.1852 Lo que confirma lo dicho anteriormente, con un p-valor de 0.1852 el cual es mayor que 0.05 lo que lleva al no rechazo de la hipotesis nula, es decir que el modelo no presenta problemas de heterocedasticidad

128 b) Realizaremos un an´alisis de regresi´on del tipo Yt = β0 + β1 Xt1 + β2 Xt2 + αIt +  Donde It =

  1 si i ∈ G. experimental 

0 si i ∈ G.Control

Ahora bien un encabezado de los datos a˜ nadiendo la v. cualitativa es y x1 x2 I 1 31 12 24 1 2 52 64 32 1 3 57 42 21 1 4 63 19 54 1 5 42 12 41 1 6 71 79 64 1 Ajustando el modelo en R tenemos la siguiente salida modg1 Call: lm(formula = y ~ x1 + x2 + I) Coefficients: (Intercept) x1 x2 65.6103 -0.1155 0.5425 entonces un modelo estimado por mco es

I -31.3141

yt = 65,6103 − −0,1155X1 + 0,5425X2 − 31,3141I Realicemos un summary para ver unas estad´ısticas importantes summary(modg1) Call: lm(formula = y ~ x1 + x2 + I)

129 Residuals: Min 1Q -39.687 -6.887

Median -2.027

3Q 9.420

Max 25.119

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 65.6103 9.2753 7.074 3.79e-06 *** x1 -0.1155 0.1827 -0.632 0.53688 x2 0.5425 0.2708 2.003 0.06354 . I -31.3141 10.3222 -3.034 0.00838 ** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.82 on 15 degrees of freedom Multiple R-squared: 0.3856,Adjusted R-squared: 0.2628 F-statistic: 3.139 on 3 and 15 DF, p-value: 0.05667 Notamos que el modelo ajusta un 38,56 % vemos que los par´ametros mas significativos son β0 la v.cualitativa(It ) y un poco el par´ametro β2 Un intervalo de confianza para el 95 % viene dado por confint(modg1) 2.5 % 97.5 % (Intercept) 45.84044944 85.3801653 x1 -0.50484528 0.2739206 x2 -0.03465486 1.1196442 I -53.31539961 -9.3128358 Donde vemos que el parametro β0 no contiene al cero asi como el parametroα que es el parametro para la v. cualitativa Realicemos el an´alisis de varianza para el modelo anova(modg1) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x1 1 42.7 42.71 0.1707 0.685328 x2 1 10.5 10.54 0.0421 0.840135 I 1 2302.5 2302.49 9.2031 0.008378 **

130 Residuals 15 3752.8 250.19 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Vemos otra ves que le par´ametro α es significativo para el modelo Ahora hagamos un an´alisis de diagnostico para estudiar la presencia de observaciones o conjuntos de observaciones sobre la estimaci´on de los par´ametros y supuestos del modelo ´ ANALISIS DE DIAGNOSTICO Uusaremos las medidas estad´ısticas propuestas para identificar y medir el conjunto de observaciones influyentes. Usaremos R-studio para ello as´ı. influence.measures(modg1) Influence measures of lm(formula = y ~ x1 + x2 + I) : dfb.1_ dfb.x1 dfb.x2 dfb.I dffit cov.r 1 -0.50963 0.29809 0.49929 -0.58791 -0.7786 1.337 2 -0.02499 0.24685 -0.23415 0.26037 0.3794 1.728 3 0.30659 0.22190 -0.74723 0.76429 0.8987 1.270 4 0.00884 -0.02752 0.01622 0.00777 0.0474 1.571 5 -0.23201 0.27738 0.06951 -0.22122 -0.4440 1.286 6 -0.42290 0.45343 0.18482 0.00656 0.6309 1.645 7 -0.01484 -0.00411 0.02878 0.05190 0.1521 1.389 8 -0.02446 0.09961 -0.07104 -0.00414 -0.1525 1.650 9 0.35607 -0.31315 -0.23204 0.05075 -0.5129 1.768 10 0.01246 0.04678 -0.07241 0.02097 -0.0991 1.824 11 -0.05507 0.00580 -0.01721 0.06746 -0.1098 1.448 12 -0.38936 -0.71580 0.20647 0.56532 -1.5049 0.104 13 -0.13214 -0.05302 0.07440 0.05965 -0.2317 1.363 14 -0.05340 -0.00961 0.03296 0.01254 -0.0774 1.509 15 -0.03517 0.00776 0.00539 0.01959 -0.0479 1.485 16 0.02727 0.05005 0.00730 -0.06716 0.1275 1.478 17 0.07896 0.14473 0.06627 -0.25166 0.4218 1.109 18 0.15841 -0.28624 0.40510 -0.46569 0.5887 1.425 19 0.95759 -0.53963 -0.28775 -0.13507 0.9900 0.641 Ahora para verificar los valores de los estad´ısticos, los podemos representar gr´aficamente para apreciar mejor los puntos

cook.d 0.148326 0.037677 0.193800 0.000602 0.049798 0.101244 0.006110 0.006192 0.068082 0.002625 0.003210 0.309202 0.014012 0.001599 0.000614 0.004322 0.044039 0.087247 0.206889

hat i 0.314 0.305 0.332 0.164 0.189 0.349 0.104 0.220 0.350 0.282 0.114 0.144 0.128 0.136 0.117 0.134 0.134 0.279 0.207

131 influyentes: si eliminamos estas observaciones veamos el mo0.35

Cook's distance

0.25 0.20

hatvalues(modg1)

0.15

0.20

19 3

0.10

0.00

0.05

0.15

0.10

Cook's distance

0.25

0.30

0.30

12

5

10

15

5

10

15

Index

1.0

Covratio

0.5

−0.5 −1.5

−1.0

DFFITS

0.0

0.5

1.5

1.0

Obs. number

5

10

15

5

10

15

Index

0.0 −0.2

DFBETAS para X1

−0.5

−0.6

−0.4

0.0

DFBETAS para el intercepto

0.5

0.2

0.4

1.0

Index

5

10

15

5

10

15

Index

0.4 0.2 −0.2

0.0

DFBETAS para It

0.0 −0.2

−0.6

−0.6

−0.4

−0.4

DFBETAS para X2

0.2

0.6

0.4

0.8

Index

5

10 Index

15

5

10

15

Index

Figura 19: Influencia Notamos que la incluci´on de las observaciones 3,12,19 disminuyen la precisi´on en las estimaciones delo nos queda

132 Call: lm(formula = y ~ x1 + x2 + I, data = dat, x = T) Coefficients: (Intercept) x1 x2 I 57.71647 0.01859 0.76918 -41.86728 Realicemos el summary para ver el grado de ajuste y algunas estadisticas importantes summary(modn) Call: lm(formula = y ~ x1 + x2 + I, data = dat, x = T) Residuals: Min 1Q -16.317 -4.051

Median -1.411

3Q 6.058

Max 11.411

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 57.71647 5.93938 9.718 4.88e-07 *** x1 0.01859 0.10740 0.173 0.865430 x2 0.76918 0.17615 4.367 0.000918 *** I -41.86728 6.68682 -6.261 4.18e-05 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.732 on 12 degrees of freedom Multiple R-squared: 0.7713,Adjusted R-squared: 0.7141 F-statistic: 13.49 on 3 and 12 DF, p-value: 0.0003762 Notamos que al eliminar estas observaciones el grado de ajuste cambia muy significativamente a 77,13 % y los par´ametros β0 , α, β2 pasan hacer muy significativos Realicemos la validaci´on de supuestos para el modelo Validaci´ on de Supuestos Linealidad Notamos que la tendencia de residuos parciales no es de tipo lineal,por lo caul podemos tener indicos de problemas con la hip´otesis de linealidad. Normalidad Observemos el gr´afico de probabilidad normal

133

Figura 20: Linealidad

−40

−30

−20

−10

e

0

10

20

para los residuales del modelo Donde observamos que la ma-

−2

−1

0

1

2

norm quantiles

Figura 21: Grafico de Normalidad de los residuos yor´ıa de residuos del modelo siguen una distribuci´on normal Prueba de Shapiro Wilk Para verificar la normalidad para los

134 errores del modelo ajustado usamos el test de Shapiro-Wilk teniendo como resultado Shapiro-Wilk normality test data: modg1$res W = 0.94711, p-value = 0.3525 Note que nuestro p-value es de (p − value = 0,3525 > 0,05) entonces no se rechaza la hip´otesis de normalidad, esto es, los errores del modelo siguen una distribuc´on normal Independencia Para veriFIcar la independecia para los errores del modelo ajustado usamos el test de Durbin -Watson durbinWatsonTest(modg1) lag Autocorrelation D-W Statistic p-value 1 0.3076459 1.15717 0.014 Alternative hypothesis: rho != 0 El p-valor de la prueba es de 0,014 < 0,05 . Se rechaza la hip´otesis nula, es decir que no se cumple el supuesto de que los errores son independientes.Ahora usamos la prueba de ra rachas para reafirmar lo antes dicho Runs Test data: X Standard Normal = -2.1184, p-value = 0.03414 alternative hypothesis: two.sided Notamos que en esta prueba el p-value resulta menor superior a 0.05 entonces podemos decir que no existe independencia en los errores Homogenidad de varianza Miremos el gr´afico de dispersi´on de los residuales del modelo En el gr´afico parece que la varianza de los errores no tiene ning´ un patr´on a crecer o decrecer, para verificar la Homogeneidad de Varianzas para los errores del modelo ajustado usamos el test de Breusch-Pagan de la libreria lmtest , teniendo como resultado bptest(modg1)

135

0 −2

−1

Residuales

1

2

Grafico de residuales

2

3

4

5

Ajustados

Figura 22: Grafico de Dispersi´on studentized Breusch-Pagan test data: modg1 BP = 2.3461, df = 3, p-value = 0.5037 que confirma lo dicho anteriormente, con un p-valor de 0.5037 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´otesis nula, es decir que el modelo no presenta problemas de heterocedasticidad C) Para el analisis del primer grupo el modelo fue el siguiente modelo Grupo 1 yˆ = 32,7327 + 0,2008X1 + 0,4063X2 Modelo de grupo 2 yˆ = 87,7411 − 0,9233X1 + 0,8222X2 Modelo con v. Dicotomica yˆ = 65,6103 − 0,1155X1 + 0,5425X2 − 31,3141 Donde el que mejor ajusto fue l modelo 3 con la incluision de la variable dic´otoma, los parametros para este fueron significativos casi todos mientras que cuando miramos los modelos

136 por separado el grado de ajuste era bajo y los parametros excepto el β0 no eran significativos d) La respueta a la pregunta ¿ Es necesario un termino de interacci´on Ii Xi1 o Ii Xi2 ? es si es necesario ya estamos analizando los puntajes obtenidos en un programa de entrenamiento en gram´atica (G), habilidades lectoras (R) y ortograf´ıa (S). esto es que la relaci´on entre los Gram´atica y habilidades lectoras y ortograf´ıa no son aditivas es decir existe existe interrelaci´on dado que la gram´atica aumenta a medida que si tienen una buenas habilidades lectoras y una buena ortograf´ıa 6. Definamos Y = R(Dolares), X = E,  1 Si es pais de Europa I1 = 0 Si es pais de America del Sur Se ajusta el modelo lineal Yi = β0 + β1 Xi + β2 Ii + i Estimando las variables en R tenemos tenemos el modelo estimado y la significancia de las variables: > mod1 mod1 Call: lm(formula = E ~ R + I) Coefficients: (Intercept) R 64.93053 0.00065

I 5.036086

> res1 res1$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 6.493053e+01 1.586391810 40.929695 1.545271e-11 R 6.511852e-04 0.000272796 2.387077 4.075219e-02 I 5.036086e+00 2.038104824 2.470965 3.551400e-02

137 as´ı llegamos al modelo: Yˆi = 64,93053 + 0,00065Xi + 5,036086Ii donde βˆ0 = 64,93053 es el intercepto cuando Ii = 0 es decir cuando los pa´ıses son de Am´erica del sur. β1 = 0,00065 es la variabilidad del modelo y β2 = 5,036086 es la diferencia entre el los pa´ıses de Europa y los pa´ıses de am´erica del sur, ademas todos los par´ametros son significativos a) ahora se ajusta el modelo lineal Yi = β0 + β1 Xi + β2 Ii + β3 Ii Xi + i Estimando las variables en R tenemos tenemos el modelo estimado y la significancia de las variables: > mod2 mod2 Call: lm(formula = E ~ R + I + I(I * R)) Coefficients: (Intercept) 62.177238

R 0.001518

> res2 res2$coefficients Estimate (Intercept) 62.177238399 R 0.001517546 I 9.653236816 I(I * R) -0.001129937

I 9.653237

I(I * R) -0.001130

Std. Error t value Pr(>|t|) 1.9229178860 32.334838 9.118731e-10 0.0004863504 3.120273 1.422281e-02 2.8689469041 3.364732 9.862320e-03 0.0005554269 -2.034357 7.633943e-02

por tanto tenemos el modelo: Yˆi = 62,17724 + 0,00152Xi + 9,65324Ii − 0,00113Ii Xi Donde notamos que el nuevo termino incluido no es significativo.

138 b) Ahora compararemos los dos modelos con varios criterios: > res1$r.squared [1] 0.732494 > AIC(mod1) [1] 64.7337 > BIC(mod1) [1] 66.67333 > require(DAAG) > press(mod1) [1] 141.2235 > res2$r.squared [1] 0.8236991 > AIC(mod2) [1] 61.73031 > BIC(mod2) [1] 64.15484 > require(DAAG) > press(mod2) [1] 115.3015 Finalmente con todos los criterios el modelo dos es mejor por tanto se recomienda trabajar con el modelo dos. 7. Primeramente definimos las variables Ii1 y Ii2 como siguen:  1 si el metodo de lavado de suelos es v1 Ii1 = . 0 otro caso

(3)

y  1 si el metodo de lavado de suelos es v2 Ii2 = . 0 otro caso Las cuales definimos respectivamente en R as´ı: I1=ifelse(VOL=="v1", 1, 0) I2=ifelse(VOL=="v2", 1, 0) Tenemos el modelo ajustado: Yˆi = βˆ0 + βˆ1 Di + βˆ2 Di2 + βˆ3 Di3 + βˆ4 Di4 + βˆ5 Ii1 + βˆ6 Ii2

(4)

139 ahora realiz´andolo en R llegamos al modelo y su respectivo resumen > modd summary(modd)

I(Dosis^2) -3.25912 I2 -5.33333

I(Dosis^3) 0.37228

Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I(Dosis^3) + I(Dosis^4) + I1 + I2) Residuals: Min 1Q -31.083 -14.514

Median -2.829

3Q 13.597

Max 48.890

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 47.91358 7.00140 6.843 3.28e-09 *** Dosis 11.70902 11.17388 1.048 0.298569 I(Dosis^2) -3.25912 5.33704 -0.611 0.543554 I(Dosis^3) 0.37228 0.84207 0.442 0.659889 I(Dosis^4) -0.01706 0.04184 -0.408 0.684822 I1 -24.32667 6.07243 -4.006 0.000161 *** I2 -5.33333 6.07243 -0.878 0.383023 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 21.04 on 65 degrees of freedom Multiple R-squared: 0.2957,Adjusted R-squared: 0.2307 F-statistic: 4.548 on 6 and 65 DF, p-value: 0.0006571

140 ahora introducimos en el modelo I1 D y I2 D Yˆi = βˆ0 + βˆ1 Di + βˆ2 Di2 + βˆ3 Di3 + βˆ4 Di4 + βˆ5 Ii1 + βˆ6 Ii2 + βˆ7 Ii1 Di + βˆ8 Ii2 Di para lo cual tenemos > modd1 summary(modd1) Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I(Dosis^3) + I(Dosis^4) + I1 + I2 + I(Dosis * I1) + I(Dosis * I2)) Residuals: Min 1Q -37.488 -15.468

Median -2.787

3Q 12.354

Max 40.095

Coefficients: (Intercept) Dosis I(Dosis^2) I(Dosis^3) I(Dosis^4)

Estimate Std. Error t value Pr(>|t|) 35.81358 8.00225 4.475 3.26e-05 *** 14.12902 10.34367 1.366 0.176808 -3.25912 4.91979 -0.662 0.510099 0.37228 0.77624 0.480 0.633181 -0.01706 0.03857 -0.442 0.659791

141 I1 -16.97881 9.92364 -1.711 0.092011 . I2 23.61881 9.92364 2.380 0.020349 * I(Dosis * I1) -1.46957 1.63883 -0.897 0.373284 I(Dosis * I2) -5.79043 1.63883 -3.533 0.000775 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 19.39 on 63 degrees of freedom Multiple R-squared: 0.4199,Adjusted R-squared: 0.3462 F-statistic: 5.7 on 8 and 63 DF, p-value: 1.873e-05 Note que el ajuste aumento y el termino I1 D resulto altamente significativo. ahora introducimos en el modelo I1 D2 y I2 D2 Yˆi =βˆ0 + βˆ1 Di + βˆ2 Di2 + βˆ3 Di3 + βˆ4 Di4 + βˆ5 Ii1 + βˆ6 Ii2 + βˆ7 Ii1 Di + βˆ8 Ii2 D + βˆ9 Ii1 D2 + βˆ10 Ii2 D2 i

i

para lo cual tenemos > modd2 modd2 Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I(Dosis^3) + I(Dosis^4) + I1 + I2 + I(Dosis * I1) + I(Dosis * I2) + I((Dosis^2) * I1) + I((Dosis^2) * I2)) Coefficients: (Intercept) 20.76730 I(Dosis^3) 0.37228 I2 46.61866 I((Dosis^2) * I1)

Dosis 25.41373 I(Dosis^4) -0.01706 I(Dosis * I1) -18.07381 I((Dosis^2) * I2)

I(Dosis^2) -4.38759 I1 5.16018 I(Dosis * I2) -23.04032

142 1.66042

1.72499

> summary(modd2) Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I(Dosis^3) + I(Dosis^4) + I1 + I2 + I(Dosis * I1) + I(Dosis * I2) + I((Dosis^2) * I1) + I((Dosis^2) * I2)) Residuals: Min 1Q -35.898 -10.371

Median -1.605

3Q 11.489

Max 35.379

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 20.76730 8.30221 2.501 0.015069 * Dosis 25.41373 9.89605 2.568 0.012692 * I(Dosis^2) -4.38759 4.50154 -0.975 0.333565 I(Dosis^3) 0.37228 0.70871 0.525 0.601290 I(Dosis^4) -0.01706 0.03521 -0.484 0.629811 I1 5.16018 11.34603 0.455 0.650867 I2 46.61866 11.34603 4.109 0.000121 *** I(Dosis * I1) -18.07381 5.33620 -3.387 0.001242 ** I(Dosis * I2) -23.04032 5.33620 -4.318 5.89e-05 *** I((Dosis^2) * I1) 1.66042 0.51221 3.242 0.001928 ** I((Dosis^2) * I2) 1.72499 0.51221 3.368 0.001317 ** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 17.7 on 61 degrees of freedom Multiple R-squared: 0.5318,Adjusted R-squared: 0.455 F-statistic: 6.929 on 10 and 61 DF, p-value: 3.925e-07 Nuevamente aumenta el ajuste del modelo, los nuevos t´erminos I1 D2 y I2 D2 resultan altamente significativos y ademas resultan significativo D y I1 D que no eran significativos en el modelo

143 ahora introducimos en el modelo I1 D3 y I2 D3 Yˆi =βˆ0 + βˆ1 Di + βˆ2 Di2 + βˆ3 Di3 + βˆ4 Di4 + βˆ5 Ii1 + βˆ6 Ii2 + βˆ7 Ii1 Di + βˆ8 Ii2 D + βˆ9 Ii1 D2 + βˆ10 Ii2 D2 + βˆ11 Ii1 D3 + βˆ12 Ii2 D3 i

i

i

i

para lo cual tenemos > modd3 modd3 Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I(Dosis^3) + I(Dosis^4) + I1 + I2 + I(Dosis * I1) + I(Dosis * I2) + I((Dosis^2) * I1) + I((Dosis^2) * I2) + I((Dosis^3) * I1) + I((Dosis^3) * I2)) Coefficients: (Intercept) 19.62182 I(Dosis^3) 0.42000 I2 50.40179 I((Dosis^2) * I1) 1.44376 I((Dosis^3) * I2) -0.15763

Dosis 28.02926 I(Dosis^4) -0.01706 I(Dosis * I1) -17.28226 I((Dosis^2) * I2) 4.08944

I(Dosis^2) -5.10352 I1 4.81351 I(Dosis * I2) -31.67845 I((Dosis^3) * I1) 0.01444

> summary(modd3) Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I(Dosis^3) + I(Dosis^4) + I1 + I2 + I(Dosis * I1) + I(Dosis * I2) + I((Dosis^2) * I1) + I((Dosis^2) * I2) + I((Dosis^3) * I1) + I((Dosis^3) * I2)) Residuals: Min

1Q

Median

3Q

Max

144 -32.205 -10.971

-0.557

10.579

34.495

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 19.62182 8.80451 2.229 0.02966 * Dosis 28.02926 11.75337 2.385 0.02032 * I(Dosis^2) -5.10352 4.84803 -1.053 0.29677 I(Dosis^3) 0.42000 0.72382 0.580 0.56395 I(Dosis^4) -0.01706 0.03552 -0.480 0.63283 I1 4.81351 12.37499 0.389 0.69870 I2 50.40179 12.37499 4.073 0.00014 *** I(Dosis * I1) -17.28226 12.01865 -1.438 0.15573 I(Dosis * I2) -31.67845 12.01865 -2.636 0.01071 * I((Dosis^2) * I1) 1.44376 2.98642 0.483 0.63057 I((Dosis^2) * I2) 4.08944 2.98642 1.369 0.17608 I((Dosis^3) * I1) 0.01444 0.19609 0.074 0.94153 I((Dosis^3) * I2) -0.15763 0.19609 -0.804 0.42471 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 17.86 on 59 degrees of freedom Multiple R-squared: 0.5392,Adjusted R-squared: 0.4455 F-statistic: 5.753 on 12 and 59 DF, p-value: 1.878e-06 el aumento del modelo es muy poco por no decir nulo, ademas los par´ametros nuevos no resultan significativos e incluso par´ametros que eran significativos ahora no lo son. Finalmente introducimos en el modelo I1 D4 y I2 D4 Yˆi =βˆ0 + βˆ1 Di + βˆ2 Di2 + βˆ3 Di3 + βˆ4 Di4 + βˆ5 Ii1 + βˆ6 Ii2 + βˆ7 Ii1 Di + βˆ8 Ii2 D + βˆ9 Ii1 D2 + βˆ10 Ii2 D2 + βˆ11 Ii1 D3 + βˆ12 Ii2 D3 + βˆ13 Ii1 D4 + βˆ14 Ii2 D4 i

i

i

i

i

i

con lo cual llegamos finalmente a: > modd4 modd4

+I((Dosis^4)*I1)+I((Dosis^4)*I2))

Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I1 + I2 + I(Dosis * I1) + I(Dosis I((Dosis^2) * I2) + I((Dosis^3) * I((Dosis^4) * I1) + I((Dosis^4) *

I(Dosis^3) + I(Dosis^4) + * I2) + I((Dosis^2) * I1) + I1) + I((Dosis^3) * I2) + I2))

Coefficients: (Intercept) 19.78185 I(Dosis^3) 0.30332 I2 47.13652 I((Dosis^2) * I1) 13.91904 I((Dosis^3) * I2) 2.22329

I(Dosis^2) -4.38672 I1 7.59869 I(Dosis * I2) -4.46789 I((Dosis^3) * I1) -2.01641 I((Dosis^4) * I2) -0.11905

Dosis 26.69568 I(Dosis^4) -0.01122 I(Dosis * I1) -40.49208 I((Dosis^2) * I2) -10.53624 I((Dosis^4) * I1) 0.10154

> summary(modd4) Call: lm(formula = Y ~ Dosis + I(Dosis^2) + I1 + I2 + I(Dosis * I1) + I(Dosis I((Dosis^2) * I2) + I((Dosis^3) * I((Dosis^4) * I1) + I((Dosis^4) * Residuals: Min 1Q -41.521 -10.343

Median -0.022

3Q 11.193

I(Dosis^3) + I(Dosis^4) + * I2) + I((Dosis^2) * I1) + I1) + I((Dosis^3) * I2) + I2))

Max 31.549

Coefficients: (Intercept) Dosis I(Dosis^2)

Estimate Std. Error t value Pr(>|t|) 19.78185 8.55745 2.312 0.02443 * 26.69568 15.77789 1.692 0.09611 . -4.38672 7.53609 -0.582 0.56280

146 I(Dosis^3) 0.30332 I(Dosis^4) -0.01122 I1 7.59869 I2 47.13652 I(Dosis * I1) -40.49208 I(Dosis * I2) -4.46789 I((Dosis^2) * I1) 13.91904 I((Dosis^2) * I2) -10.53624 I((Dosis^3) * I1) -2.01641 I((Dosis^3) * I2) 2.22329 I((Dosis^4) * I1) 0.10154 I((Dosis^4) * I2) -0.11905 --Signif. codes: 0 ‘***’ 0.001

1.18904 0.05908 12.10206 12.10206 22.31331 22.31331 10.65763 10.65763 1.68155 1.68155 0.08355 0.08355

0.255 -0.190 0.628 3.895 -1.815 -0.200 1.306 -0.989 -1.199 1.322 1.215 -1.425

0.79957 0.85000 0.53259 0.00026 *** 0.07483 . 0.84201 0.19679 0.32703 0.23544 0.19139 0.22923 0.15965

‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.15 on 57 degrees of freedom Multiple R-squared: 0.5895,Adjusted R-squared: 0.4887 F-statistic: 5.847 on 14 and 57 DF, p-value: 6.859e-07 Note que el ajuste del modelo aumenta sin embargo los par´ametros en su mayor´ıa resultan muy significativos por lo cual se recomienda trabajar con el modelo: Yˆi =20,767 + 25,414Di − 4,388Di2 + 0,373Di3 − 0,017Di4 + 5,160Ii1 + 46,619Ii2 − 18,074Ii1 Di − 23,040Ii2 D + 1,660Ii1 Di2 + 1,725Ii2 Di2 el cual a pesar de tener un menor ajuste, sus par´ametros son en su mayor´ıa significativos

Modelo de Regresi´ on No Lineal 4

a) Del modelo Yi = β0 + β1 exp(β2 Xi ) + i A Partir de los par´ametros podemos decir que este modelo es no lineal en los par´ametros y en la variable Xi , ademas no se puede linealizar el modelo, por lo que el modelo es intr´ınsecamente no lineal. b) para hallar las ecuaciones de estimaci´on se tiene que:

147

f (β0 , β1 , β2 , X) = β0 + β1 exp(β2 Xi ) entonces i = Yi − β0 − β1 exp(β2 Xi ) y por tanto la funci´on de log-verosimilitud viene dada por: n n 1 X l(β0 , β1 , β2 , X) = − ·ln(2π)−n·ln(σ)− 2 (Yi − β0 − β1 exp(β2 Xi ))2 , 2 2σ i=1

para i = 1, 2, ..., n. Entonces las funciones de scores son: n 1 X ∂l(β0 , β1 , β2 , X) = 2 (Yi − β0 − β1 exp(β2 Xi )) ∂β0 σ i=1 n ∂l(β0 , β1 , β2 , X) 1 X (Yi − β0 − β1 exp(β2 Xi )) exp(β2 Xi ) = 2 ∂β1 σ i=1 n ∂l(β0 , β1 , β2 , X) 1 X =− 2 (Yi − β0 − β1 exp(β2 Xi )) β1 Xi exp(β2 Xi ) ∂β2 σ i=1 n ∂l(β0 , β1 , β2 , X) n 1 X =− + 3 (Yi − β0 − β1 exp(β2 Xi ))2 ∂σ σ σ i=1

As´ı tenemos jβ0 β0

n = 2 σ

jβ1 β1

n 1 X = 2 [exp(β2 Xi )]2 σ i=1

jβ2 β1

jβ2 β2

jβ1 β0

n 1 X = 2 exp(β2 Xi ) σ i=1

jβ2 β0

n 1 X = 2 β1 Xi exp(β2 Xi ) σ i=1

n 2 X = 2 [β1 Xi exp(β2 Xi )]2 + β1 Xi [exp(β2 Xi )]2 σ i=1

n n 1 X 2 X 2 = 2 [β1 Xi exp(β2 Xi )] jσβ0 = 3 (Yi − β0 − β1 exp(β2 Xi )) σ i=1 σ i=1

148

jσβ1

n 2 X = 3 (Yi − β0 − β1 exp(β2 Xi )) exp(β2 Xi ) σ i=1

jσβ2 =

n 2 X (Yi − β0 − β1 exp(β2 Xi )) β1 Xi exp(β2 Xi ) σ 3 i=1

y jσσ

n 2 X n (Yi − β0 − β1 exp(β2 Xi ))2 =− 2 + 4 σ σ i=1

c) As´ı tenemos que los elementos de la matriz de informaci´on de Fisher vienen dados por: n

iβ0 β0

X 1 exp(2β2 Xi ) = 2 (n), iβ1 β1 = σ i=1

iβ2 β0 =

n 1 X β1 Xi exp(β2 Xi ) σ 2 i=1

iσβ0 = 0 = iσβ1 , iσσ 5

iβ1 β0

2n = 2 σ

n 1 X exp(β2 Xi ) = 2 σ i=1

iβ2 β1 = − y iβ2 β2

n 1 X β1 Xi exp(2β2 Xi ) σ 2 i=1

n 1 X [β1 Xi exp(β2 Xi )]2 = 2 σ i=1

a) Tenemos que el modelo de regresi´on no lineal: Yi = β0 − β1 β2Xi el cual lo podemos rescribir como Yi − β0 = −β1 β2Xi ln(Yi − β0 ) = ln(−β1 ) + ln(β2 )Xi Yi∗ = β1∗ + β2∗ Xi Donde tenemos Yi∗ = ln(Yi − β0 ), β1∗ = ln(−β2 ) y β2∗ = ln(β2 ), entonces el modelo es intr´ınsicamente lineal si β0 es conocido y β0 ≤ Yi , ademas β1 ≤ −1 y β2 ≥ 1 b) Para este caso tenemos f (β0 , β1 , β2 , X) = β0 − β1 β2X entonces i = Yi − β0 + β1 β2Xi

149 y por tanto la funci´on de log-verosimilitud viene dada por: n
2 n 1 X l(β0 , β1 , β2 , X) = − ·ln(2π)−n·ln(σ)− 2 Yi − β0 + β1 β2Xi , 2 2σ i=1

para i = 1, 2, ..., n. Entonces las funciones de scores son: n
∂l(β0 , β1 , β2 , X) 1 X Yi − β0 + β1 β2Xi = 2 ∂β0 σ i=1 n
∂l(β0 , β1 , β2 , X) 1 X Yi − β0 + β1 β2Xi β2Xi =− 2 ∂β1 σ i=1 n
∂l(β0 , β1 , β2 , X) 1 X =− 2 Yi − β0 + β1 β2Xi β1 β2Xi ln(β2 ) ∂β2 σ i=1 n
2 n 1 X ∂l(β0 , β1 , β2 , X) =− + 3 Yi − β0 + β1 β2Xi ∂σ σ σ i=1

As´ı tenemos jβ0 β0 =

n σ2

jβ1 β0 = −

n 1 X Xi β σ 2 i=1 2

jβ1 β1

n 1 X 2Xi = 2 β σ i=1 2

jβ2 β0

n 1 X =− 2 β1 β2Xi ln(β2 ) σ i=1

jβ2 β1

n 2 X = 2 β1 β22Xi ln(β2 ) σ i=1

150

jβ2 β2

n      1 X = 2 β1 β2Xi −1 (Yi − β0 ) β2 ln(β2 )2 + 1 + β1 2β2Xi +1 ln(β2 )2 + β2Xi σ i=1

jσβ0 =

n
2 X Xi Y − β + β β i 0 1 2 σ 3 i=1

jσβ1

n
2 X Yi − β0 + β1 β2Xi β2Xi =− 3 σ i=1

jσβ2

n
2 X =− 3 Yi − β0 + β1 β2Xi β1 β2Xi ln(β2 ) σ i=1

y jσσ

n
2 n 2 X =− 2 + 4 Yi − β0 + β1 β2Xi σ σ i=1

Por lo tanto, los elementos de la matriz de informaci´on son: iβ1 β0

n 1 X Xi β =− 2 σ i=1 2

iβ2 β1

n 1 X Xi =− 2 β [1 + β1 ln(β2 )] σ i=1 2

n

iβ0 β0

X 1 = 2 (n), iβ1 β1 = β22Xi σ i=1

iβ2 β0

n 1 X =− 2 β1 β2Xi ln(β2 ) σ i=1

iσβ0 = 0 = iσβ1 , iσσ =

2n σ2

y i β2 β2 =

n 1 X [β1 β2Xi ln(β2 )]2 σ 2 i=1

10. Aplicando logaritmo natural al modelo tenemos ln(Ti ) = βo + β1 Xi + ln(εi ). Ahora tomando Ti∗ = ln(Ti ) y ε∗i = ln(εi ), entoces el modelo linealizado es: Ti∗ = βo + β1 Xi + ε∗i . Ahora El modelo a estudiar ser´a Ti∗ = βo + β1 Xi + ε∗i . Ajustando el modelo por MCO, quedar´ıa

151 reg|t|) (Intercept) 301.98 74.11 4.074 0.00114 ** X -25.35 7.89 -3.212 0.00626 ** ---

153 Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 42.17 on 14 degrees of freedom Multiple R-squared: 0.4243,Adjusted R-squared: 0.3832 F-statistic: 10.32 on 1 and 14 DF, p-value: 0.006264

Del anava parcial se observa que los par´ametros estimados βb0 y βb1 resultan significativos ya que no se acepta la hip´otesis nula H0 : βb0 = 0 y H0 : βb1 = 0 por que el p-valor de las prueba es menor que el nivel de significancia 0.05 reafirmando los resultados del an´alisis de varianza con un error est´andar de 74.11 y 7.89. Ademas Podemos concluir que el ajuste del modelo explica en un 42.43 % la variabilidad de los datos, y que es muy peque´ na la diferencia con el valor del coeficiente de determinaci´on ajustado. En base al coeficiente de determinaci´on obtenido, el ajuste es muy bueno. Asimismo, dicha bondad viene corroborada estad´sticamente con un p-valor de 0.006264. Luego |tc1 | = 3,212 Ahora se realizara la validaci´on de supuestos del modelo ajustado Normalidad Para verificar la normalidad para los errores del modelo ajustado usamos el test de Shapiro-Wilk teniendo como resultado: shapiro.test(reg$res) Shapiro-Wilk normality test data: reg$res W = 0.95328, p-value = 0.5433 Para el juego de hip´otesis: H0 : los residuales se distribuyen normal vs H1 : los residuales no se distribuyen normal a un nivel de significancia del 0.05 no se rechaza la hip´otesis nula ya que el p-valor = 0.0.5433 es mayor que 0.05. Por tanto los residuales del modelo siguen una distribuci´0n normal.

154 Independencia Para verificar la independencia para los errores del modelo ajustado usamos el test de Durbin-Watson de la libreria library(car), teniendo como resultado: durbinWatsonTest(reg) lag Autocorrelation D-W Statistic p-value 1 -0.3072194 2.426917 0.504 Alternative hypothesis: rho != 0 La hip´otesis nula en ´esta prueba es H0 :Los residuales son independientes. El p-valor de la prueba es de 0,504 > 0,05 por tanto no se rechaza la hip´otesis nula, es decir que se cumple el supuesto de que los errores son independientes. Homogeneidad de Varianzas bptest(reg) studentized Breusch-Pagan test data: reg BP = 0.10871, df = 1, p-value = 0.7416 La hip´otesis nula en ´esta prueba es son H0 :Los residuales tienen varianza constante. De acuerdo con los resultados de la prueba, el p-valor es 0.7416 el cual es mayor que 0.05 lo que lleva al no rechazo de la hip´o1tesis nula, es decir que el modelo no presenta problemas de heterocedasticidad. y por ultimo obtengo el modelo estimado para los datos originales. reg mod4 = glm(y~x,family = poisson());mod4 Call:

glm(formula = y ~ x, family = poisson())

Coefficients: (Intercept) 1.8893

x 0.6698

Degrees of Freedom: 8 Total (i.e. Null); Null Deviance: 18.42 Residual Deviance: 2.939 AIC: 41.05

7 Residual

> summary(mod4) Call: glm(formula = y ~ x, family = poisson()) Deviance Residuals: Min 1Q Median -0.8472 -0.2601 -0.2137

3Q 0.5214

Max 0.8788

Coefficients: Estimate Std. Error z value Pr(>|z|)

156

0.0 −0.5

resid(mod4)

0.5

Normalidad resi

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

norm quantiles

Figura 23: Gr´aficos de diagn´ostico

(Intercept) 1.8893 0.1421 13.294 < 2e-16 *** x 0.6698 0.1787 3.748 0.000178 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 18.4206 Residual deviance: 2.9387 AIC: 41.052

on 8 on 7

degrees of freedom degrees of freedom

Number of Fisher Scoring iterations: 4 Las estimaciones de los par´ametros del modelo log-lineal de Poisson para los datos son dadas en la tabla ??. En el Figura?? se evidencia que el modelo esta bien ajustado El modelo estimado viene dado por: µ b(x) = e1,8893+0,6698x

157 Par´ametro β0 β1 Desv´ıo

Estimaci´on 1.8893 0.6698 2.939

Cuadro 12: Estimativas para el modelo log-lineal de Poisson.

Note que el desv´ıo es de 2.939 con 8 grados de libertad, entonces el p − valor = 0,4123 > 0,05, es decir para este caso tendr´ıamos un ajuste de R2 = 1 − pchisq(2,939, 8) = 0,9381 por tanto el modelo tiene un buen ajuste. Ahora se realiza un an´alisis diagn´osticos para analizar si existen observaciones at´ıpicas o muy influyentes sobre las estimativas de los par´ametros del modelo, las cuales tambi´en pueden influenciar en las pruebas de validaci´on de los errores del modelo. An´ alisis diagnostico para el modelo: Un an´alisis de diagnostico graficamente y en R es: > influence.measures(mod4) Influence measures of glm(formula = y ~ x, family = poisson()) :

1 2 3 4 5 6 7 8 9

dfb.1_ -0.8597 -0.1889 -0.1481 0.0899 0.3325 0.6346 -0.3297 -0.0828 0.1919

dfb.x 0.9170 0.2015 0.0827 -0.0502 -0.1858 -0.3546 -0.7745 -0.1946 0.4508

dffit -1.0081 -0.2215 -0.1481 0.0899 0.3325 0.6346 -1.2022 -0.3020 0.6997

cov.r 0.887 1.792 1.499 1.544 1.253 0.763 0.824 1.835 1.400

cook.d 0.14601 0.01130 0.00508 0.00200 0.02582 0.07655 0.21121 0.02108 0.10653

hat inf 0.273 0.273 0.134 0.134 0.134 0.134 0.307 0.307 0.307

De lo anterior se concluye que no existen observaciones at´ıpicas o de alta influencia sobre las estimativas de los par´ametros del modelo pues

158

0.30

Cook's distance

DFFITS

0.0

0.25

−0.5

hatvalues(mod4)

0.10

9

0.20

1

0.00

0.15

−1.0

0.05

Cook's distance

0.15

0.20

0.5

7

2

4

6

8

2

4

6

8

2

4

Index

6

8

6

8

Index

1.4 Çovratio 1.2

0.0

DFBETAS para X1

0.0

0.8

−0.5

1.0

−0.5

DFBETAS para el intercepto

0.5

1.6

0.5

1.8

Obs. number

2

4

6

8

2

4

Index

6

8

2

Index

4 Index

Figura 24: Gr´aficos de diagn´ostico

las observaciones mas influyentes no afectan mucho la variabilidad del modelo. 4

a) Ajustar el modelo de regresi´on log´ıstico  log

pi 1 − pi

 = β0 + β1 Xi

Se tiene que en forma general el modelo de regresi´on log´ıstico para las variables explicativas X1 , X2 , · · · Xp−1 , viene dado por:

P (Yi = 1) =

1 1 + exp {−(β0 + β1 Xi1 + β2 Xi2 + · · · + βp−1 Xi(p−1) )}

donde β0 , β1 , · · · , βp−1 son los par´ametros del modelo. En este caso p=1, de modo que P (Yi = 1) =

1 1 + exp {−(β0 + β1 Xi )}

159 usando las funciones en R para el ajuste de un modelo log-lineal se tienen las siguientes estimativas > mod anova(mod) Analysis of Deviance Table Model: binomial, link: logit Response: Estado Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev

160 NULL Ingresos

19 18

1 0.033813

27.526 27.492

Donde se observa que el par´ametro β1 no es signi cativo en el modelo, por lo tanto no se rechaza la hip´otesis nula H0 : β1 = 0. c) Verifique un buen ajuste del modelo. Mediante la funci´on summary, podemos corroborar el ajuste del modelo, en efecto: > summary(mod) Call: glm(formula = Estado ~ Ingresos, family = binomial(link = logit), data = datos) Deviance Residuals: Min 1Q Median -1.355 -1.251 1.040

3Q 1.103

Max 1.106

Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.435e-01 5.462e-01 0.263 0.793 Ingresos 6.699e-07 3.672e-06 0.182 0.855 (Dispersion parameter for binomial family taken to be 1) Null deviance: 27.526 Residual deviance: 27.492 AIC: 31.492

on 19 on 18

degrees of freedom degrees of freedom

Number of Fisher Scoring iterations: 4 d) Ajuste ahora el modelo log´ıstico   pi = β0 + β1 Xi + β11 Xi2 log 1 − pi usando las funciones en R para el ajuste de un modelo log-lineal se tienen las siguientes estimativas

161 > mod1 Anova(mod1) Analysis of Deviance Table (Type II tests) Response: Estado LR Chisq Df Pr(>Chisq) Ingresos 0.40574 1 0.5241 I(Ingresos^2) 0.44142 1 0.5064 Donde se observa que el par´ametro β11 no es significativo en el modelo, por lo tanto no se rechaza la hip´otesis nula H0 : β11 = 0. f) Compare el AIC y BIC para los dos modelos, concluya.   pi = β0 + β1 Xi , el BIC es: Notemos que para el modelo log 1−p i > BIC(mod) [1] 33.4832 Mientras que para el modelo log datos por



pi 1−pi



= β0 +β1 Xi +β11 Xi2 viene

162 > BIC(mod1) [1] 36.03751 De modo que para comparar los modelos, tenemos: Modelo

AIC

BIC

pi log 1−p = β0 + β1 Xi i  pi log 1−p = β0 + β1 Xi + β11 Xi2 i

31.492

33.4832

33.05

36.03751







Como ambos criterios buscan el menor valor para un mejor ajuste del modelo, podemos notar que el modelo de regresi´on log´ıstico que mejor ajusta viene dado por la expresi´on:  log

pi 1 − pi

 = β0 + β1 Xi

dando a entender que el modelo de regresi´on lineal log´ıstico tiene un mejor ajuste a los datos.