Глава 6: Нормални линеарни модели

Задача 1

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

period	refined	not_refined
1936–39	32.0	16.3
1946–49	31.2	23.1
1956–59	27.0	23.6
1966–69	21.0	27.7
1976–79	14.9	34.6
1986–89	8.8	33.9

Со следниов код ги подготвуваме податоците во R.

period = c(1939, 1949, 1959, 1969, 1979, 1989)
refined = c(32, 31.2, 27, 21, 14.9, 8.8)
not_refined = c(16.3, 23.1, 23.6, 27.7, 34.6, 33.9)
sugar_df = data.frame(period, refined, not_refined)
sugar_df

(a)

Пред се, ја плотираме потрошувачката на шеќер низ времето. Може да се забележи дека иако податоците не се конкретно линеарни, има одредена линеарна зависност.

with(sugar_df, plot(period, refined))

with(sugar_df, plot(period, not_refined))

Креираме два едноставни линеарни модели, ги печатиме информациите за нив (оценети коефициенти, стандардни грешки, p-вредности, R2) и ги плотираме одредените прави на регресијата со податоците:

refined_lm = lm(refined ~ period, data=sugar_df)
not_refined_lm = lm(not_refined ~ period, data=sugar_df)

summary(refined_lm)

Call:
lm(formula = refined ~ period, data = sugar_df)

Residuals:
      1       2       3       4       5       6 
-2.6905  1.3924  2.0752  0.9581 -0.2590 -1.4762 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 981.47648   95.78387   10.25 0.000511 ***
period       -0.48829    0.04877  -10.01 0.000559 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.04 on 4 degrees of freedom
Multiple R-squared:  0.9616,    Adjusted R-squared:  0.952 
F-statistic: 100.2 on 1 and 4 DF,  p-value: 0.0005593

summary(not_refined_lm)

Call:
lm(formula = not_refined ~ period, data = sugar_df)

Residuals:
      1       2       3       4       5       6 
-1.1905  1.9924 -1.1248 -0.6419  2.6410 -1.6762 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -683.87352   96.35750  -7.097  0.00208 **
period         0.36171    0.04906   7.373  0.00180 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.052 on 4 degrees of freedom
Multiple R-squared:  0.9315,    Adjusted R-squared:  0.9143 
F-statistic: 54.36 on 1 and 4 DF,  p-value: 0.001804

with(sugar_df, plot(period, refined))
abline(refined_lm$coefficients[1], refined_lm$coefficients[2])

with(sugar_df, plot(period, not_refined))
abline(not_refined_lm$coefficients[1], not_refined_lm$coefficients[2])

Годишните интервали на доверба се:

confint(refined_lm, "period")

            2.5 %     97.5 %
period -0.6236872 -0.3528842

confint(not_refined_lm, "period")

           2.5 %    97.5 %
period 0.2255019 0.4979267

(b)

Пред се, ги пресметуваме просеците и ги плотираме:

sugar_df["average"] = (sugar_df$refined + sugar_df$not_refined) / 2
sugar_df

with(sugar_df, plot(period, average))

Креираме едноставен модел за просекот, за да тестираме дали коефициентот пред period е 0.

average_lm = lm(average ~ period, data=sugar_df)
summary(average_lm)

Call:
lm(formula = average ~ period, data = sugar_df)

Residuals:
      1       2       3       4       5       6 
-1.9405  1.6924  0.4752  0.1581  1.1910 -1.5762 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 148.80148   77.08724   1.930    0.126
period       -0.06329    0.03925  -1.612    0.182

Residual standard error: 1.642 on 4 degrees of freedom
Multiple R-squared:  0.3939,    Adjusted R-squared:  0.2424 
F-statistic:   2.6 on 1 and 4 DF,  p-value: 0.1822

with(sugar_df, plot(period, average))
abline(average_lm$coefficients[1], average_lm$coefficients[2])

Коефициентот пред period е близок до 0, но можеме да ја орфрлиме хипотезата дека не е 0, бидејќи p-вредноста е голема.

Задача 2

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

K	Yield
0	1753.9
15	3107.7
10	2400.0
40	4923.1
30	4415.4
5	2861.6
50	5246.2
50	4938.4
40	3723.0
5	3184.6
5	3046.2
30	4892.3
10	3538.5
0	2553.8
40	4784.6
30	4000.0
10	3323.1
20	3184.6
15	4184.6
40	4461.5
0	2723.1
40	4692.3
20	4215.4
50	4784.6
20	3600.0
40	4153.9
15	3169.3

Со следниов код ги подготвуваме податоците во R.

K = c(0, 15, 10, 40, 30, 5, 50, 50, 40, 5, 5, 30, 10, 0, 40, 30, 10, 20, 15, 40, 0, 40, 20, 50, 20, 40, 15)
yield = c(1753.9, 3107.7, 2400, 4923.1, 4415.4, 2861.6, 5246.2, 4938.4, 3723, 3184.6, 3046.2, 4892.3, 3538.5, 2553.8, 4784.6, 4000, 3323.1, 3184.6, 4184.6, 4461.5, 2723.1, 4692.3, 4215.4, 4784.6, 3600, 4153.9, 3169.3)
grass_df = data.frame(K, yield)
grass_df

(a) (b)

Креираме модели, со yield кренато на секој степен од 2 до 10, со пресметан логаритам со сите основи од 2 до 10 и со коренуван yield, и за секој од нив го правиме и линеарниот модел што го има и самото К. Од нив, го бираме моделот со највисок R2.

r2 = c()
model_name = c()

generate_model <- function(predictor, dataframe) {
  model = lm(as.formula(paste("yield ~", predictor)), data=dataframe)
  r2 <<- append(r2, summary(model)$r.squared)
  model_name <<- append(model_name, predictor)
  
  model = lm(as.formula(paste("yield ~", predictor, " + K")), data=dataframe)
  r2 <<- append(r2, summary(model)$r.squared)
  model_name <<- append(model_name, paste(predictor, "+K", sep=""))
}

generate_model("K", grass_df)

grass_df["K_sqrt"] = grass_df$K ^ 0.5
generate_model("K_sqrt", grass_df)

for (exp in 2:10) {
  grass_df[paste("K_", exp, sep="")] = grass_df$K ^ exp
  generate_model(paste("K_", exp, sep=""), grass_df)
}

for (base in 2:10) {
  grass_df[paste("K_log_", base, sep="")] = log(grass_df$K + 0.01, base=base)
  generate_model(paste("K_log_", base, sep=""), grass_df)
}

best_idx = which.max(r2)
best_model_name = model_name[best_idx]
best_model = lm(as.formula(paste("yield ~", best_model_name)), data=grass_df)
best_predictors = strsplit(best_model_name, "\\+")
print(best_predictors)

[[1]]
[1] "K_sqrt" "K"     

summary(best_model)

Call:
lm(formula = as.formula(paste("yield ~", best_model_name)), data = grass_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-838.96 -241.29   92.13  281.15  731.22 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2351.06     240.58   9.772 7.69e-10 ***
K_sqrt        206.93     135.73   1.525    0.140    
K              22.55      17.67   1.277    0.214    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 433.8 on 24 degrees of freedom
Multiple R-squared:  0.7961,    Adjusted R-squared:  0.7791 
F-statistic: 46.85 on 2 and 24 DF,  p-value: 5.171e-09

(c)

Ги плотираме резидуалите за најдобриот модел. Делува дека тие не следат некоја шема.

grass_df["residuals"] = best_model$residuals
plot(grass_df[c("K", "residuals")])

Задача 3

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

Carbohydrate	Age	Weight	Protein
33	33	100	14
40	47	92	15
37	49	135	18
27	35	144	12
30	46	140	15
43	52	101	15
34	62	95	14
48	23	101	17
30	32	98	15
38	42	105	14
50	31	108	17
51	61	85	19
30	63	130	19
36	40	127	20
41	50	109	15
42	64	107	16
46	56	117	18
24	61	100	13
35	48	118	18
37	28	102	14

Со следниов код ги подготвуваме податоците во R.

carbo = c(33, 40, 37, 27, 30, 43, 34, 48, 30, 38, 50, 51, 30, 36, 41, 42, 46, 24, 35, 37)
age = c(33, 47, 49, 35, 46, 52, 62, 23, 32, 42, 31, 61, 63, 40, 50, 64, 56, 61, 48, 28)
weight = c(100, 92, 135, 144, 140, 101, 95, 101, 98, 105, 108, 85, 130, 127, 109, 107, 117, 100, 118, 102)
protein = c(14, 15, 18, 12, 15, 15, 14, 17, 15, 14, 17, 19, 19, 20, 15, 16, 18, 13, 18, 14)
carbo_df = data.frame(carbo, age, weight, protein)
carbo_df

(a)

for (predictor in c("age", "weight", "protein")) {
  plot(carbo_df[c(predictor, "carbo")])
}

Возраста делува линеарно поврзана, но има и доста outliers.

Кај тежината можеме да забележиме хетероскедастичност (помали тежини делуваат ко да имаат повисока варијанса).

Протеините делуваат линеарно поврзани, но не може да се потврди само визуелно.

(b)

Како што расте вистинската вредност, расте и грешката, што сугерира дека постои линеарна зависност помеѓу carbo и резидуалот.

carbo_all_lm = lm(carbo ~ age + weight + protein, data=carbo_df)
summary(carbo_all_lm)

Call:
lm(formula = carbo ~ age + weight + protein, data = carbo_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.3424  -4.8203   0.9897   3.8553   7.9087 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 36.96006   13.07128   2.828  0.01213 * 
age         -0.11368    0.10933  -1.040  0.31389   
weight      -0.22802    0.08329  -2.738  0.01460 * 
protein      1.95771    0.63489   3.084  0.00712 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.956 on 16 degrees of freedom
Multiple R-squared:  0.4805,    Adjusted R-squared:  0.3831 
F-statistic: 4.934 on 3 and 16 DF,  p-value: 0.01297

carbo_df["residuals_all_lm"] = carbo_all_lm$residuals
plot(carbo_df[c("carbo", "residuals_all_lm")])

(c)

Со пресметување на F статистиката, согледуваме дека вредноста не е доволно различна од центрираната F дистрибуција за да се отфрли \(H_0: \beta_1 = 0\).

carbo_age_protein_lm = lm(carbo ~ age + protein, data=carbo_df)
age_protein_anova = anova(carbo_age_protein_lm)
carbo_protein_lm = lm(carbo ~ protein, data=carbo_df)
protein_anova = anova(carbo_protein_lm)

dfs = c(age_protein_anova$Df[3], protein_anova$Df[2])
Ds = c(age_protein_anova$`Sum Sq`[3], protein_anova$`Sum Sq`[2])
delta_D = Ds[2] - Ds[1]
delta_dfs = dfs[2] - dfs[1]
F_statistic = (delta_D / delta_dfs) / (Ds[1] / dfs[1])
real_F_statistic = pf(F_statistic, delta_dfs, dfs[1])
c(F_statistic, real_F_statistic, F_statistic - real_F_statistic)

[1]  0.511473248  0.515791133 -0.004317885

Задача 4

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

CHOL	Age	BMI
5.94	52	20.7
6.48	65	26.3
4.71	46	21.3
8.83	76	22.7
5.86	51	25.4
5.10	47	21.5
6.52	44	22.7
5.81	43	20.7
6.80	70	23.9
4.65	30	18.9
5.23	33	24.3
6.82	58	23.9
4.97	21	22.2
6.28	78	24.3
8.78	63	26.2
5.15	49	23.8
5.13	56	23.3
2.92	36	19.6
6.74	54	29.2
9.27	67	24.3
5.95	44	22.7
5.57	42	22.0
5.83	71	21.9
4.92	29	22.5
5.74	39	22.4
6.72	33	24.1
4.92	58	20.2
5.57	42	22.7
6.69	58	24.4
6.25	66	27.3

Со следниов код ги подготвуваме во R:

chol = c(5.94, 6.48, 4.71, 8.83, 5.86, 5.1, 6.52, 5.81, 6.8, 4.65, 5.23, 6.82, 4.97, 6.28, 8.78, 5.15, 5.13, 2.92, 6.74, 9.27, 5.95, 5.57, 5.83, 4.92, 5.74, 6.72, 4.92, 5.57, 6.69, 6.25)
age = c(52, 65, 46, 76, 51, 47, 44, 43, 70, 30, 33, 58, 21, 78, 63, 49, 56, 36, 54, 67, 44, 42, 71, 29, 39, 33, 58, 42, 58, 66)
bmi = c(20.7, 26.3, 21.3, 22.7, 25.4, 21.5, 22.7, 20.7, 23.9, 18.9, 24.3, 23.9, 22.2, 24.3, 26.2, 23.8, 23.3, 19.6, 29.2, 24.3, 22.7, 22, 21.9, 22.5, 22.4, 24.1, 20.2, 22.7, 24.4, 27.3)
chol_df = data.frame(chol, age, bmi)
chol_df

Креираме модел со и без bmi и пресметуваме F статистика за да ја тестираме хипотезата \(H_0:\beta_{BMI}=0\). Високото отстапување од F дистрибуцијата ни укажува дека нултата хипотеза се отфрла, т.е. холестеролот зависи од bmi кога е вклучено age.

with_bmi_lm = lm(chol ~ age + bmi, data=chol_df)
with_bmi_anova = anova(with_bmi_lm)
no_bmi_lm = lm(chol ~ age)
no_bmi_anova = anova(no_bmi_lm)

dfs = c(with_bmi_anova$Df[3], no_bmi_anova$Df[2])
Ds = c(with_bmi_anova$`Sum Sq`[3], no_bmi_anova$`Sum Sq`[2])
delta_D = Ds[2] - Ds[1]
delta_dfs = dfs[2] - dfs[1]
F_statistic = (delta_D / delta_dfs) / (Ds[1] / dfs[1])
real_F_statistic = pf(F_statistic, delta_dfs, dfs[1])
c(F_statistic, real_F_statistic, F_statistic - real_F_statistic)

[1] 5.1473853 0.9685125 4.1788728

Задача 5

Податоците од задачата ги сетираме во R. Мораме да внимаваме да дефинираме фактор променлива.

hyper = c(2.3, 4.1, 4.2, 4.0, 4.6, 4.6, 3.8, 5.2, 3.1, 3.7, 3.8)
non_hyper = c(3.0, 4.1, 3.9, 3.1, 3.3, 2.9, 3.3, 3.9)
control = c(3.0, 2.6, 3.1, 2.2, 2.1, 2.4, 2.8, 3.4, 2.9, 2.6, 3.1, 3.2)

plasma = c(hyper, non_hyper, control)
group = factor(c(
  rep("hyper", length(hyper)), 
  rep("non_hyper", length(non_hyper)), 
  rep("control", length(control))
))

plasma_df = data.frame(plasma, group)
plasma_df

(a)

Користиме анова за да пресметаме F статистика. Ниската вредност на Pr(>F) ни укажува дека средините на групите не се еднакви.

plasma_lm = lm(plasma~group)
anova(plasma_lm)

Analysis of Variance Table

Response: plasma
          Df Sum Sq Mean Sq F value    Pr(>F)    
group      2 7.8083  3.9041  11.651 0.0002082 ***
Residuals 28 9.3827  0.3351                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(b)

Прво пресметуваме стандардна девијација низ групите.

pooled_std = sqrt(sum(c(
  var(hyper) * (length(hyper) - 1), 
  var(non_hyper) * (length(non_hyper) - 1),
  var(control) * (length(control) - 1)
)) / (length(plasma) - 3))

Интервалот на доверба е:

estimated_diff = mean(hyper) - mean(non_hyper)
t_alpha_half = qt(0.975, 28)
CI = c(estimated_diff - t_alpha_half * pooled_std, estimated_diff + t_alpha_half * pooled_std)
CI

[1] -0.6778168  1.6937258

(c)

Резидуалите линеарно растат во секоја група.

plasma_df["residuals"] = plasma_lm$residuals
plot(hyper, plasma_df[plasma_df$group == "hyper",]$residuals, ylab="residuals")

plot(non_hyper, plasma_df[plasma_df$group == "non_hyper",]$residuals, ylab="residuals")

plot(control, plasma_df[plasma_df$group == "control",]$residuals, ylab="residuals")

Задача 6

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

weight = c(35.7, 38.4, 34.9, 37.1, 37.1, 37.2, 34.3, 35.5, 36.7, 38.1, 34.5, 36.5, 37.7, 36.9, 33.7, 36, 35.3, 37.2, 36.2, 33.8, 34.7, 36.9, 32, 35.8, 35.2, 38.5, 35.2, 32.9, 34.6, 36.4, 33.5, 35.7, 36.4, 37.8, 32.9, 38, 35.2, 36.1, 33.3, 36.1)
worker = factor(rep(c(1, 2, 3, 4), 10))
day = factor(c(rep(1, 20), rep(2, 20)))

weight_df = data.frame(weight, worker, day)
weight_df

Од двостраната анова заклучуваме дека има значајни разлики помеѓу работниците и помеѓу деновите (F статистика значајно отскокнува од вредноста на дистрибуцијата, обележано и со ѕвездички). Нема доволно показатели за постоечка интеракција.

weight_lm = lm(weight ~ worker + day + worker * day, data=weight_df)
anova(weight_lm)

Analysis of Variance Table

Response: weight
           Df Sum Sq Mean Sq F value    Pr(>F)    
worker      3 54.622 18.2073 14.4948 3.895e-06 ***
day         1  6.084  6.0840  4.8435   0.03508 *  
worker:day  3  2.958  0.9860  0.7850   0.51117    
Residuals  32 40.196  1.2561                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Задача 7

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

samples = c(6.8, 6.6, 5.3, 6.1, 7.5, 7.4, 7.2, 6.5, 7.8, 9.1, 8.8, 9.1)
b = factor(rep(c(-1, -1, 1, 1), 3))
a = factor(c(rep(-1, 4), rep(1, 4), rep(0, 4)))

samples_df = data.frame(samples, a, b)
samples_df

(a)

Ова е матрицата X:

X = matrix(
  c(
    1, −1, −1, −1, 1, 1,
    1, −1, −1, −1, 1, 1, 
    1, −1, −1, 1, −1, −1, 
    1, −1, −1, 1, −1, −1, 
    1, 1, 0, −1, −1, 0, 
    1, 1, 0, −1, −1, 0, 
    1, 1, 0, 1, 1, 0, 
    1, 1, 0, 1, 1, 0, 
    1, 0, 1, −1, 0, −1, 
    1, 0, 1, −1, 0, −1, 
    1, 0, 1, 1, 0, 1, 
    1, 0, 1, 1, 0, 1
  ),
  nrow=12,
  ncol=6,
  byrow=TRUE
)
X

      [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1   -1   -1   -1    1    1
 [2,]    1   -1   -1   -1    1    1
 [3,]    1   -1   -1    1   -1   -1
 [4,]    1   -1   -1    1   -1   -1
 [5,]    1    1    0   -1   -1    0
 [6,]    1    1    0   -1   -1    0
 [7,]    1    1    0    1    1    0
 [8,]    1    1    0    1    1    0
 [9,]    1    0    1   -1    0   -1
[10,]    1    0    1   -1    0   -1
[11,]    1    0    1    1    0    1
[12,]    1    0    1    1    0    1

Пресметуваме XtX. Производот е блок дијагонална матрица (3х3 блоковите по дијагонала се ненулти, другите се нулти).

XtX = t(X) %*% X
XtX

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   12    0    0    0    0    0
[2,]    0    8    4    0    0    0
[3,]    0    4    8    0    0    0
[4,]    0    0    0   12    0    0
[5,]    0    0    0    0    8    4
[6,]    0    0    0    0    4    8

Ги пресметуваме моделите 6.9, 6.10 и 6.12 со оваа спецификација и можеме да се увериме дека навистина се исти како во табелата во текстот.

samples_glm_a_b_ab = glm(samples ~ a * b, family=gaussian, data=samples_df)
samples_glm_a_b = glm(samples ~ a + b, family=gaussian, data=samples_df)
samples_glm_b = glm(samples ~ b, family=gaussian, data=samples_df)

anova(samples_glm_a_b_ab)

Analysis of Deviance Table

Model: gaussian, link: identity

Response: samples

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev
NULL                    11    15.8300
a     2  12.7400         9     3.0900
b     1   0.4033         8     2.6867
a:b   2   1.2067         6     1.4800

anova(samples_glm_a_b)

Analysis of Deviance Table

Model: gaussian, link: identity

Response: samples

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev
NULL                    11    15.8300
a     2  12.7400         9     3.0900
b     1   0.4033         8     2.6867

anova(samples_glm_b)

Analysis of Deviance Table

Model: gaussian, link: identity

Response: samples

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev
NULL                    11     15.830
b     1  0.40333        10     15.427

(b)

Доволно е да ги видиме коефициентите на моделите и да заклучиме дека оценетиот просек ќе е ист.

samples_glm_a_b_ab$coefficients

(Intercept)          a0          a1          b1       a0:b1       a1:b1 
       6.70        1.75        0.75       -1.00        1.50        0.40 

samples_glm_a_b$coefficients

(Intercept)          a0          a1          b1 
  6.3833333   2.5000000   0.9500000  -0.3666667 

Задача 8

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

samples = c(5, 3, 4, 6, 4, 4, 3, 7, 6, 8)
a = factor(c(1, 1, 1, 2, 2, 2, 2, 3, 3, 3))
b = factor(c(1, 2, 2, 1, 1, 2, 2, 1, 2, 2))

samples_df = data.frame(samples, a, b)
samples_df

(a)

За да провериме дали има интеракција, треба да гo тестираме сатурираниот модел. Oд тестот следи дека не можеме да ја отфрлиме хипотезата дека има интеракција.

saturated_lm = glm(samples ~ a * b, family=gaussian, data=samples_df)
anova(saturated_lm, test="Chisq")

Analysis of Deviance Table

Model: gaussian, link: identity

Response: samples

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
NULL                     9    26.0000            
a     2  17.2500         7     8.7500 0.001008 **
b     1   2.6786         6     6.0714 0.143235   
a:b   2   1.0714         4     5.0000 0.651439   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(b)

За да провериме дали има ефект поради А, треба да гo тестираме адитивниот модел. Oд тестот следи дека можеме да ја орфрлиме хипотезата дека има ефект поради А.

additive_lm = glm(samples ~ a + b, family=gaussian, data=samples_df)
anova(additive_lm, test="Chisq")

Analysis of Deviance Table

Model: gaussian, link: identity

Response: samples

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                     9    26.0000              
a     2  17.2500         7     8.7500 0.0001987 ***
b     1   2.6786         6     6.0714 0.1037417    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1