Taller 3
Alejandro Pinto
2024-11-01
DATA Blocks
Para iniciar con el modelo de regresión generalizado, se inicia con la data blocks de la librería:
library(GLMsData)
data(blocks)
head(blocks)## Child Number Time Trial Shape Age
## 1 A 11 30.0 1 Cube 4.67
## 2 B 9 19.0 1 Cube 5.00
## 3 C 8 18.6 1 Cube 4.42
## 4 D 9 23.0 1 Cube 4.33
## 5 E 10 29.0 1 Cube 4.33
## 6 F 13 178.0 1 Cube 4.83Si revisamos las columnas y sus datos tenemos que:
table(blocks$Shape)##
## Cube Cylinder
## 50 50Observamos que la data tiene 6 variables:
Child: Identificador del niño, de A a Y. Number: Número de bloques que el niño logró apilar con éxito; un vector numérico. Time: Tiempo en segundos que los niños tardaron en apilar sus bloques; un vector numérico. Trial: Número de prueba en la que se recopilaron los datos; un factor con niveles 1 y 2. Shape: Forma de los bloques apilados; un factor con niveles Cubo y Cilindro. Age: Edad del niño en años completos; un vector numérico.
Para poder realizar el modelo lineal generalizado, necesitamos cambiar las variables categóricas a variables dummies:
blocks$Shape = ifelse(blocks$Shape == "Cube",1,0)
head(blocks)## Child Number Time Trial Shape Age
## 1 A 11 30.0 1 1 4.67
## 2 B 9 19.0 1 1 5.00
## 3 C 8 18.6 1 1 4.42
## 4 D 9 23.0 1 1 4.33
## 5 E 10 29.0 1 1 4.33
## 6 F 13 178.0 1 1 4.83En este caso, como la variable de interes es el tiempo de duran los niños en construir una torre, se realiza una prueba para determinar si se usa el modelo lineal generalizado gamma o inversa gausiana:
hist(blocks$Time)
modelblockg=glm(Time ~ Number + Trial + Shape + Age, data = blocks,family=Gamma(link="log"))
modelblockig=glm(Time ~ Number + Trial + Shape + Age, data = blocks,family=inverse.gaussian(link="log"))
log(logLik(modelblockig)/logLik(modelblockg))## 'log Lik.' -0.01018081 (df=6)Si el logaritmo del cociente de las versomilitudes del modelo son mayores a cero. Use un modelo inversa Gaussiana. En caso que este valor sea menor a cero un modelo gamma.
En este caso se observa que el valor es negativo, luego se hace un modelo gamma.
A. Seleccion del mejor subconjunto de variables predictoras:
Para esto recurrimos a la prueba de Akaike (AIC):
step(modelblockg)## Start: AIC=797.99
## Time ~ Number + Trial + Shape + Age
##
## Df Deviance AIC
## - Trial 1 27.581 796.01
## - Shape 1 28.166 797.39
## <none> 27.574 797.99
## - Age 1 32.579 807.77
## - Number 1 44.888 836.74
##
## Step: AIC=796.02
## Time ~ Number + Shape + Age
##
## Df Deviance AIC
## - Shape 1 28.175 795.43
## <none> 27.581 796.02
## - Age 1 32.582 805.90
## - Number 1 44.905 835.17
##
## Step: AIC=796.25
## Time ~ Number + Age##
## Call: glm(formula = Time ~ Number + Age, family = Gamma(link = "log"),
## data = blocks)
##
## Coefficients:
## (Intercept) Number Age
## 2.9451 0.2219 -0.2944
##
## Degrees of Freedom: 99 Total (i.e. Null); 97 Residual
## Null Deviance: 51.05
## Residual Deviance: 28.18 AIC: 796.2Como se puede observar, las variables Age y Number explican bien el el tiempo de construccion de una torre, con un valor AIC de 796.25.
Se eliminaron los predictores Trial y Shape, asi que se usarà este modelo:
bestmodelblockg =glm(formula = Time ~ Number + Age, family = Gamma(link = "log"),
data = blocks)para evaluar la calidad del ajuste al modelo
library(DescTools)
PseudoR2(bestmodelblockg, "Nagelkerke")## Nagelkerke
## 0.4685197Existe una buena explicacion del tiempo empleado de los niños para construir torres con bloques cùbicos y cilindricos, por el modelo de regresion gamma. Ahora del summary:
summary(bestmodelblockg)##
## Call:
## glm(formula = Time ~ Number + Age, family = Gamma(link = "log"),
## data = blocks)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.94508 0.37460 7.862 5.26e-12 ***
## Number 0.22193 0.03111 7.134 1.76e-10 ***
## Age -0.29443 0.09315 -3.161 0.0021 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.4303746)
##
## Null deviance: 51.050 on 99 degrees of freedom
## Residual deviance: 28.175 on 97 degrees of freedom
## AIC: 796.25
##
## Number of Fisher Scoring iterations: 6Podemos realizar la prueba Deviance:
1-pchisq(51.050-28.175,2)## [1] 1.078343e-05mayor a 0.05, no rechace ho, existe evidencia estadistica de que el modelo de regresion gamma no se ajusta a los dats.
Observandoo los coeficientes estimados: son menor a 0.05, rechace ho, existe evidencia estadistica de que el numero de bloques la edad del niño influye en el tiempo en el que deben cosntruir una torre.
coef(bestmodelblockg)## (Intercept) Number Age
## 2.9450795 0.2219259 -0.2944274Para la variable Number:
exp(0.2219259)## [1] 1.248479Por cada bloque adicional, el tiempo que tarda el niño en construir la torre aumenta 1.24 segundos, en promedio. Para la variable Age:
1/exp(-0.2944274)## [1] 1.342358Por cada semana edad adicional del niño, el tiempo que tarda el niño en costruir la torre disminuye 1.34 segundos, en promedio.
B. Prediccioncon number = 10, trail = 2, Shape = Cube and Age =6:
testblock = data.frame(Trial = 2, Number=10, Shape=1, Age=5)
predict(bestmodelblockg,testblock, type="response")## 1
## 40.1331Un niño con 5 años se espera que tenga arpoximadamente 40.1331 segundos en construir una torre con bloques, independientemente si son cubicos o no.
B. Datos atipicos e influyentes:
En el caso del modelo lineal generalizado, podemos revisar los datos atipicos residuales de Deviance y Pearson y comparamos sus resultados. Si graficamos los residuos:
par(mfrow=c(1,2))
plot(abs(residuals(bestmodelblockg)), main = "Prueba Deviance")
abline(h=2,col="red")
plot(abs(residuals(bestmodelblockg,type="pearson")), main = "Prueba Pearson")
abline(h=2,col="red")
No se observan datos atipicos en la primera prubea, pero se registran 4
datos atipicos en la segunda prueba.
atipicos=data.frame(deviance=abs(residuals(bestmodelblockg)), pearson=abs(residuals(bestmodelblockg,type="pearson")))
atipicos[atipicos$deviance>2 & atipicos$pearson>2,]## [1] deviance pearson
## <0 rows> (or 0-length row.names)La tabla anerior nos muestra los datos mayores a 2 para la prueba Deviance y la prueba Pearson. Es decir, hay 0 datos que se consideran datos atìpicos.
Pasando a los datos influyentes, podemos revisar los StudResd, Hat-values y la distancia de Cook:
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:DescTools':
##
## Recodepar(mfrow=c(1,1))
influencePlot(bestmodelblockg)
## StudRes Hat CookD
## 6 1.798213385 0.08802096 1.118017e-01
## 47 3.330975397 0.06905116 3.839142e-01
## 56 0.006874005 0.08802096 1.038929e-06
## 76 4.008848835 0.01829075 1.704523e-01Se observan 3 posibles datos influyentes: 1, 2 y 7. Si analizamos La distancia de Cook, ninguno pasa la prueba. Por ende, no se observan datos influyentes.
DATA Leukwbc
Para iniciar con el modelo de regresión generalizado, se inicia con la data Leeukwbc de la librería:
library(GLMsData)
data("leukwbc")
dim(leukwbc)## [1] 33 3head(leukwbc)## WBC Time AG
## 1 2300 65 1
## 2 750 156 1
## 3 4300 100 1
## 4 2600 134 1
## 5 6000 16 1
## 6 10500 108 1Si revisamos las columnas y sus datos tenemos que:
hist(leukwbc$Time)
Observamos que la data con 33 observaciones sobre las siguientes 3 variables:
WBC: Conteo de glóbulos blancos; un vector numérico. Time: Tiempo hasta la muerte en semanas; un vector numérico. AG: Variable morfológica, el factor ag; un vector numérico donde 1 significa ag-positivo y 2 significa ag-negativo.
En este caso, como la variable de interes es el tiempo hasta la muerte y conteo de glóbulos blancos en dos grupos de pacientes con leucemia, se realiza una prueba para determinar si se usa el modelo lineal generalizado gamma o inversa gausiana:
modeleukg=glm(Time ~ WBC + AG, data = leukwbc,family=Gamma(link="log"))
modeleukig=glm(Time ~ WBC + AG, data = leukwbc,family=inverse.gaussian(link="log"))## Warning: glm.fit: algorithm did not convergelog(logLik(modeleukig)/logLik(modeleukg))## 'log Lik.' 0.03151353 (df=4)Si el logaritmo del cociente de las versomilitudes del modelo son mayores a cero. Use un modelo inversa Gaussiana. En caso que este valor sea menor a cero un modelo gamma.
En este caso se observa que el valor es negativo, luego se hace un modelo inversa gausiana.
A. Seleccion del mejor subconjunto de variables predictoras:
Para esto recurrimos a la prueba de Akaike (AIC):
step(modeleukg)## Start: AIC=305.18
## Time ~ WBC + AG
##
## Df Deviance AIC
## - WBC 1 46.198 304.99
## <none> 44.292 305.18
## - AG 1 53.085 311.52
##
## Step: AIC=304.85
## Time ~ AG
##
## Df Deviance AIC
## <none> 46.198 304.85
## - AG 1 58.138 314.67##
## Call: glm(formula = Time ~ AG, family = Gamma(link = "log"), data = leukwbc)
##
## Coefficients:
## (Intercept) AG
## 5.382 -1.248
##
## Degrees of Freedom: 32 Total (i.e. Null); 31 Residual
## Null Deviance: 58.14
## Residual Deviance: 46.2 AIC: 304.9Como se puede observar, la variable AG explica bien el el tiempo de muerte de pacientes con leucemia, con un valor AIC de 304.85.
Se elimin`aron los predictores Trial y Shapeò el predictor WBC, asi que se usarà este modelo:
bestmodeleukig =glm(formula = Time ~ AG, family = Gamma(link = "log"), data = leukwbc)para evaluar la calidad del ajuste al modelo
library(DescTools)
PseudoR2(bestmodeleukig, "Nagelkerke")## Nagelkerke
## 0.2454965Existe una aceptable explicacion del tiempo de muerte de pacientes con leucemia, por el modelo de regresion inversa gausiana. Ahora del summary:
summary(bestmodeleukig)##
## Call:
## glm(formula = Time ~ AG, family = Gamma(link = "log"), data = leukwbc)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3825 0.5486 9.812 5.04e-11 ***
## AG -1.2478 0.3502 -3.564 0.00121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.010594)
##
## Null deviance: 58.138 on 32 degrees of freedom
## Residual deviance: 46.198 on 31 degrees of freedom
## AIC: 304.85
##
## Number of Fisher Scoring iterations: 6Podemos realizar la prueba Deviance:
1-pchisq(58.138-46.198,1)## [1] 0.0005494148menor a 0.05, rechace ho, existe evidencia estadistica de que el modelo de regresion inversa gausiana se ajusta a los dats.
Observandoo los coeficientes estimados:menor a0.05, rechace ho, existe evidencia estadistica de que la variable morfològica AG influye en el tiempo de muerte de pacientes con leucemia.
coef(bestmodeleukig)## (Intercept) AG
## 5.382498 -1.2478021/exp(coef(bestmodeleukig))## (Intercept) AG
## 0.004596325 3.482680877Por cada valor AG, el tiempo de muerte de un paciente con leucemia disminuye 3.48, en promedio.
B. Prediccion con WBG = 1000, AG = 1:
testleuk = data.frame(WBC = 1000, AG=1)
predict(bestmodeleukig,testleuk, type="response")## 1
## 62.47059Un paciente con leucemia, Ag positivo y 1000 celulas blancas, se espera que el tiempo de muerte sea de 62.47 semanas.
B. Datos atipicos e influyentes:
En el caso del modelo lineal generalizado, podemos revisar los datos atipicos residuales de Deviance y Pearson y comparamos sus resultados. Si graficamos los residuos:
par(mfrow=c(1,2))
plot(abs(residuals(bestmodeleukig)), main = "Prueba Deviance")
abline(h=2,col="red")
plot(abs(residuals(bestmodeleukig,type="pearson")), main = "Prueba Pearson")
abline(h=2,col="red")
Se observan dos datos atipicos
residuos=data.frame(deviance=abs(residuals(bestmodeleukig)),
pearson=abs(residuals(bestmodeleukig,
type="pearson")))
residuos[residuos$deviance>2& residuos$pearson>2,]## [1] deviance pearson
## <0 rows> (or 0-length row.names)La tabla anerior nos muestra los datos mayores a 2 para la prueba Deviance y la prueba Pearson. Es decir, hay o datos que se consideran datos atìpicos.
Pasando a los datos influyentes, podemos revisar los StudResd, Hat-values y la distancia de Cook:
library(car)
par(mfrow=c(1,1))
influencePlot(bestmodeleukig)
## StudRes Hat CookD
## 14 -2.198060 0.05882353 0.03181161
## 15 -2.198060 0.05882353 0.03181161
## 18 1.241887 0.06250000 0.15841688
## 19 1.472112 0.06250000 0.24219034Se observan 4 posibles datos influyentes: 14, 15, 18 y 19. Si analizamos La distancia de Cook, ninguno pasa la prueba. Por ende, no se observan datos influyentes.
DATA Affairs
Para iniciar con el modelo de regresión generalizado, se inicia con la data Affairs de la librería:
library(AER)## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Loading required package: survivaldata("Affairs")
dim(Affairs)## [1] 601 9head(Affairs)## affairs gender age yearsmarried children religiousness education occupation
## 4 0 male 37 10.00 no 3 18 7
## 5 0 female 27 4.00 no 4 14 6
## 11 0 female 32 15.00 yes 1 12 1
## 16 0 male 57 15.00 yes 5 18 6
## 23 0 male 22 0.75 no 2 17 6
## 29 0 female 32 1.50 no 2 17 5
## rating
## 4 4
## 5 4
## 11 4
## 16 5
## 23 3
## 29 5Si revisamos las columnas y sus datos tenemos que:
table(Affairs$children)##
## no yes
## 171 430Affairs$gender = ifelse(Affairs$gender=="male",1,0)
Affairs$children = ifelse(Affairs$children=="yes",1,0)Observamos que la data con 33 observaciones sobre las siguientes 9 variables:
affairs: Frecuencia de participación en relaciones sexuales extramaritales en el último año (variable numérica): 0 = ninguna, 1 = una vez, 2 = dos veces, 3 = tres veces, 7 = de 4 a 10 veces, 12 = mensual, semanal o diario. gender: Género; variable categórica. age: Edad en años (variable numérica): 17.5 = menor de 20, 22 = 20–24, 27 = 25–29, 32 = 30–34, 37 = 35–39, 42 = 40–44, 47 = 45–49, 52 = 50–54, 57 = 55 o más. yearsmarried: Años de matrimonio (variable numérica): 0.125 = 3 meses o menos, 0.417 = 4–6 meses, 0.75 = de 6 meses a 1 año, 1.5 = de 1 a 2 años, 4 = de 3 a 5 años, 7 = de 6 a 8 años, 10 = de 9 a 11 años, 15 = 12 años o más. children: Indica si hay hijos en el matrimonio; variable categórica. religiousness: Nivel de religiosidad (variable numérica): 1 = anti-religioso, 2 = nada religioso, 3 = algo religioso, 4 = moderadamente religioso, 5 = muy religioso. education: Nivel educativo (variable numérica): 9 = escuela primaria, 12 = graduado de secundaria, 14 = algo de universidad, 16 = graduado de universidad, 17 = estudios de posgrado, 18 = maestría, 20 = doctorado, medicina u otro título avanzado. occupation: Ocupación según la clasificación de Hollingshead (con numeración inversa); variable numérica. rating: Autoevaluación de la felicidad en el matrimonio (variable numérica): 1 = muy infeliz, 2 = algo infeliz, 3 = promedio, 4 = más feliz que el promedio, 5 = muy feliz.
En este caso, como la variable de interes es la cantidad de infidelidades, se realiza una prueba de Cameron para determinar si sobredispersion o no:
modelaffair=glm(affairs~.,family=poisson,data=Affairs)
barplot(Affairs$affairs)
dispersiontest(modelaffair,trafo=1)##
## Overdispersion test
##
## data: modelaffair
## z = 5.2595, p-value = 7.224e-08
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 5.873703menor a 0.05. REchace Ho. Existe evidencia estadistica de que hay sobredispersion. En este caso se realiza el modelo con la distribuciòn binomial negativa
library(MASS)
modelaffairnb=glm.nb(affairs~.,data=Affairs)
summary(modelaffairnb) ##
## Call:
## glm.nb(formula = affairs ~ ., data = Affairs, init.theta = 0.1427050721,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3381919 1.0425812 2.243 0.024916 *
## gender -0.0186320 0.2776900 -0.067 0.946505
## age 0.0002843 0.0207056 0.014 0.989046
## yearsmarried 0.0803867 0.0380018 2.115 0.034401 *
## children 0.1161738 0.3288787 0.353 0.723907
## religiousness -0.4257717 0.1042682 -4.083 4.44e-05 ***
## education -0.0260332 0.0591538 -0.440 0.659869
## occupation 0.0807710 0.0828678 0.975 0.329711
## rating -0.4152284 0.1087853 -3.817 0.000135 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1427) family taken to be 1)
##
## Null deviance: 391.17 on 600 degrees of freedom
## Residual deviance: 339.61 on 592 degrees of freedom
## AIC: 1476.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1427
## Std. Err.: 0.0160
##
## 2 x log-likelihood: -1456.2010A. Seleccion del mejor subconjunto de variables predictoras:
Para esto recurrimos a la prueba de Akaike (AIC):
step(modelaffairnb)## Start: AIC=1474.2
## affairs ~ gender + age + yearsmarried + children + religiousness +
## education + occupation + rating
##
## Df Deviance AIC
## - age 1 339.61 1472.2
## - gender 1 339.62 1472.2
## - children 1 339.75 1472.3
## - education 1 339.79 1472.4
## - occupation 1 340.52 1473.1
## <none> 339.61 1474.2
## - yearsmarried 1 343.91 1476.5
## - rating 1 353.45 1486.0
## - religiousness 1 354.75 1487.3
##
## Step: AIC=1472.2
## affairs ~ gender + yearsmarried + children + religiousness +
## education + occupation + rating
##
## Df Deviance AIC
## - gender 1 339.63 1470.2
## - children 1 339.76 1470.3
## - education 1 339.80 1470.4
## - occupation 1 340.53 1471.1
## <none> 339.62 1472.2
## - yearsmarried 1 348.65 1479.2
## - rating 1 353.85 1484.4
## - religiousness 1 355.15 1485.7
##
## Step: AIC=1470.21
## affairs ~ yearsmarried + children + religiousness + education +
## occupation + rating
##
## Df Deviance AIC
## - children 1 339.75 1468.3
## - education 1 339.81 1468.4
## - occupation 1 340.56 1469.2
## <none> 339.61 1470.2
## - yearsmarried 1 348.70 1477.3
## - rating 1 353.89 1482.5
## - religiousness 1 355.56 1484.2
##
## Step: AIC=1468.35
## affairs ~ yearsmarried + religiousness + education + occupation +
## rating
##
## Df Deviance AIC
## - education 1 339.75 1466.5
## - occupation 1 340.46 1467.2
## <none> 339.59 1468.3
## - yearsmarried 1 350.85 1477.6
## - rating 1 355.19 1482.0
## - religiousness 1 355.43 1482.2
##
## Step: AIC=1466.51
## affairs ~ yearsmarried + religiousness + occupation + rating
##
## Df Deviance AIC
## - occupation 1 340.28 1465.2
## <none> 339.57 1466.5
## - yearsmarried 1 351.02 1476.0
## - religiousness 1 355.42 1480.4
## - rating 1 357.10 1482.0
##
## Step: AIC=1465.22
## affairs ~ yearsmarried + religiousness + rating
##
## Df Deviance AIC
## <none> 339.33 1465.2
## - yearsmarried 1 352.35 1476.2
## - religiousness 1 356.10 1480.0
## - rating 1 356.21 1480.1##
## Call: glm.nb(formula = affairs ~ yearsmarried + religiousness + rating,
## data = Affairs, init.theta = 0.141882111, link = log)
##
## Coefficients:
## (Intercept) yearsmarried religiousness rating
## 2.36869 0.08531 -0.43515 -0.42355
##
## Degrees of Freedom: 600 Total (i.e. Null); 597 Residual
## Null Deviance: 389.6
## Residual Deviance: 339.3 AIC: 1467Como se puede observar, las variables yearsmarried, religiousness y rating explican bien la cantidad de infidelidades en 1969, con un valor AIC de 304.85.
Se eliminò el predictor age, gender, children, education, occupation, asi que se usarà este modelo:
bestmodelaffairnb=glm.nb(formula = affairs ~ yearsmarried + religiousness + rating,
data = Affairs, init.theta = 0.141882111, link = log)para evaluar la calidad del ajuste al modelo
library(DescTools)
PseudoR2(bestmodelaffairnb, "Nagelkerke")## Nagelkerke
## 0.07881018Existe una mala explicacion de la cantidad de infidelidades, por el modelo de regresion binomial negativa. . Ahora del summary:
summary(bestmodelaffairnb)##
## Call:
## glm.nb(formula = affairs ~ yearsmarried + religiousness + rating,
## data = Affairs, init.theta = 0.1418818828, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.36871 0.55987 4.231 2.33e-05 ***
## yearsmarried 0.08531 0.02242 3.806 0.000141 ***
## religiousness -0.43516 0.10408 -4.181 2.90e-05 ***
## rating -0.42355 0.10748 -3.941 8.13e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1419) family taken to be 1)
##
## Null deviance: 389.62 on 600 degrees of freedom
## Residual deviance: 339.33 on 597 degrees of freedom
## AIC: 1467.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1419
## Std. Err.: 0.0158
##
## 2 x log-likelihood: -1457.2180Podemos realizar la prueba Deviance:
1-pchisq(389.62-339.33,3)## [1] 6.930079e-11menor a 0.05, rechace ho, existe evidencia estadistica de que el modelo de regresion binomial inversa se ajusta a los datos.
Observandoo los coeficientes estimados:menor a0.05, rechace ho, existe evidencia estadistica de que las variables yearsmarried, religiousness y rating influyen en el numero de infidelidades.
Interpretacion de los coeficientes estimados del modelo
coef(bestmodelaffairnb)## (Intercept) yearsmarried religiousness rating
## 2.36870641 0.08530847 -0.43516099 -0.42355410exp(0.08530847)## [1] 1.089053Si el numero de años casados aumenta, el numero de infidelidades aumenta en promedio 1.089.
1/exp(-0.43516099)## [1] 1.545212Si la religiosidad disminuye, el numero de infidelidades aumenta en promedio 0.435.
1/exp(-0.42355410)## [1] 1.52738Si el matrimonio va mal, el numero de infidelidades aumenta en promedio 1.527.
B. Datos atipicos e influyentes:
En el caso del modelo lineal generalizado, podemos revisar los datos atipicos residuales de Deviance y Pearson y comparamos sus resultados. Si graficamos los residuos:
par(mfrow=c(1,2))
plot(abs(residuals(bestmodelaffairnb)), main = "Prueba Deviance")
abline(h=2,col="red")
plot(abs(residuals(bestmodelaffairnb,type="pearson")), main = "Prueba Pearson")
abline(h=2,col="red")
Se observa una amplia data de datos atipìcos
residuos=data.frame(deviance=abs(residuals(bestmodelaffairnb)),
pearson=abs(residuals(bestmodelaffairnb,
type="pearson")))
residuos[residuos$deviance>2& residuos$pearson>2,]## deviance pearson
## 1622 2.445133 9.185573
## 1674 2.053562 6.715776La tabla anerior nos muestra los datos mayores a 2 para la prueba Deviance y la prueba Pearson. Es decir, hay 2 datos que se consideran datos atìpicos.
Pasando a los datos influyentes, podemos revisar los StudResd, Hat-values y la distancia de Cook:
library(car)
par(mfrow=c(1,1))
influencePlot(bestmodelaffairnb)
## StudRes Hat CookD
## 751 -1.2364689 0.019173523 0.0006734711
## 1622 3.3610697 0.004006078 0.0851840990
## 1669 -0.4795884 0.022144568 0.0003780502
## 1674 2.8263272 0.006039737 0.0689306672Se observan 4 posibles datos influyentes: 751, 1622, 1669 y 1674. Si analizamos La distancia de Cook, ninguno pasa la prueba. Por ende, no se observan datos influyentes.
C. Prediccion con gender = 1, Age = 37, Yearsmarried = 8, Children = 0, religiousness = 3, education = 188, rating = 5 y occupation = 8:
testaffair = data.frame(yearsmarried = 8, religiousness=3, rating = 5)
predict(bestmodelaffairnb,testaffair, type="response")## 1
## 0.6892996El promedio esperado de infidelidad de un hombre casado 8 años, poco religioso y muy feliz en el matrimnio, es de 0.68.