predicción de un modelo simple en glm
#En este experimento se está probando la eficacia de dosis de un componente para generar mortalidad en ratas.
dose<-c(0:10)#suponga que las unidades son mg
status<-c("lived","lived","lived","died","died", "lived","died","died", "died","died","died")
data.frame(dose,status)->rats
attach(rats)
## The following objects are masked _by_ .GlobalEnv:
##
## dose, status
boxplot(dose~status)
#Modelo
rats.lr<-glm(status~dose,data=rats,family="binomial")
summary(rats.lr)
##
## Call:
## glm(formula = status ~ dose, family = "binomial", data = rats)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3553 -0.3775 -0.1260 0.3877 1.7726
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0295 2.0812 1.456 0.1455
## dose -0.8735 0.5132 -1.702 0.0887 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 14.4206 on 10 degrees of freedom
## Residual deviance: 7.1148 on 9 degrees of freedom
## AIC: 11.115
##
## Number of Fisher Scoring iterations: 6
anova(rats.lr, test="Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: status
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 10 14.4206
## dose 1 7.3058 9 7.1148 0.006873 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Plot
library(effects)
plot(allEffects(rats.lr))
#predecir
dose.new<-data.frame(dose=seq(0.5,10,1))
dose.new
## dose
## 1 0.5
## 2 1.5
## 3 2.5
## 4 3.5
## 5 4.5
## 6 5.5
## 7 6.5
## 8 7.5
## 9 8.5
## 10 9.5
p<-predict(rats.lr,dose.new,type="response")
cbind(dose.new,p)#se denota entonces que a dosis menores la probabilidad de sobrevivencia es mayor , que a dosis altas
## dose p
## 1 0.5 0.930393652
## 2 1.5 0.848032032
## 3 2.5 0.699673858
## 4 3.5 0.493060926
## 5 4.5 0.288791340
## 6 5.5 0.144950779
## 7 6.5 0.066095830
## 8 7.5 0.028699072
## 9 8.5 0.012185179
## 10 9.5 0.005123513
#Forma gráfica de la predicción.
fit.pre <- predict(rats.lr, type='response', se.fit=T)
plot(rats$dose, fit.pre$fit, xlab='Dosis', ylab='Probabilidad de sobrevivencia', ylim=c(0,1), type="l")
#Si mi interés es buscar ciertas valores específicos
dose.new=data.frame(dose=c(5,10,15))
dose.new
## dose
## 1 5
## 2 10
## 3 15
p=predict(rats.lr,dose.new,type="response")
cbind(dose.new,p)
## dose p
## 1 5 2.078370e-01
## 2 10 3.316478e-03
## 3 15 4.220006e-05
Predicción de un modelo con dos factores
library(MASS)
head(birthwt)
## low age lwt race smoke ptl ht ui ftv bwt
## 85 0 19 182 2 0 0 0 1 0 2523
## 86 0 33 155 3 0 0 0 0 3 2551
## 87 0 20 105 1 1 0 0 0 1 2557
## 88 0 21 108 1 1 0 0 1 2 2594
## 89 0 18 107 1 1 0 0 1 0 2600
## 91 0 21 124 3 0 0 0 0 0 2622
nacimiento <- glm(low ~ age + race, family = binomial, data = birthwt)
summary(nacimiento)
##
## Call:
## glm(formula = low ~ age + race, family = binomial, data = birthwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1054 -0.8820 -0.7607 1.3466 1.8182
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.33286 0.86012 -0.387 0.699
## age -0.04355 0.03219 -1.353 0.176
## race 0.28600 0.17413 1.642 0.101
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.21 on 186 degrees of freedom
## AIC: 235.21
##
## Number of Fisher Scoring iterations: 4
pi.hat<- predict.glm(nacimiento, data.frame(age=25, race=3), type="response", se.fit=TRUE)
pi.hat$fit
## 1
## 0.3627181
library(effects)
plot(allEffects(nacimiento))
