POLECENIE: oszacuj uporządkowany model logitowy dla zmiennej zależnej opisującej stan zdrowia. Zinterpretuj model z wykorzystaniem ilorazu szans i wykresów pokazujących oszacowane prawdopodobieństwa dla różnych wartości zmiennych objaśniających. Oblicz prognozę dla przykładowych wartości zm. objaśniającej/objaśniających. Wykorzystaj poznane miary dopasowania (AIC, pseudoR2) przy porównywaniu modeli z różnym zestawem zmiennych objaśniających.

library(AER)
## Warning: pakiet 'AER' został zbudowany w wersji R 4.1.3
## Ładowanie wymaganego pakietu: car
## Warning: pakiet 'car' został zbudowany w wersji R 4.1.3
## Ładowanie wymaganego pakietu: carData
## Warning: pakiet 'carData' został zbudowany w wersji R 4.1.3
## Ładowanie wymaganego pakietu: lmtest
## Warning: pakiet 'lmtest' został zbudowany w wersji R 4.1.3
## Ładowanie wymaganego pakietu: zoo
## Warning: pakiet 'zoo' został zbudowany w wersji R 4.1.3
## 
## Dołączanie pakietu: 'zoo'
## Następujące obiekty zostały zakryte z 'package:base':
## 
##     as.Date, as.Date.numeric
## Ładowanie wymaganego pakietu: sandwich
## Ładowanie wymaganego pakietu: survival
library(MASS)
## Warning: pakiet 'MASS' został zbudowany w wersji R 4.1.3
library(effects)
## Warning: pakiet 'effects' został zbudowany w wersji R 4.1.3
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(reshape2) 
## Warning: pakiet 'reshape2' został zbudowany w wersji R 4.1.3
library(ggplot2)
library(pscl)
## Warning: pakiet 'pscl' został zbudowany w wersji R 4.1.3
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
data("NMES1988")
is.ordered(NMES1988$health)
## [1] FALSE
mydata=data.frame(NMES1988)
mydata$health_ord = factor(mydata$health, ordered = TRUE, levels = c("poor", "average", "excellent"))
is.ordered((mydata$health_ord))
## [1] TRUE
table(mydata$health_ord)
## 
##      poor   average excellent 
##       554      3509       343
model1=polr(health_ord ~ chronic +  insurance, data = mydata, Hess = TRUE)
summary(model1)
## Call:
## polr(formula = health_ord ~ chronic + insurance, data = mydata, 
##     Hess = TRUE)
## 
## Coefficients:
##                Value Std. Error t value
## chronic      -0.6584    0.02986 -22.049
## insuranceyes  0.6582    0.09225   7.136
## 
## Intercepts:
##                   Value    Std. Error t value 
## poor|average       -2.7628   0.1060   -26.0762
## average|excellent   2.2708   0.1008    22.5282
## 
## Residual Deviance: 5041.257 
## AIC: 5049.257
results = data.frame(summary(model1)$coefficients)
results$pval = (pnorm(abs(results$t.value), lower.tail = FALSE) * 2)
results$OR = exp(results$Value)
results
##                        Value Std..Error    t.value          pval         OR
## chronic           -0.6583636 0.02985893 -22.049136 9.736553e-108 0.51769779
## insuranceyes       0.6582204 0.09224582   7.135505  9.643280e-13 1.93135232
## poor|average      -2.7628023 0.10595118 -26.076181 6.792861e-150 0.06311466
## average|excellent  2.2707876 0.10079771  22.528166 2.198735e-112 9.68702710
1-0.51769779
## [1] 0.4823022

Obie zmienne objasniajace sa statyscznie istone na standardowym poziomie istotnosci. Ilorazy szans: dla chronic iloraz szans wynosi ok 0.52, czyli wraz ze wwzrostem 1. przewleklych chorob o 1 szanse na lepsze samopoczucie/zdorwie (na przejscie do wyzszej kategorii tego zdrowia) maleją o ok 48.23% ceteris paribus.

Dla insurance ilośc szans wynosi 1.93, czyli przy posiadaniu prywatnego uzbezpiecznia szanse na lepsze samopoczucie wzrastają o ok 93%, ceteris paribus.

PREDYJCA

example = data.frame(insurance = "yes", chronic = 3)
predict(model1, newdata = example, type = "probs")
##      poor   average excellent 
## 0.1906285 0.7824532 0.0269184

MIARY DOPASOWANIA: AIC” 5049.257

pR2(model1) #library(pscl)
## fitting null model for pseudo-r2
##           llh       llhNull            G2      McFadden          r2ML 
## -2520.6285211 -2823.2025838   605.1481254     0.1071741     0.1283317 
##          r2CU 
##     0.1776496

McFadden

Effect(focal.predictors = "insurance", mod = model1) # library(effects)
## 
## insurance effect (probability) for poor
## insurance
##         no        yes 
## 0.14834848 0.08272892 
## 
## insurance effect (probability) for average
## insurance
##        no       yes 
## 0.8155961 0.8498978 
## 
## insurance effect (probability) for excellent
## insurance
##         no        yes 
## 0.03605542 0.06737331
plot(Effect(focal.predictors = "chronic", mod = model1))

plot(Effect(focal.predictors = c("chronic", "insurance"),model1))

eff=Effect(focal.predictors = c("chronic", "insurance"),model1)
plotdata = data.frame(eff$x,eff$prob)
plotdata2 = melt(plotdata, id.vars = c("chronic", "insurance"),
                 variable.name = "health", value.name="Probability") #library(reshape2)
ggplot(data = plotdata2, aes(x = chronic, y = Probability, color = health)) +       
  geom_line(aes(group = health))  + facet_wrap(vars(insurance), labeller="label_both")

# dla dwoch zm. jakosciowych np. insurance i … skorzystaj z facet_grid zamiast facet_wrap

eff=Effect(focal.predictors = c(“chronic”, “insurance”, “gender”),model1) plotdata = data.frame(eff\(x,eff\)prob) plotdata2 = melt(plotdata, id.vars = names(plotdata)[1:3], variable.name = “health”, value.name=“Probability”) ggplot(data = plotdata2, aes(x = chronic, y = Probability, color = health)) +
geom_line(aes(group = health)) + facet_grid(insurance ~ gender, labeller=“label_both”)