Глава 7: Бинарни променливи и логистичка регресија

Задача 1

Ова се податоците од задачата во R:

leukemia_deaths = c(13, 5, 5, 3, 4, 18)
other_deaths = c(378, 200, 151, 47, 31, 33)
deaths = c(leukemia_deaths, other_deaths)
radiation_dose = c(0, 1, 10, 50, 100, 200)
radiation = rep(radiation_dose, 2)
leukemia = c(rep(1, 6), rep(0, 6))

leukemia_df = data.frame(deaths, radiation, leukemia)
leukemia_df

Во R треба да се подготват податоците за да се креира binomial модел.

proportion = cbind(leukemia_deaths, other_deaths)
leukemia_proportion_df = data.frame(proportion, radiation_dose)
leukemia_proportion_df

(a)

proportion треба да е веројатност, па затоа користиме logit и probit link функции.

summary(glm(
  proportion~radiation_dose, 
  family=binomial(link=logit), 
  data=leukemia_proportion_df
))

Call:
glm(formula = proportion ~ radiation_dose, family = binomial(link = logit), 
    data = leukemia_proportion_df)

Deviance Residuals: 
       1         2         3         4         5         6  
 0.41428  -0.48994  -0.13991   0.02835   0.00048   0.00269  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -3.488973   0.204062 -17.098  < 2e-16 ***
radiation_dose  0.014410   0.001817   7.932 2.15e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 54.35089  on 5  degrees of freedom
Residual deviance:  0.43206  on 4  degrees of freedom
AIC: 26.097

Number of Fisher Scoring iterations: 4

summary(glm(
  proportion~radiation_dose, 
  family=binomial(link=probit), 
  data=leukemia_proportion_df
))

Call:
glm(formula = proportion ~ radiation_dose, family = binomial(link = probit), 
    data = leukemia_proportion_df)

Deviance Residuals: 
       1         2         3         4         5         6  
 0.49118  -0.44422  -0.15683  -0.12140  -0.21161   0.09856  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -1.895922   0.088925 -21.321  < 2e-16 ***
radiation_dose  0.007504   0.001002   7.488 6.97e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 54.35089  on 5  degrees of freedom
Residual deviance:  0.53242  on 4  degrees of freedom
AIC: 26.197

Number of Fisher Scoring iterations: 3

(b)

Вредноста на AIC e пониска за logit линк функцијата. Доволно е ниска за да можеме да заклучиме дека моделот добро ги научил податоците.

Проблематична е поделбата на самите податоци во задачата. Пример, имаме податоци за смртноста од леукемија за дози на радијација од 1 до 9 (9 посебни дози) и имаме дози на радијација од 100 до 199 (100 посебни дози). Двете се претставени со една точка и тоа малку ја иместува правата на регресија, давајќи поголема значајност на точките што се подесно. Доколку податоците беа поделени на подеднаков број на дози на радијација (1-9, па 10-19, па 20-29, …, 190-199, 200+), регресијата ќе беше доста поточна.

(c)

Моделот ја крериа следнава функција за резултатите: \(logit(-3.48 + 0.014*radiation\_dose)\). Резултатите секогаш ќе бидат помеѓу 0 и 1, благодарение на логит функцијата.

Како што расте дозата на радијација, расте и пропорцијата на умрени од леукемија. За зголемување од една радијациска доза, порастот на смртност од леукемија е очекувано да биде отприлика \(logit(\frac{radiation\_dose + 1}{radiation\_dose})\).

Задача 3

Ова се податоците за мажите:

men_st = c(
  18, 22, 16, 30, 9, 14, 10, 16, 
  16, 23, 13, 22, 9, 12, 7, 11, 
  7, 17, 11, 25, 12, 19, 12, 15, 
  12, 25, 12, 14, 12, 15, 8, 9, 
  24, 50, 8, 12, 20, 28, 5, 7, 
  16, 21, 11, 20, 16, 21, 1, 2, 
  22, 32, 4, 10, 25, 31, 16, 22, 
  12, 14, 4, 12, 32, 38, 19, 25, 
  22, 34, 0, 0, 4, 5, 0, 0, 
  28, 37, 13, 23, 25, 31, 25, 35
)
men_survival = men_st[seq(1, length(men_st), 2)]
men_death = men_st[seq(2, length(men_st), 2)]
men_sd = cbind(men_survival, men_death)

men_year = c(
  as.vector(
    mapply(function (x) { 
      return(rep(x, 4)) 
    }, 
    seq(1938, 1947))
  )
)
men_faculty = factor(c(rep(c('M', 'A', 'S', 'E'), 10)))

men_df = data.frame(men_year, men_faculty, men_sd)
men_df

Ова се податоците за жените:

women_st = c(
  14, 19, 1, 1, 
  11, 16, 4, 4, 
  15, 18, 6, 7, 
  15, 21, 3, 3, 
  8, 9, 4, 4, 
  13, 13, 8, 9, 
  18, 22, 5, 5, 
  18, 22, 16, 17, 
  1, 1, 1, 1, 
  13, 16, 10, 10
)
women_survival = women_st[seq(1, length(women_st), 2)]
women_death = women_st[seq(2, length(women_st), 2)]
women_sd = cbind(women_survival, women_death)

women_year = as.vector(
  mapply(function (x) { 
    return(rep(x, 2)) 
  }, 
  seq(1938, 1947))
)
women_faculty = factor(rep(c('A', 'S'), 10))

women_df = data.frame(women_year, women_faculty, women_sd)
women_df

Вкупните податоци по година и факултет се:

total_survival = men_df$men_survival
total_death = men_df$men_death
total_sd = cbind(total_survival, total_death)
total_year = men_df$men_year
total_faculty = men_df$men_faculty
total_df = data.frame(total_year, total_faculty, total_sd)

for (faculty in c('A', 'S')) {
  total_df[total_df$total_faculty == faculty,]$total_survival = 
    total_df[total_df$total_faculty == faculty,]$total_survival + 
    women_df[women_df$women_faculty == faculty,]$women_survival
  
  total_df[total_df$total_faculty == faculty,]$total_death = 
    total_df[total_df$total_faculty == faculty,]$total_death + 
    women_df[women_df$women_faculty == faculty,]$women_death
}

total_df

Вкупните податоци по година се:

total_survival_year = rep(0, 10)
total_death_year = rep(0, 10)
for (faculty in c('S', 'M', 'A', 'E')) {
  total_survival_year = 
    total_survival_year + 
    total_df[total_df$total_faculty == faculty,]$total_survival
  
  total_death_year = 
    total_death_year + 
    total_df[total_df$total_faculty == faculty,]$total_death
}

total_year_year = as.vector(seq(1938, 1947))
total_year_sd = cbind(total_survival_year, total_death_year)

total_year_df = data.frame(total_year_year, total_year_sd)
total_year_df

Вкупните податоци по година за мажи се:

total_survival_men = rep(0, 10)
total_death_men = rep(0, 10)
for (faculty in c('S', 'M', 'A', 'E')) {
  total_survival_men = 
    total_survival_men + 
    men_df[men_df$men_faculty == faculty,]$men_survival
  
  total_death_men = 
    total_death_men + 
    men_df[men_df$men_faculty == faculty,]$men_death
}

total_men_year = as.vector(seq(1938, 1947))
total_men_sd = cbind(total_survival_men, total_death_men)

total_men_df = data.frame(total_men_year, total_men_sd)
total_men_df

Вкупните податоци по година за жени се:

total_survival_women = rep(0, 10)
total_death_women = rep(0, 10)
for (faculty in c('S', 'A')) {
  total_survival_women = 
    total_survival_women + 
    women_df[women_df$women_faculty == faculty,]$women_survival
  
  total_death_women = 
    total_death_women + 
    women_df[women_df$women_faculty == faculty,]$women_death
}

total_women_year = as.vector(seq(1938, 1947))
total_women_sd = cbind(total_survival_women, total_death_women)

total_women_df = data.frame(total_women_year, total_women_sd)
total_women_df

Конечно, податоците за стапката на смртност на мажи и жени за Arts и Sciences се:

as_year = c(
  men_df[men_df$men_faculty == 'A', ]$men_year,
  men_df[men_df$men_faculty == 'S', ]$men_year,
  women_df[women_df$women_faculty == 'A', ]$women_year,
  women_df[women_df$women_faculty == 'S', ]$women_year
)

as_survival = c(
  men_df[men_df$men_faculty == 'A', ]$men_survival,
  men_df[men_df$men_faculty == 'S', ]$men_survival,
  women_df[women_df$women_faculty == 'A', ]$women_survival,
  women_df[women_df$women_faculty == 'S', ]$women_survival
)

as_death = c(
  men_df[men_df$men_faculty == 'A', ]$men_death,
  men_df[men_df$men_faculty == 'S', ]$men_death,
  women_df[women_df$women_faculty == 'A', ]$women_death,
  women_df[women_df$women_faculty == 'S', ]$women_death
)

as_faculty = as.factor(c(
  rep('A', length(men_df[men_df$men_faculty == 'A', ]$men_faculty)),
  rep('S', length(men_df[men_df$men_faculty == 'S', ]$men_faculty)),
  rep('A', length(women_df[women_df$women_faculty == 'A', ]$women_faculty)),
  rep('S', length(women_df[women_df$women_faculty == 'S', ]$women_faculty))
))

as_proportion = cbind(as_survival, as_death)
as_df = data.frame(as_year, as_faculty, as_proportion)
as_df

(a)

Прво, ги плотираме пропроциите секоја година, за да видиме како се движи пропорцијата. Речиси е јасно дека нема никаква промена низ годините.

total_year_df["proportion"] = 
  total_year_df$total_survival_year / (
    total_year_df$total_survival_year + 
      total_year_df$total_death_year
    )

with(total_year_df, plot(total_year_year, proportion))

Моделот креиран за вкупните податоци по години добро ги предвидува податоците, Slope-от е доста низок, a и p-вредноста е многу голема, па не ја отфрламе хипотезата дека slope-от e нула т.е. дека низ годините е иста пропорцијата.

summary(glm(
  total_year_sd ~ total_year_year, 
  family=binomial(link="logit"), 
  data=total_year_df
))

Call:
glm(formula = total_year_sd ~ total_year_year, family = binomial(link = "logit"), 
    data = total_year_df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7031  -0.3309   0.1485   0.2171   0.5067  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)     -35.55903   32.50748  -1.094    0.274
total_year_year   0.01813    0.01673   1.083    0.279

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2.5983  on 9  degrees of freedom
Residual deviance: 1.4238  on 8  degrees of freedom
AIC: 60.695

Number of Fisher Scoring iterations: 3

До истиот заклучок доаѓаме и со прост модел кој ја предвидува пропорцијата. Повторно не можеме да ја отфрлиме хипотезата дека slope-от е нула, т.е., дека низ годините нема промена на пропорцијата. Дополнително, и R2 е доста ниско.

summary(lm(
  proportion ~ total_year_year, 
  data=total_year_df
))

Call:
lm(formula = proportion ~ total_year_year, data = total_year_df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.024817 -0.014108  0.005449  0.008852  0.019581 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)     -6.823890   3.541230  -1.927   0.0901 .
total_year_year  0.003725   0.001823   2.043   0.0753 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01656 on 8 degrees of freedom
Multiple R-squared:  0.3429,    Adjusted R-squared:  0.2608 
F-statistic: 4.175 on 1 and 8 DF,  p-value: 0.07528

(b)

Прво, ги плотираме пропроциите секоја година, за да видиме како се движи пропорцијата. Од самата слика се забележува дека и да има некоја ситна позитивна промена, таа е премногу мала.

total_men_df["proportion"] = 
  total_men_df$total_survival_men / (
    total_men_df$total_survival_men + 
      total_men_df$total_death_men
    )

with(total_men_df, plot(total_men_year, proportion))

Моделот креиран за вкупните податоци по години добро ги предвидува податоците, Slope-от е доста низок, a и p-вредноста е многу голема, па не ја отфрламе хипотезата дека slope-от e нула т.е. дека низ годините е иста пропорцијата.

summary(glm(
  total_men_sd ~ total_men_year, 
  family=binomial(link="logit"), 
  data=total_men_df
))

Call:
glm(formula = total_men_sd ~ total_men_year, family = binomial(link = "logit"), 
    data = total_men_df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7412  -0.2076   0.1373   0.2980   0.4545  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)
(Intercept)    -36.87246   36.84660  -1.001    0.317
total_men_year   0.01877    0.01897   0.990    0.322

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2.6682  on 9  degrees of freedom
Residual deviance: 1.6873  on 8  degrees of freedom
AIC: 58.303

Number of Fisher Scoring iterations: 3

До истиот заклучок доаѓаме и со прост модел кој ја предвидува пропорцијата. Повторно не можеме да ја отфрлиме хипотезата дека slope-от е нула, т.е., дека низ годините нема промена на пропорцијата. Дополнително, и R2 е доста ниско.

summary(lm(proportion ~ total_men_year, data=total_men_df))

Call:
lm(formula = proportion ~ total_men_year, data = total_men_df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.033779 -0.010077  0.006668  0.012150  0.019612 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)    -7.433755   4.249935  -1.749    0.118
total_men_year  0.004033   0.002188   1.843    0.103

Residual standard error: 0.01987 on 8 degrees of freedom
Multiple R-squared:  0.2981,    Adjusted R-squared:  0.2103 
F-statistic: 3.397 on 1 and 8 DF,  p-value: 0.1025

(c)

Прво, ги плотираме пропроциите секоја година, за да видиме како се движи пропорцијата. Од самата слика се забележува дека е можно да има некаква позитивна промена, но дека не е многу силна.

total_women_df["proportion"] = 
  total_women_df$total_survival_women / (
    total_women_df$total_survival_women + 
      total_women_df$total_death_women
    )

with(total_women_df, plot(total_women_year, proportion))

Моделот креиран за вкупните податоци по години добро ги предвидува податоците, Slope-от е доста низок, a и p-вредноста е многу голема, па не ја отфрламе хипотезата дека slope-от e нула т.е. дека низ годините е иста пропорцијата. Сепак, доказите против хипотезата се посилни овде.

summary(glm(
  total_women_sd ~ total_women_year, 
  family=binomial(link="logit"), 
  data=total_women_df
))

Call:
glm(formula = total_women_sd ~ total_women_year, family = binomial(link = "logit"), 
    data = total_women_df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.27766  -0.12300  -0.06279   0.14367   0.37813  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)
(Intercept)      -38.35152   69.62339  -0.551    0.582
total_women_year   0.01965    0.03584   0.548    0.583

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 0.67289  on 9  degrees of freedom
Residual deviance: 0.37190  on 8  degrees of freedom
AIC: 44.174

Number of Fisher Scoring iterations: 3

До истиот заклучок доаѓаме и со прост модел кој ја предвидува пропорцијата. Повторно не можеме да ја отфрлиме хипотезата дека slope-от е нула, т.е., дека низ годините нема промена на пропорцијата. Дополнително, и R2 е доста ниско. Сепак, доказите против хипотезата се посилни овде.

summary(lm(proportion ~ total_women_year, data=total_women_df))

Call:
lm(formula = proportion ~ total_women_year, data = total_women_df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.022766 -0.010391 -0.007051  0.016238  0.024718 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)  
(Intercept)      -11.502437   3.978943  -2.891   0.0202 *
total_women_year   0.006159   0.002048   3.007   0.0169 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01861 on 8 degrees of freedom
Multiple R-squared:  0.5305,    Adjusted R-squared:  0.4718 
F-statistic:  9.04 on 1 and 8 DF,  p-value: 0.0169

(d)

Од генералниот линеарен модел забележуваме дека коефициентот пред факторот факултет е повисок од 0. Заклучокот е дека најверојатно има разлика помеѓу факултетите, т.е., дека поголема е смртноста кај тие од Sciences.

summary(glm(
  as_proportion ~ as_year + as_faculty, 
  family=binomial(link="logit"), 
  data=as_df
))

Call:
glm(formula = as_proportion ~ as_year + as_faculty, family = binomial(link = "logit"), 
    data = as_df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3686  -0.2279   0.1518   0.3336   0.9718  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -31.74089   43.94010  -0.722    0.470
as_year       0.01614    0.02263   0.713    0.476
as_facultyS   0.16444    0.12905   1.274    0.203

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12.5073  on 38  degrees of freedom
Residual deviance:  9.8258  on 36  degrees of freedom
AIC: 151.16

Number of Fisher Scoring iterations: 3