FINAL REGRESION’
Ricardo Enciso
2023-08-26
#PUNTO 1
#La base de datos danishlc contiene el n´umero de incidentes de c´ancer de pulm´on de 1968 a 1971 en cuatro ciudades danesas, con el registrado por grupo de edad y el n´umero de sujetos en cada grupo de edad. #(a) Ajuste el modelo de regresi´on que permita estimar el n´umero de sujetos con c´ancer de pulm´on. Muestre las meidas de bondad de ajuste. #(b) Estime el n´umero medio hombres con c´ancer de pulm´on en la ciudad de V ejle en un grupo de 1000 hombres de 40 a˜nos. #(c) Muestre los residuos de pearson y la deviance. ¿Existe alg´un dato at´ıpico?
remove(list=ls(all=TRUE))
library(car)## Warning: package 'car' was built under R version 4.2.3## Loading required package: carDatalibrary(GLMsData)
library(carData)
library(betareg)## Warning: package 'betareg' was built under R version 4.2.3library(AER)## Warning: package 'AER' was built under R version 4.2.3## Loading required package: lmtest## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.2.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.2.3## Loading required package: survivaldata("danishlc")
danishlctable (danishlc$Age)##
## >74 40-54 55-59 60-64 65-69 70-74
## 4 4 4 4 4 4table (danishlc$City)##
## Fredericia Horsens Kolding Vejle
## 6 6 6 6danishlc$Age40a54 = ifelse(danishlc$Age=="40-54", 1,0)
danishlc$Age55a59 = ifelse(danishlc$Age=="55-59", 1,0)
danishlc$Age60a64 = ifelse(danishlc$Age=="60-64", 1,0)
danishlc$Age65a69 = ifelse(danishlc$Age=="65-69", 1,0)
danishlc$Age70a74 = ifelse(danishlc$Age=="70-74", 1,0)danishlc$Fredericia = ifelse(danishlc$City=="Fredericia", 1,0)
danishlc$Horsens = ifelse(danishlc$City=="Horsens", 1,0)
danishlc$Kolding = ifelse(danishlc$City=="Kolding", 1,0)
danishlcdanishlc = danishlc[, -(3:4)]
danishlcmod1 = glm(Cases~.,data=danishlc , family = poisson)
dispersiontest(mod1, trafo = 1)##
## Overdispersion test
##
## data: mod1
## z = -0.87126, p-value = 0.8082
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## -0.2092929#No rechace Ho. existe evidencia estadística de que no hay sobredispersión por lo tanto se utilizará el modelo de regresión poisson
#Selección “mejor modelo”
step(mod1)## Start: AIC=136.83
## Cases ~ Pop + Age40a54 + Age55a59 + Age60a64 + Age65a69 + Age70a74 +
## Fredericia + Horsens + Kolding
##
## Df Deviance AIC
## - Kolding 1 20.496 134.88
## - Age55a59 1 20.603 134.99
## - Pop 1 20.673 135.06
## - Age40a54 1 20.705 135.09
## - Horsens 1 21.122 135.51
## - Age70a74 1 21.186 135.57
## - Fredericia 1 22.023 136.41
## <none> 20.445 136.83
## - Age60a64 1 22.576 136.96
## - Age65a69 1 23.122 137.51
##
## Step: AIC=134.88
## Cases ~ Pop + Age40a54 + Age55a59 + Age60a64 + Age65a69 + Age70a74 +
## Fredericia + Horsens
##
## Df Deviance AIC
## - Age55a59 1 20.615 133.00
## - Pop 1 20.673 133.06
## - Age40a54 1 20.706 133.09
## - Horsens 1 21.267 133.66
## - Age70a74 1 21.315 133.70
## - Fredericia 1 22.227 134.62
## <none> 20.496 134.88
## - Age60a64 1 22.585 134.97
## - Age65a69 1 23.153 135.54
##
## Step: AIC=133
## Cases ~ Pop + Age40a54 + Age60a64 + Age65a69 + Age70a74 + Fredericia +
## Horsens
##
## Df Deviance AIC
## - Pop 1 20.689 131.08
## - Age40a54 1 20.710 131.10
## - Horsens 1 21.293 131.68
## - Age70a74 1 21.334 131.72
## - Fredericia 1 22.458 132.85
## <none> 20.615 133.00
## - Age60a64 1 23.084 133.47
## - Age65a69 1 23.331 133.72
##
## Step: AIC=131.08
## Cases ~ Age40a54 + Age60a64 + Age65a69 + Age70a74 + Fredericia +
## Horsens
##
## Df Deviance AIC
## - Age40a54 1 20.735 129.12
## - Horsens 1 21.293 129.68
## - Age70a74 1 22.058 130.45
## - Fredericia 1 22.678 131.07
## <none> 20.689 131.08
## - Age60a64 1 23.107 131.50
## - Age65a69 1 23.947 132.34
##
## Step: AIC=129.12
## Cases ~ Age60a64 + Age65a69 + Age70a74 + Fredericia + Horsens
##
## Df Deviance AIC
## - Horsens 1 21.340 127.73
## - Age70a74 1 22.097 128.49
## - Fredericia 1 22.724 129.11
## <none> 20.735 129.12
## - Age60a64 1 23.226 129.61
## - Age65a69 1 24.141 130.53
##
## Step: AIC=127.73
## Cases ~ Age60a64 + Age65a69 + Age70a74 + Fredericia
##
## Df Deviance AIC
## - Age70a74 1 22.701 127.09
## - Fredericia 1 22.819 127.21
## <none> 21.340 127.73
## - Age60a64 1 23.830 128.22
## - Age65a69 1 24.745 129.13
##
## Step: AIC=127.09
## Cases ~ Age60a64 + Age65a69 + Fredericia
##
## Df Deviance AIC
## - Fredericia 1 24.180 126.57
## - Age60a64 1 24.429 126.82
## <none> 22.701 127.09
## - Age65a69 1 25.232 127.62
##
## Step: AIC=126.57
## Cases ~ Age60a64 + Age65a69
##
## Df Deviance AIC
## - Age60a64 1 25.908 126.30
## <none> 24.180 126.57
## - Age65a69 1 26.711 127.10
##
## Step: AIC=126.3
## Cases ~ Age65a69
##
## Df Deviance AIC
## - Age65a69 1 27.704 126.09
## <none> 25.908 126.30
##
## Step: AIC=126.09
## Cases ~ 1##
## Call: glm(formula = Cases ~ 1, family = poisson, data = danishlc)
##
## Coefficients:
## (Intercept)
## 2.234
##
## Degrees of Freedom: 23 Total (i.e. Null); 23 Residual
## Null Deviance: 27.7
## Residual Deviance: 27.7 AIC: 126.1es de resaltar que la selección de modelo por pasos resulto sin predictores
modfinal1_danish=glm(Cases ~ 1, family = poisson, data = danishlc)
round(exp(modfinal1_danish$coefficients),2)## (Intercept)
## 9.33summary(modfinal1_danish)##
## Call:
## glm(formula = Cases ~ 1, family = poisson, data = danishlc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9163 -0.7994 0.2157 0.5304 1.7031
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.23359 0.06682 33.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 27.704 on 23 degrees of freedom
## Residual deviance: 27.704 on 23 degrees of freedom
## AIC: 126.09
##
## Number of Fisher Scoring iterations: 41-pchisq(27.704-27.704, 23-23)## [1] 1library(DescTools)## Warning: package 'DescTools' was built under R version 4.2.3##
## Attaching package: 'DescTools'## The following object is masked from 'package:car':
##
## RecodePseudoR2(modfinal1_danish,"Nagelkerke")## Nagelkerke
## 0par(mfrow=c(1,2))
plot(abs(residuals(modfinal1_danish)))
abline(h=2,col="red")
plot(abs(residuals(modfinal1_danish, type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(modfinal1_danish)),abs(residuals(modfinal1_danish,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]influence.measures(modfinal1_danish)## Influence measures of
## glm(formula = Cases ~ 1, family = poisson, data = danishlc) :
##
## dfb.1_ dffit cov.r cook.d hat inf
## 1 0.1012 0.1012 1.079 0.01350 0.0417
## 2 0.1012 0.1012 1.079 0.01350 0.0417
## 3 0.1012 0.1012 1.079 0.01350 0.0417
## 4 0.0410 0.0410 1.089 0.00216 0.0417
## 5 0.1012 0.1012 1.079 0.01350 0.0417
## 6 0.0410 0.0410 1.089 0.00216 0.0417
## 7 0.2203 0.2203 1.038 0.06535 0.0417
## 8 -0.2277 -0.2277 1.035 0.05401 0.0417
## 9 0.3425 0.3425 0.972 0.15609 0.0417
## 10 0.0410 0.0410 1.089 0.00216 0.0417
## 11 0.1607 0.1607 1.062 0.03457 0.0417
## 12 -0.6714 -0.6714 0.741 0.26141 0.0417 *
## 13 -0.4051 -0.4051 0.931 0.13827 0.0417
## 14 -0.0853 -0.0853 1.083 0.00864 0.0417
## 15 -0.1536 -0.1536 1.065 0.02647 0.0417
## 16 0.1012 0.1012 1.079 0.01350 0.0417
## 17 -0.0208 -0.0208 1.090 0.00054 0.0417
## 18 0.1607 0.1607 1.062 0.03457 0.0417
## 19 -0.3101 -0.3101 0.991 0.09128 0.0417
## 20 -0.1536 -0.1536 1.065 0.02647 0.0417
## 21 0.0410 0.0410 1.089 0.00216 0.0417
## 22 0.2806 0.2806 1.008 0.10586 0.0417
## 23 -0.0853 -0.0853 1.083 0.00864 0.0417
## 24 -0.1536 -0.1536 1.065 0.02647 0.0417influencePlot(modfinal1_danish)
#Existe una mala explicación del número de incidentes de cáncer de
pulmón por el modelo de regresión poisson para esta base de datos. #La
observación 12 puede considerarse atípica #Ala luz de todos los
criterios no se observan datos que se puedan considerar influyentes
pred=data.frame(Pop = 1000,Age40a54= 1, City= 0 )
predict( modfinal1_danish,pred,type="response")## 1
## 9.333333#El número medio de hombres con edad de 40 años en un grupo de 1000 con cáncer de pulmón en la ciudad de V ejle es de 9.33
#PUNTO 2 #2. En la base de datos trout se busca determinar el n´umero de huevos de trucha muertos por una toxina aplicada en los animales. #(a) Ajuste el ”mejor” modelo de regresi´on con para estimar el n´umero de huevos da˜nados como variable de int´eres. Muestre las medidas de bondad de ajuste. #(b) Realice una estimaci´on con Con = 90, When = later y Number = 234. #(c) Existe datos at´ıpicos o influyentes?. Muestre los residuales y las medidas de influencia.
remove(list=ls(all=TRUE))
library(car)
library(GLMsData)
library(carData)
library(betareg)
library(AER)
library(MASS)## Warning: package 'MASS' was built under R version 4.2.3data("trout")
head(trout)trout$ahora=ifelse(trout$When=="Now",1,0)
trout$despues=ifelse(trout$When=="Later",1,0)
trouttrout=trout[,-(2)]
troutmod2=glm(Dead~., family = poisson,data = trout)
mod2##
## Call: glm(formula = Dead ~ ., family = poisson, data = trout)
##
## Coefficients:
## (Intercept) Conc Number ahora despues
## 1.0715035 0.0007977 0.0089936 -0.0262859 NA
##
## Degrees of Freedom: 47 Total (i.e. Null); 44 Residual
## Null Deviance: 1538
## Residual Deviance: 383.9 AIC: 597.2dispersiontest(mod2,trafo=1)##
## Overdispersion test
##
## data: mod2
## z = 2.9176, p-value = 0.001764
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 7.807809#Existe evidencia estadística de que hay sobre dispersión en el modelo de regresión poisson, por lo tanto se utilizará un modelo binomial negativo
mod2.1=glm.nb(Dead~.,data = trout)#selección del “mejor modelo”
step(mod2.1)## Start: AIC=368.18
## Dead ~ Conc + Number + ahora + despues
##
##
## Step: AIC=368.18
## Dead ~ Conc + Number + ahora
##
## Df Deviance AIC
## - ahora 1 54.141 366.30
## <none> 54.020 368.18
## - Number 1 59.600 371.76
## - Conc 1 119.993 432.15
##
## Step: AIC=366.3
## Dead ~ Conc + Number
##
## Df Deviance AIC
## <none> 54.031 366.30
## - Number 1 59.487 369.75
## - Conc 1 119.757 430.02##
## Call: glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012342403,
## link = log)
##
## Coefficients:
## (Intercept) Conc Number
## 0.5933659 0.0008476 0.0123867
##
## Degrees of Freedom: 47 Total (i.e. Null); 45 Residual
## Null Deviance: 135.2
## Residual Deviance: 54.03 AIC: 368.3modelofinal2= glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012342403,
link = log)
modelofinal2##
## Call: glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012337563,
## link = log)
##
## Coefficients:
## (Intercept) Conc Number
## 0.5933400 0.0008476 0.0123869
##
## Degrees of Freedom: 47 Total (i.e. Null); 45 Residual
## Null Deviance: 135.2
## Residual Deviance: 54.03 AIC: 368.3summary(modelofinal2)##
## Call:
## glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012337563,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5356 -0.8736 -0.2645 0.3977 2.4899
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5933400 0.5972457 0.993 0.3205
## Conc 0.0008476 0.0001101 7.700 1.36e-14 ***
## Number 0.0123869 0.0051871 2.388 0.0169 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.0123) family taken to be 1)
##
## Null deviance: 135.249 on 47 degrees of freedom
## Residual deviance: 54.031 on 45 degrees of freedom
## AIC: 368.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.012
## Std. Err.: 0.480
##
## 2 x log-likelihood: -360.2981-pchisq(135.249- 54.031,2)## [1] 0#Rechace Ho. existe evidencia estadística de que el modelo binimial negativo se ajusta alos datos #Rechace Ho. exisite evidencia estadística de que la variable Conc(concentración de cianato de potasio), influye sobre el número de huevos de trucha muertos #Rechace Ho. exisite evidencia estadística de que la variable Number(número de huevos utilizados), influye sobre el número de huevos de trucha muertos
par(mfrow=c(1,2))
plot(abs(residuals(modelofinal2)))
abline(h=2,col="red")
plot(abs(residuals(modelofinal2,type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(modelofinal2)),abs(residuals(modelofinal2,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]#Los datos 38 y 39 se pueden considerar datos qtípicos
library(car)
influence.measures(modelofinal2)## Influence measures of
## glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012337563, link = log) :
##
## dfb.1_ dfb.Conc dfb.Nmbr dffit cov.r cook.d hat inf
## 1 0.002152 -0.00520 -0.000128 0.00780 1.109 2.54e-05 0.0353
## 2 0.076216 -0.06010 -0.051472 0.11422 1.095 6.93e-03 0.0426
## 3 0.026082 -0.04496 -0.008376 0.06970 1.099 2.38e-03 0.0355
## 4 -0.000357 -0.00274 0.001387 0.00402 1.116 6.69e-06 0.0411
## 5 -0.071023 0.01499 0.062971 -0.07431 1.188 2.02e-03 0.1024
## 6 -0.021911 0.03780 0.006207 -0.06174 1.098 1.31e-03 0.0330
## 7 0.024539 -0.04233 -0.006951 0.06915 1.096 2.35e-03 0.0330
## 8 -0.003355 0.05619 -0.019136 -0.08639 1.093 2.39e-03 0.0351
## 9 -0.048069 0.02683 0.034057 -0.06871 1.101 1.60e-03 0.0366
## 10 -0.040292 0.06120 0.010988 -0.11715 1.067 3.86e-03 0.0286
## 11 0.023041 0.03871 -0.039720 -0.07208 1.108 1.77e-03 0.0423
## 12 0.092953 -0.07283 -0.056359 0.16157 1.047 1.61e-02 0.0318
## 13 -0.334594 0.01281 0.296240 -0.36488 1.001 2.45e-02 0.0632
## 14 -0.242625 0.01544 0.210220 -0.27284 1.034 1.60e-02 0.0534
## 15 0.000687 -0.00634 0.004103 0.02080 1.095 1.89e-04 0.0239
## 16 -0.035089 0.00774 0.026288 -0.05004 1.098 8.83e-04 0.0308
## 17 -0.120220 0.02856 0.137705 0.19019 1.088 2.12e-02 0.0563
## 18 0.026191 0.05749 -0.012550 0.14550 1.042 1.30e-02 0.0264
## 19 0.125030 0.09158 -0.109999 0.22292 1.009 3.42e-02 0.0343
## 20 -0.012514 -0.01519 0.009360 -0.03701 1.097 5.00e-04 0.0279
## 21 0.107827 -0.15244 -0.096544 -0.21575 1.206 1.43e-02 0.1312 *
## 22 0.028425 -0.07948 -0.021953 -0.09794 1.196 3.41e-03 0.1102
## 23 -0.003280 -0.02159 0.005271 -0.02446 1.200 2.37e-04 0.1092 *
## 24 -0.006504 -0.04281 0.010453 -0.04851 1.199 8.99e-04 0.1092
## 25 0.026977 0.05711 -0.047902 -0.08824 1.118 2.61e-03 0.0523
## 26 -0.234750 -0.12117 0.274525 0.29824 1.425 4.50e-02 0.2620 *
## 27 0.029528 0.09641 -0.065269 -0.14437 1.089 5.96e-03 0.0460
## 28 0.073775 0.17087 -0.136564 -0.26177 1.033 1.48e-02 0.0506
## 29 0.013694 -0.04882 0.006162 0.07606 1.094 2.90e-03 0.0333
## 30 0.258476 0.20770 -0.333676 -0.41444 1.073 3.56e-02 0.1006
## 31 0.010683 0.11110 -0.054643 -0.17132 1.056 7.40e-03 0.0374
## 32 -0.224772 0.26763 0.111546 -0.46537 0.750 1.96e-02 0.0343 *
## 33 -0.087191 0.06346 0.055067 -0.14458 1.060 5.52e-03 0.0326
## 34 -0.084883 -0.06334 0.110367 0.13941 1.147 9.94e-03 0.0810
## 35 -0.520389 0.05108 0.477306 -0.53156 1.079 5.33e-02 0.1280
## 36 -0.008051 0.11951 -0.046473 -0.21307 1.001 9.33e-03 0.0303
## 37 0.182151 -0.03666 -0.139105 0.25176 0.972 4.70e-02 0.0320
## 38 0.158065 -0.09078 -0.076686 0.37643 0.753 1.40e-01 0.0236 *
## 39 0.268235 -0.05399 -0.204846 0.37075 0.843 1.24e-01 0.0320
## 40 -0.847424 -0.01304 0.784249 -0.87843 0.778 6.84e-02 0.1029 *
## 41 -0.060836 -0.03887 0.055079 -0.09619 1.095 2.93e-03 0.0385
## 42 0.070452 -0.05882 -0.092225 -0.21256 1.011 9.62e-03 0.0327
## 43 0.136488 -0.06431 -0.165073 -0.28783 0.975 1.51e-02 0.0397
## 44 -0.070997 -0.12853 0.041500 -0.32125 0.859 1.35e-02 0.0267
## 45 0.029076 -0.03779 -0.026328 -0.05498 1.237 1.16e-03 0.1362 *
## 46 0.017748 0.05313 -0.022853 0.06141 1.214 1.67e-03 0.1205 *
## 47 0.039155 -0.05536 -0.035058 -0.07834 1.228 2.27e-03 0.1312 *
## 48 0.009185 0.02580 -0.011676 0.02995 1.218 3.82e-04 0.1224 *influencePlot(modelofinal2)
A la luz de todos los criterio spropuestos en la literatura no existen datos que se puedan considerar influyentes en este conjunto de datos.
coef(modelofinal2)## (Intercept) Conc Number
## 0.5933399751 0.0008475838 0.0123868569round(exp(coef(modelofinal2)),2)## (Intercept) Conc Number
## 1.81 1.00 1.01round(ifelse(coef(modelofinal2)>0,exp(coef(modelofinal2)),1/exp(coef(modelofinal2))),2)## (Intercept) Conc Number
## 1.81 1.00 1.01#Si la concentración de de cianato de potasio se incrementa el número de huevos muertos aumenta a razón de 1 #Si se aument ala cantidad de huevos utilizados el número de huevos muertos aumenta a razón de 1.01
#Multicolinealidad
cor(data.frame(trout[1:5]))## Conc Number Dead ahora despues
## Conc 1.0000000 0.1584168 0.878001228 0.000000000 0.000000000
## Number 0.1584168 1.0000000 0.308213822 0.181479694 -0.181479694
## Dead 0.8780012 0.3082138 1.000000000 0.002010753 -0.002010753
## ahora 0.0000000 0.1814797 0.002010753 1.000000000 -1.000000000
## despues 0.0000000 -0.1814797 -0.002010753 -1.000000000 1.000000000#No se onservan correlaciones significativas en el modelo de el número de huevos de trucha muertos
#Predicción #Realice una estimaci´on con Con = 90, When = later y Number = 234.
huevos=data.frame(Conc=90,
despues=1,
Number=234)
predict(modelofinal2,huevos,type="response")## 1
## 35.45084#El número de huevos muertos por la toxina en concentración de 90 mg/litro, con aplicación después y utilizando 234 huevos es en promedio de 35.4.
#PUNTO3 #Un grupo de 24 aficionados al ciclismo de nivel principiante se dieron cita en un parque para una carrera de ciclismo que consiste en andar en bicicleta durante 30 minutos a lo largo de un circuito de 1.3 millas. #El ganador es el que pedalea la distancia m´as larga. La informaci´on registrada para cada participante es el g´enero (M/F), si tuvo una experiencia previa en carreras as´ı (s´ı/no), autoevaluaci´on de habilidades para terminar la carrera y hacerlo bien (en una escala de 10 puntos, siendo 10 la confianza m´as alta) y la distancia recorrida en bicicleta (en millas). [?] Los datos se presentan a continuaci´on:
library(readxl)## Warning: package 'readxl' was built under R version 4.2.3CICLISTAS <- read_excel("C:/Users/Maira/Desktop/CICLISTAS.xlsx")
View(CICLISTAS)CICLISTAS$fem=ifelse(CICLISTAS$gender=="F",1,0)
CICLISTASCICLISTAS=CICLISTAS[,-(1)]
CICLISTASCICLISTAS$exp=ifelse(CICLISTAS$priorexpr=="yes",1,0)
CICLISTASCICLISTAS=CICLISTAS[,-(1)]
CICLISTASmod3=lm(distance~.,data = CICLISTAS)
mod3##
## Call:
## lm(formula = distance ~ ., data = CICLISTAS)
##
## Coefficients:
## (Intercept) selfeval fem exp
## 1.5844 0.5740 -0.7378 0.4217summary(mod3)##
## Call:
## lm(formula = distance ~ ., data = CICLISTAS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0607 -0.6091 -0.1115 0.5757 3.4535
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.5844 1.1515 1.376 0.184054
## selfeval 0.5740 0.1299 4.418 0.000265 ***
## fem -0.7378 0.7280 -1.013 0.322987
## exp 0.4217 0.6353 0.664 0.514406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.453 on 20 degrees of freedom
## Multiple R-squared: 0.6187, Adjusted R-squared: 0.5615
## F-statistic: 10.82 on 3 and 20 DF, p-value: 0.0001941#coeficiente de correlación lineal simple
sqrt(0.6187)## [1] 0.78657491/0.7378## [1] 1.355381#Exisite una buena relación entre la distancia recorrida por los ciclistas y las demás variables #el modelo exoplica el 61.876 de las distancias recorridas por los ciclistas
#Prueba anova
#rechace Ho. exisite evidencia estadística de que el modelo de regresión lineal simple se ajusta a los datos
#Prueba de coeficientes estimados
#Rechace Ho. existe e videncia estadística de que la variable selfeval(autoevaluación) influye sobre la distancia recorrida # No rechace Ho. existe e videncia estadística de que la variable fem(género femenino) no influye sobre la distancia recorrida #No Rechace Ho. existe evidencia estadística de que la variable exp(experiencias previas)no influye sobre la distancia recorrida
#Por cada unidad de autoevaluación el recorrido aumenta a razón de 0.57 millas #Si el ciclista es de género femenino el recorrido disminuye a razón de 1.35millas #si el ciclista tiene experiencias previas en estas carreras el recorrido aumenta a razón de 0.42 millas
#Predicción
carrera=data.frame(selfeval=6,
fem=0,
exp=0)
predict(mod3,carrera,type="response")## 1
## 5.028669#La distancia recorrida por un hombre sin experiencia, con una auto evaluación de 6(1-10) durante 30 munutos es de 5 millas en promedio
par(mfrow=c(1,2))
plot(abs(residuals(mod3)))
abline(h=2,col="red")
plot(abs(residuals(mod3,type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(mod3)),abs(residuals(mod3,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]Los datos 3,10,22 se pueden considerar como datos atípicos para este conjunto de datos
library(car)
influence.measures(mod3)## Influence measures of
## lm(formula = distance ~ ., data = CICLISTAS) :
##
## dfb.1_ dfb.slfv dfb.fem dfb.exp dffit cov.r cook.d hat inf
## 1 -0.020557 0.01948 0.00208 0.0112 -0.02976 1.454 2.33e-04 0.1562
## 2 0.022845 -0.02165 -0.00231 -0.0124 0.03308 1.453 2.88e-04 0.1562
## 3 0.578707 -0.46650 -0.57386 -0.4763 -0.82808 0.897 1.58e-01 0.2006
## 4 0.012949 0.04704 -0.07747 -0.1243 -0.17174 1.391 7.69e-03 0.1480
## 5 -0.101796 0.04236 0.09914 0.0736 -0.11985 1.657 3.77e-03 0.2653 *
## 6 -0.550571 0.52173 0.45835 0.2991 0.68754 1.479 1.18e-01 0.3228
## 7 0.081900 -0.06863 0.00266 -0.0622 0.13925 1.350 5.06e-03 0.1179
## 8 -0.072693 0.04431 -0.02262 0.0881 -0.17164 1.281 7.64e-03 0.0957
## 9 0.039430 0.00630 0.04926 -0.1077 0.19647 1.242 9.95e-03 0.0894
## 10 0.631003 -0.36109 -0.79164 0.3540 1.18698 0.358 2.62e-01 0.1509 *
## 11 -0.003355 0.03022 0.03332 -0.0516 0.09608 1.339 2.42e-03 0.0991
## 12 -0.201725 0.04279 0.33476 -0.1647 -0.55253 0.970 7.30e-02 0.1377
## 13 0.029297 -0.06733 -0.05123 0.0623 -0.12624 1.370 4.17e-03 0.1248
## 14 0.003738 -0.01735 0.01241 0.0253 0.03723 1.486 3.65e-04 0.1749
## 15 0.089256 -0.27593 0.14085 -0.1094 -0.50778 1.224 6.42e-02 0.1943
## 16 0.008340 0.03030 -0.04990 -0.0801 -0.11061 1.420 3.21e-03 0.1480
## 17 -0.006678 0.06015 0.06632 -0.1028 0.19127 1.271 9.46e-03 0.0991
## 18 -0.117159 -0.37632 0.36212 0.4859 -0.78073 1.323 1.49e-01 0.3024
## 19 -0.005841 0.05261 0.05801 -0.0899 0.16729 1.292 7.27e-03 0.0991
## 20 -0.050560 0.04658 0.04359 -0.0159 -0.06449 1.716 1.09e-03 0.2862 *
## 21 0.029297 -0.06733 -0.05123 0.0623 -0.12624 1.370 4.17e-03 0.1248
## 22 -0.276366 0.85436 -0.43611 0.3386 1.57224 0.271 4.23e-01 0.1943 *
## 23 0.000325 -0.00151 0.00108 0.0022 0.00324 1.488 2.76e-06 0.1749
## 24 -0.189654 0.04023 0.31473 -0.1548 -0.51947 1.013 6.52e-02 0.1377influencePlot(mod3)
#A la luz de todos los criterios para es te conjunto de datos no se observan datos influyentes.
#PUNTO 4
#Los datos siguientes son de un estudio de 60 sujetos: 30 diab´eticos y 30 no diab´eticos que ten´ıa como objetivo determinar el efecto de la obesidad, los antecedentes familiares y practicar alg´un tipo de actividad f´ısica en el desarrollo de la diabetes. IMC hace referencia al ´ındice de masa corporal, ATF hace referencia a realizar (1=si) alg´un tipo de actividad f´ısica y HF hace referencia a los antecedentes familiares de la enfermedad.
remove(list=ls(all=TRUE))
library(readxl)
library(GLMsData)
library(MASS)
library(AER)
library(DescTools)
library(readxl)library(readxl)
Estudio_diabéticos <- read_excel("C:/Users/Maira/Desktop/bosque/Estudio diabéticos.xlsx")## New names:
## • `` -> `...1`View(Estudio_diabéticos)Estudio_diabéticos=Estudio_diabéticos[-(1)]
Estudio_diabéticosmod4=glm(HF1~IMC1+ATF1,data = Estudio_diabéticos,family = binomial())
mod4##
## Call: glm(formula = HF1 ~ IMC1 + ATF1, family = binomial(), data = Estudio_diabéticos)
##
## Coefficients:
## (Intercept) IMC1 ATF1
## -1.94987 0.07033 0.05450
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 41.05
## Residual Deviance: 40.01 AIC: 46.01step(mod4)## Start: AIC=46.01
## HF1 ~ IMC1 + ATF1
##
## Df Deviance AIC
## - ATF1 1 40.009 44.009
## - IMC1 1 40.783 44.783
## <none> 40.006 46.006
##
## Step: AIC=44.01
## HF1 ~ IMC1
##
## Df Deviance AIC
## - IMC1 1 41.054 43.054
## <none> 40.009 44.009
##
## Step: AIC=43.05
## HF1 ~ 1##
## Call: glm(formula = HF1 ~ 1, family = binomial(), data = Estudio_diabéticos)
##
## Coefficients:
## (Intercept)
## 0.2683
##
## Degrees of Freedom: 29 Total (i.e. Null); 29 Residual
## Null Deviance: 41.05
## Residual Deviance: 41.05 AIC: 43.05modelofinal4= glm(formula = HF1 ~ 1, family = binomial(), data = Estudio_diabéticos)
modelofinal4##
## Call: glm(formula = HF1 ~ 1, family = binomial(), data = Estudio_diabéticos)
##
## Coefficients:
## (Intercept)
## 0.2683
##
## Degrees of Freedom: 29 Total (i.e. Null); 29 Residual
## Null Deviance: 41.05
## Residual Deviance: 41.05 AIC: 43.05#predicción
diabetico=data.frame(HF1=1,
IMC1=26,
ATF1=1)
predict(modelofinal4,diabetico,type="response")## 1
## 0.5666667par(mfrow=c(1,2))
plot(abs(residuals(modelofinal4)))
abline(h=2,col="red")
plot(abs(residuals(modelofinal4,type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(modelofinal4)),abs(residuals(modelofinal4,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]library(car)
influence.measures(modelofinal4)## Influence measures of
## glm(formula = HF1 ~ 1, family = binomial(), data = Estudio_diabéticos) :
##
## dfb.1_ dffit cov.r cook.d hat inf
## 1 0.169 0.169 1.04 0.0273 0.0333
## 2 -0.206 -0.206 1.03 0.0466 0.0333
## 3 0.169 0.169 1.04 0.0273 0.0333
## 4 0.169 0.169 1.04 0.0273 0.0333
## 5 0.169 0.169 1.04 0.0273 0.0333
## 6 0.169 0.169 1.04 0.0273 0.0333
## 7 0.169 0.169 1.04 0.0273 0.0333
## 8 0.169 0.169 1.04 0.0273 0.0333
## 9 0.169 0.169 1.04 0.0273 0.0333
## 10 0.169 0.169 1.04 0.0273 0.0333
## 11 -0.206 -0.206 1.03 0.0466 0.0333
## 12 -0.206 -0.206 1.03 0.0466 0.0333
## 13 -0.206 -0.206 1.03 0.0466 0.0333
## 14 -0.206 -0.206 1.03 0.0466 0.0333
## 15 -0.206 -0.206 1.03 0.0466 0.0333
## 16 0.169 0.169 1.04 0.0273 0.0333
## 17 -0.206 -0.206 1.03 0.0466 0.0333
## 18 -0.206 -0.206 1.03 0.0466 0.0333
## 19 -0.206 -0.206 1.03 0.0466 0.0333
## 20 0.169 0.169 1.04 0.0273 0.0333
## 21 0.169 0.169 1.04 0.0273 0.0333
## 22 0.169 0.169 1.04 0.0273 0.0333
## 23 -0.206 -0.206 1.03 0.0466 0.0333
## 24 0.169 0.169 1.04 0.0273 0.0333
## 25 0.169 0.169 1.04 0.0273 0.0333
## 26 -0.206 -0.206 1.03 0.0466 0.0333
## 27 -0.206 -0.206 1.03 0.0466 0.0333
## 28 -0.206 -0.206 1.03 0.0466 0.0333
## 29 0.169 0.169 1.04 0.0273 0.0333
## 30 0.169 0.169 1.04 0.0273 0.0333influencePlot(modelofinal4)
#PUNTO 5 #El departamento de Administraci´on de la Salud est´a analizando la utilizaci´on de los recursos hospitalarios. Se recuperaron datos sobre 30 hospitales elegidos al azar en EEUU y se est´a analizando el porcentaje de pacientes de la sala de emergencias (ER) que fueron hospitalizados, la ubicaci´on del hospital (urbano o rural), tipo de hospital (p´ublico o privado) y el tama˜no del hospital por n´umero de camas. Los datos son:
remove(list=ls(all=TRUE))
library(car)
library(GLMsData)
library(carData)
library(betareg)
library(AER)library(readxl)
hospitales <- read_excel("C:/Users/Maira/Desktop/bosque/hospitales.xlsx")
View(hospitales)
head(hospitales)hospitales$rural=ifelse(hospitales$tip=="public",1,0)
hospitales$public=ifelse(hospitales$loc=="rural",1,0)
hospitaleshospitales=hospitales[,-(2:3)]
hospitaleshospitales$propor=hospitales$percent/hospitales$beds
hospitales$propor## [1] 0.30357143 0.27083333 0.62295082 0.25806452 0.22727273 0.04244482
## [7] 0.09433962 0.05479452 0.31168831 0.10526316 0.08176101 0.42857143
## [13] 0.15217391 0.19653179 0.49206349 0.04395604 0.07792208 0.16455696
## [19] 0.73214286 1.04651163 0.20312500 0.21761658 0.07713499 0.05166667
## [25] 0.10256410 0.13183280 0.13846154 0.29545455 0.09185804 0.22222222hospitales[20,5]=0.99999999
hospitaleslibrary(betareg)
mod5=betareg(propor ~ percent + beds + public + rural, data = hospitales)
mod5.1=betareg(propor~1,data = hospitales)## Warning in betareg.fit(X, Y, Z, weights, offset, link, link.phi, type,
## control): no valid starting value for precision parameter found, using 1
## insteadmod5.2=betareg(propor ~ percent, data = hospitales)## Warning in betareg.fit(X, Y, Z, weights, offset, link, link.phi, type,
## control): no valid starting value for precision parameter found, using 1
## insteadmod5.3=betareg(propor ~ beds, data = hospitales)## Warning in betareg.fit(X, Y, Z, weights, offset, link, link.phi, type,
## control): no valid starting value for precision parameter found, using 1
## insteadmod5.4=betareg(propor ~ public, data = hospitales)## Warning in betareg.fit(X, Y, Z, weights, offset, link, link.phi, type,
## control): no valid starting value for precision parameter found, using 1
## insteadmod5.5=betareg(propor ~ rural, data = hospitales)## Warning in betareg.fit(X, Y, Z, weights, offset, link, link.phi, type,
## control): no valid starting value for precision parameter found, using 1
## insteadAIC(mod5)## [1] -32.14938AIC(mod5.1)## [1] -9.821081AIC(mod5.2)## [1] -13.36709AIC(mod5.3)## [1] -12.34189AIC(mod5.4)## [1] -10.73227AIC(mod5.5)## [1] -8.849339modelofinal5=betareg(propor ~ percent + beds + public + rural, data = hospitales)
modelofinal5##
## Call:
## betareg(formula = propor ~ percent + beds + public + rural, data = hospitales)
##
## Coefficients (mean model with logit link):
## (Intercept) percent beds public rural
## -1.843594 0.075885 -0.008097 0.094683 0.496392
##
## Phi coefficients (precision model with identity link):
## (phi)
## 4.217summary(modelofinal5)##
## Call:
## betareg(formula = propor ~ percent + beds + public + rural, data = hospitales)
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -1.8983 -0.5022 -0.1803 0.3164 8.5381
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.843594 0.447840 -4.117 3.84e-05 ***
## percent 0.075885 0.014598 5.198 2.01e-07 ***
## beds -0.008097 0.001304 -6.211 5.25e-10 ***
## public 0.094683 0.374890 0.253 0.801
## rural 0.496392 0.348283 1.425 0.154
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 4.217 1.087 3.879 0.000105 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 22.07 on 6 Df
## Pseudo R-squared: 0.3087
## Number of iterations: 23 (BFGS) + 1 (Fisher scoring)summary(mod5.1)##
## Call:
## betareg(formula = propor ~ 1, data = hospitales)
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -0.8024 -0.5024 -0.2258 -0.0033 6.9605
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4011 0.2303 -1.742 0.0815 .
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 1.2021 0.2523 4.764 1.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 6.911 on 2 Df
## Number of iterations: 9 (BFGS) + 1 (Fisher scoring)#Existe una aceptable explicación de la proporcion de los pacientes hospitalizados por el modelo #Rechace Ho. Exisite evidencia estadística de que el porcentaje de pacientes en sala influye sobre la proporción de pacientes hospitalizados #Rechace Ho. Existe evidencia estadística de que el número de camas no influye sobre la proporción de pacientes hospitalizados #No rechace Ho.Exisite evidencia estadística de que la ubicación del hospital (rural) influye sobre la proporción de pacientes hospitalizados #No rechacce Ho.Exisite evidencia estadística de que eltipo hospital (pblic) influye sobre la proporción de pacientes hospitalizados
par(mfrow=c(1,2))
plot(abs(residuals(modelofinal5)))
abline(h=2,col="red")
plot(abs(residuals(modelofinal5,type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(modelofinal5)),abs(residuals(modelofinal5,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]#no se observan datos atípicos para este conjunto de datos
1-pchisq(22.07-6.911,4)## [1] 0.004382551#Rechace Ho. Exisit evidencia estadística de que el modelo de regresión beta se ajusta a los datos
exp(coef(modelofinal5))## (Intercept) percent beds public rural (phi)
## 0.1582476 1.0788380 0.9919353 1.0993108 1.6427835 67.81442031/0.9919353## [1] 1.00813#La proporción de pacientes hospitalizados aumenta a razón de 1.07 cuando se incrementa el porcentaje de pacientes en sala de emergencias. #Si el el hospital es de menor tamaño medido por su número de camas la proporción de pacientes hospitalizados disminuye e razón de 1.00813 #Si el hospital es de caracter público la cantidad de pacientes hospitalizados aumenta a razón de 1.09 #si el hospital tiene ubicación rural la cantidad de pacientes hospitalizados aumenta a razón de 1.64
#Predicción
paciente=data.frame(beds=50,
public=1,
rural=1,
percent=0)
predict(modelofinal5,paciente,type="response")## 1
## 0.160113#La proporción de pacientes hospitalizados en un hospital público con 50 camas y ubicación rural es de 16.01%
#PUNTO 6
#Una persona est´a pensando en comprar un carro nuevo. Realiza una b´usqueda en l´ınea y recopila informaci´on sobre las marcas y modelos que le gustan. Conjetura que las siguientes caracter´ısticas del autom´ovil pueden influir potencialmente en su precio: estilo de carrocer´ıa (cup´e, hatchback o sed´an), pa´ıs de fabricaci´on (EE. UU., Alemania o Jap´on), kilometraje en carretera (en mpg), n´umero de puertas (2 o 4) y si el interior es de cuero o no. [?] Los datos de estas variables y el precio (en d´olares estadounidenses) se dan a continuaci´on para 27 autom´oviles.
remove(list=ls(all=TRUE))
library(GLMsData)
library(MASS)
library(AER)
library(ISLR)## Warning: package 'ISLR' was built under R version 4.2.3library(car)
library(carData)
library(readxl)library(readxl)
AUTOS <- read_excel("C:/Users/Maira/Desktop/bosque/AUTOS.xlsx")
head(AUTOS)AUTOS$sit=ifelse(AUTOS$CUERO=="YES",1,0)
AUTOSAUTOS=AUTOS[,-(5)]
AUTOSAUTOS$paijap=ifelse(AUTOS$PAIS=="JAPAN",1,0)
AUTOS$paigerm=ifelse(AUTOS$PAIS=="GERMANY",1,0)
AUTOSAUTOS=AUTOS[,-(2)]
AUTOSAUTOS$estcup=ifelse(AUTOS$ESTILO=="COUPE",1,0)
AUTOS$estsed=ifelse(AUTOS$ESTILO=="SEDAN",1,0)
AUTOSAUTOS=AUTOS[,-(1)]
AUTOSmod6=glm(PRECIO~.,data = AUTOS)
mod6##
## Call: glm(formula = PRECIO ~ ., data = AUTOS)
##
## Coefficients:
## (Intercept) HWY PUERTAS sit paijap paigerm
## 8652.7 152.9 1413.0 10350.9 -2252.7 13.3
## estcup estsed
## 1146.8 6582.3
##
## Degrees of Freedom: 26 Total (i.e. Null); 19 Residual
## Null Deviance: 1.357e+09
## Residual Deviance: 256800000 AIC: 528.5step(mod6)## Start: AIC=528.46
## PRECIO ~ HWY + PUERTAS + sit + paijap + paigerm + estcup + estsed
##
## Df Deviance AIC
## - paigerm 1 256841336 526.46
## - estcup 1 260377556 526.83
## - HWY 1 272416014 528.05
## - paijap 1 275180021 528.32
## <none> 256840651 528.46
## - PUERTAS 1 300230384 530.68
## - estsed 1 412212540 539.24
## - sit 1 594834534 549.14
##
## Step: AIC=526.46
## PRECIO ~ HWY + PUERTAS + sit + paijap + estcup + estsed
##
## Df Deviance AIC
## - estcup 1 260619846 524.86
## - HWY 1 273379429 526.15
## <none> 256841336 526.46
## - paijap 1 277765319 526.58
## - PUERTAS 1 302231147 528.86
## - estsed 1 414569884 537.39
## - sit 1 615336350 548.05
##
## Step: AIC=524.86
## PRECIO ~ HWY + PUERTAS + sit + paijap + estsed
##
## Df Deviance AIC
## - HWY 1 273493888 524.16
## - paijap 1 279137426 524.71
## <none> 260619846 524.86
## - PUERTAS 1 305164932 527.12
## - estsed 1 445218368 537.32
## - sit 1 635535536 546.92
##
## Step: AIC=524.16
## PRECIO ~ PUERTAS + sit + paijap + estsed
##
## Df Deviance AIC
## - paijap 1 288853797 523.63
## <none> 273493888 524.16
## - PUERTAS 1 314561004 525.94
## - estsed 1 456876606 536.01
## - sit 1 642825514 545.23
##
## Step: AIC=523.63
## PRECIO ~ PUERTAS + sit + estsed
##
## Df Deviance AIC
## <none> 288853797 523.63
## - PUERTAS 1 334047035 525.56
## - estsed 1 501564304 536.53
## - sit 1 699180648 545.50##
## Call: glm(formula = PRECIO ~ PUERTAS + sit + estsed, data = AUTOS)
##
## Coefficients:
## (Intercept) PUERTAS sit estsed
## 13885 1402 8835 6361
##
## Degrees of Freedom: 26 Total (i.e. Null); 23 Residual
## Null Deviance: 1.357e+09
## Residual Deviance: 288900000 AIC: 523.6modelofinal6=glm(formula = PRECIO ~ PUERTAS + sit + estsed, data = AUTOS)
modelofinal6##
## Call: glm(formula = PRECIO ~ PUERTAS + sit + estsed, data = AUTOS)
##
## Coefficients:
## (Intercept) PUERTAS sit estsed
## 13885 1402 8835 6361
##
## Degrees of Freedom: 26 Total (i.e. Null); 23 Residual
## Null Deviance: 1.357e+09
## Residual Deviance: 288900000 AIC: 523.6summary(modelofinal6)##
## Call:
## glm(formula = PRECIO ~ PUERTAS + sit + estsed, data = AUTOS)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5953.2 -2135.3 -199.7 1894.4 8950.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13885.2 2507.9 5.537 1.24e-05 ***
## PUERTAS 1401.8 738.9 1.897 0.070458 .
## sit 8834.7 1545.6 5.716 8.05e-06 ***
## estsed 6360.9 1545.6 4.115 0.000422 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 12558861)
##
## Null deviance: 1357075463 on 26 degrees of freedom
## Residual deviance: 288853797 on 23 degrees of freedom
## AIC: 523.63
##
## Number of Fisher Scoring iterations: 21-pchisq(1357075463 -288853797,3)## [1] 0#Rechace Ho exisite evidencia estadística de que el modelo de regresión limeal ajusta a los datos
library(DescTools)
PseudoR2(modelofinal6,"Nagelkerke")## Nagelkerke
## 0.7871498#el modelo de regresión explica de manera excelente el comportamiento de la variable de interés
#datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modelofinal6)))
abline(h=2,col="red")
plot(abs(residuals(modelofinal6,type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(modelofinal6)),abs(residuals(modelofinal6,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]library(car)
influence.measures(modelofinal6)## Influence measures of
## glm(formula = PRECIO ~ PUERTAS + sit + estsed, data = AUTOS) :
##
## dfb.1_ dfb.PUER dfb.sit dfb.ests dffit cov.r cook.d hat inf
## 1 0.0444 -0.1094 0.0732 0.0732 -0.1940 1.237 0.009679 0.095
## 2 -0.0890 0.2192 -0.1467 -0.1467 0.3888 1.024 0.037090 0.095
## 3 0.2972 -0.2269 -0.0620 -0.0620 0.3285 1.190 0.027260 0.124
## 4 -0.0866 0.2134 -0.1428 -0.1428 0.3785 1.038 0.035262 0.095
## 5 -0.1034 0.2547 -0.1705 -0.1705 0.4518 0.941 0.049027 0.095
## 6 -0.4470 0.3413 0.0932 0.0932 -0.4941 1.008 0.059161 0.124
## 7 -0.3481 0.2658 0.0726 0.0726 -0.3848 1.132 0.036937 0.124
## 8 0.0874 -0.2154 0.1442 0.1442 -0.3821 1.033 0.035888 0.095
## 9 -0.0279 0.0420 0.1103 -0.0874 0.1479 1.438 0.005691 0.184
## 10 -0.4313 0.3294 0.0900 0.0900 -0.4767 1.029 0.055362 0.124
## 11 0.0411 -0.1012 0.0678 0.0678 -0.1796 1.249 0.008312 0.095
## 12 0.0614 -0.0469 -0.0128 -0.0128 0.0679 1.356 0.001202 0.124
## 13 0.0366 -0.0902 0.0604 0.0604 -0.1599 1.263 0.006613 0.095
## 14 -0.0658 0.1622 -0.1086 -0.1086 0.2877 1.147 0.020890 0.095
## 15 -0.0174 0.0429 -0.0287 -0.0287 0.0760 1.307 0.001507 0.095
## 16 0.2048 -0.1962 0.2455 -0.1352 0.3362 1.514 0.029103 0.256
## 17 -0.2574 0.2466 -0.3086 0.1699 -0.4225 1.463 0.045584 0.256
## 18 -0.0266 0.0401 0.1054 -0.0835 0.1413 1.440 0.005199 0.184
## 19 1.2931 -1.2388 -0.8536 1.5503 2.1229 0.248 0.738667 0.256 *
## 20 0.0145 -0.0111 -0.0291 -0.0291 -0.0622 1.405 0.001010 0.153
## 21 0.0385 -0.0294 -0.0773 -0.0773 -0.1653 1.372 0.007095 0.153
## 22 -0.0228 0.0219 0.0151 -0.0274 -0.0375 1.605 0.000368 0.256 *
## 23 -0.0391 0.0298 0.0784 0.0784 0.1678 1.371 0.007308 0.153
## 24 0.1764 -0.2656 0.5533 -0.6984 -0.9361 0.763 0.194608 0.184
## 25 -0.0602 0.0459 0.1208 0.1208 0.2586 1.319 0.017185 0.153
## 26 0.0934 -0.0713 -0.1876 -0.1876 -0.4014 1.203 0.040471 0.153
## 27 0.0343 -0.0517 0.1077 -0.1360 -0.1823 1.425 0.008624 0.184influencePlot(modelofinal6)
#A la luz de todos los criterios propuestos en la literatura no se observan datos que puedan considerarse influyentes
AUTOS#Predicicón
pric=data.frame( HWY = 0,
PUERTAS=4,
sit=1,
paijap=1,
estsed=1)
predict(modelofinal6,pric,type="response")## 1
## 34687.85#el precio del vehículo cero kilómetros, ensamblado en Japon,cojinería en cuero, estilo sedan, de 4 puertas cuesta en promedio 34687.85 dólares.
#PUNTO 7
#En una cl´ınica infantil se llev´o a cabo un estudio de intervenci´on y control sobre la obesidad infantil. Una cohorte de 36 ni˜nos obesos, de 6 a 16 a˜nos, fueron seguidos durante 9 meses. La intervenci´on consisti´o en sesiones educativas. para padres y actividades de ejercicio vigoroso para ni˜nos. Los participantes del grupo de control recibieron cursos sobre otros estilos de vida activos y saludables. Se tiene las variables g´enero (M/F), edad (en a˜nos), grupo (Tx de intervenci´on o Cx de control) y los IMC pre y post se muestran en la tabla a continuaci´on:
remove(list=ls(all=TRUE))
library(GLMsData)
library(MASS)
library(AER)
library(ISLR)
library(car)
library(carData)
library(readxl)library(readxl)
OBESIDAD_INF <- read_excel("C:/Users/Maira/Desktop/OBESIDAD INF.xlsx")
View(OBESIDAD_INF)
OBESIDAD_INFOBESIDAD_INF$GENF=ifelse(OBESIDAD_INF$Gender=="F",1,0)
OBESIDAD_INFOBESIDAD_INF=OBESIDAD_INF[,-(2)]
OBESIDAD_INFOBESIDAD_INF$control=ifelse(OBESIDAD_INF$Group=="Cx",1,0)
OBESIDAD_INFOBESIDAD_INF=OBESIDAD_INF[,-(3)]
OBESIDAD_INFOBESIDAD_INF=OBESIDAD_INF[,-(1)]
OBESIDAD_INFmod7=glm(DIF~ Age,data = OBESIDAD_INF,family = inverse.gaussian(link = "log"))
mod7##
## Call: glm(formula = DIF ~ Age, family = inverse.gaussian(link = "log"),
## data = OBESIDAD_INF)
##
## Coefficients:
## (Intercept) Age
## 0.28885 0.09073
##
## Degrees of Freedom: 35 Total (i.e. Null); 34 Residual
## Null Deviance: 9.857
## Residual Deviance: 9.216 AIC: 164.8mod7.1=glm(DIF~ Age,data = OBESIDAD_INF,family = Gamma(link = "log"))sd(OBESIDAD_INF$DIF)/abs(mean(OBESIDAD_INF$DIF))## [1] 0.9029153según metodo descriptivo se debe utilizar modelo inverse gaussian
log(logLik(mod7)/ logLik(mod7.1))## 'log Lik.' -0.001390006 (df=3)#al se r inferior a 0 se recomienda usar un modelo de regresión Gamma
step(mod7.1)## Start: AIC=165.06
## DIF ~ Age
##
## Df Deviance AIC
## <none> 22.288 165.06
## - Age 1 24.361 165.95##
## Call: glm(formula = DIF ~ Age, family = Gamma(link = "log"), data = OBESIDAD_INF)
##
## Coefficients:
## (Intercept) Age
## 0.37590 0.08255
##
## Degrees of Freedom: 35 Total (i.e. Null); 34 Residual
## Null Deviance: 24.36
## Residual Deviance: 22.29 AIC: 165.1modelofinal7=glm(formula = DIF ~ Age, family = Gamma(link = "log"), data = OBESIDAD_INF)
modelofinal7##
## Call: glm(formula = DIF ~ Age, family = Gamma(link = "log"), data = OBESIDAD_INF)
##
## Coefficients:
## (Intercept) Age
## 0.37590 0.08255
##
## Degrees of Freedom: 35 Total (i.e. Null); 34 Residual
## Null Deviance: 24.36
## Residual Deviance: 22.29 AIC: 165.1round(exp(modelofinal7$coefficients),2)## (Intercept) Age
## 1.46 1.09library(DescTools)
PseudoR2(modelofinal7,"Nagelkerke")## Nagelkerke
## 0.09455913summary(modelofinal7)##
## Call:
## glm(formula = DIF ~ Age, family = Gamma(link = "log"), data = OBESIDAD_INF)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6412 -0.7383 -0.3533 0.1797 1.6385
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.37590 0.52676 0.714 0.480
## Age 0.08255 0.04591 1.798 0.081 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.7162977)
##
## Null deviance: 24.361 on 35 degrees of freedom
## Residual deviance: 22.288 on 34 degrees of freedom
## AIC: 165.06
##
## Number of Fisher Scoring iterations: 51-pchisq(24.361 -22.288,1)## [1] 0.1499264#No rechace HO. exisite evidencia estadística de que el mmodelo de regresión Gamma no se ajusta a los datos
par(mfrow=c(1,2))
plot(abs(residuals(modelofinal7)))
abline(h=2,col="red")
plot(abs(residuals(modelofinal7,type="pearson")))
abline(h=2,col="red")
residuos=data.frame(abs(residuals(modelofinal7)),abs(residuals(modelofinal7,type="pearson")))
residuos[residuos[,1]>2&residuos[,2]>2,]OBESIDAD_INFestimacion=data.frame(Age=10,
control=1 )
predict(modelofinal7,estimacion,type="response")## 1
## 3.324685#la diferencia esperada para una niña de 10 años es de 3.32 en promedio