Мастицкий, глава 6
ЛИНЕЙНЫЕ МОДЕЛИ В ДИСПЕРСИОННОМ АНАЛИЗЕ
Формула
Факторы и контрасты
| факторы | is.factor(), as.factor(), уровни (levels), ординальные и номинальные факторы
Usage
factor(x = character(), levels, labels = levels, ordered = is.ordered(x)) Arguments
|
#============================================================================# --- небольшой поясняющий пример
#============================================================================ tmp.s1<-c("a","b","b","a","a","c","b","a","c")
tmp.f1<-factor(tmp.s1)
Str(tmp.f1)
# Factor w/ 3 levels "a","b","c": 1 2 2 1 1 3 2 1 3
tmp.s2<-c("XX","XY","XX","XY","XX","XY","XY","XY","XY")
tmp.f2<-factor(tmp.s2)
Str(tmp.f2)
# Factor w/ 2 levels "XX","XY": 1 2 1 2 1 2 2 2 2
tmp.tab1<-model.matrix(~tmp.f1)
Tmp.tab1
# (Intercept) tmp.fb tmp.fc
#1 1 0 0 <-- a
#2 1 1 0 <-- b
#3 1 1 0 <-- b
tmp.tab2<-model.matrix(~tmp.f2-1)
Tmp.tab2
# tmp.f2XX tmp.f2XY
#1 1 0 <-- XX
#2 0 1 <-- XY
#3 1 0
tmp.tab3<-model.matrix(~tmp.f1+tmp.f2-1)
Tmp.tab3
# tmp.f1a tmp.f1b tmp.f1c tmp.f2XY
#1 1 0 0 0
#2 0 1 0 1
#3 0 1 0 0
#4 1 0 0 1
tmp.tab4<-cbind(model.matrix(~tmp.f1-1),model.matrix(~tmp.f2-1))
Tmp.tab4
# tmp.f1a tmp.f1b tmp.f1c tmp.f2XX tmp.f2XY
#1 1 0 0 1 0
#2 0 1 0 0 1
#3 0 1 0 1 0
Контрасты
Номинальные и ординальные переменные
contr.treatment(n, base = 1, contrasts = TRUE, sparse = FALSE)contr.helmert(n, contrasts = TRUE, sparse = FALSE)contr.poly(n, scores = 1:n, contrasts = TRUE, sparse = FALSE)contr.sum(n, contrasts = TRUE, sparse = FALSE) Виды контрастов| Тип | Описание |
| treatment | Сравнение каждого уровня с базой |
| helmert | Сравнение уровня со средним для предыдущих уровней |
| sum | Сравнение уровня с суммой предыдущих уровней |
| poly | Для ординальных, позволяет использовать нелинейные тренды заданного порядка |
Contr.treatment(4)
2 3 4
1 0 0 0
2 1 0 0
3 0 1 0
4 0 0 1
Contr.helmert(4)
[,1] [,2] [,3]
1 -1 -1 -1
2 1 -1 -1
3 0 2 -1
4 0 0 3
Contr.sum(4)
[,1] [,2] [,3]
1 1 0 0
2 0 1 0
3 0 0 1
4 -1 -1 -1
contr.poly(4)
.L.Q.C
[1,] -0.6708204 0.5 -0.2236068
[2,] -0.2236068 -0.5 0.6708204
[3,] 0.2236068 -0.5 -0.6708204
[4,] 0.6708204 0.5 0.2236068
Параметры
| Параметр | Описание |
| scores | Значения аргумента для вычисления полинома |
| contrasts | Нужно ли создавать реальную матрицу (или фиктивную, т.е. линейно зависимую) |
| sparse | Нужно ли создавать разреженную матрицу |
| base | Какое значение использовать как базу сравнения |
Contr.poly(3)
.L.Q
[1,] -7.071068e-01 0.4082483
[2,] -7.850462e-17 -0.8164966
[3,] 7.071068e-01 0.4082483
contrasts(tmp.f1ord, how.many = 1) <- contr.poly(3)
model.matrix(~tmp.f1ord)
(Intercept) tmp.f1ord.L
1 1 -7.071068e-01
2 1 -7.850462e-17
3 1 -7.850462e-17
4 1 -7.071068e-01
5 1 -7.071068e-01
6 1 7.071068e-01
contrasts(tmp.f1ord, how.many = 2) <- contr.poly(3)
model.matrix(~tmp.f1ord)
(Intercept) tmp.f1ord.L tmp.f1ord.Q
1 1 -7.071068e-01 0.4082483
2 1 -7.850462e-17 -0.8164966
3 1 -7.850462e-17 -0.8164966
4 1 -7.071068e-01 0.4082483
5 1 -7.071068e-01 0.4082483
6 1 7.071068e-01 0.4082483
lm(formula, data, subset, weights, na.action,
method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, offset,...)
model.matrix(~ tmp.f1 + tmp.f2)
(Intercept) tmp.f1b tmp.f1c tmp.f2XY
1 1 0 0 0
2 1 1 0 1
3 1 1 0 0
4 1 0 0 1
5 1 0 0 0
6 1 0 1 1
7 1 1 0 1
8 1 0 0 1
9 1 0 1 1
model.matrix(~ tmp.f1 + tmp.f2,
contrasts = list(tmp.f1 = "contr.sum", tmp.f2 = "contr.helmert"))
(Intercept) tmp.f11 tmp.f12 tmp.f21
1 1 1 0 -1
2 1 0 1 1
3 1 0 1 -1
4 1 1 0 1
5 1 1 0 -1
6 1 -1 -1 1
7 1 0 1 1
8 1 1 0 1
9 1 -1 -1 1
Пример: изучается эффект двух методов лечения, method1, method2, и эффект от их совместного применения method12 по сравнению с суммой частных эффектов. Результат оценивается по показателю y, чем он больше, тем лучше.
Set.seed(123)
beta0 <-5
beta1 <-1
beta2 <-2
beta12<-(-0.7)
beta.effect<-c(0,beta1,beta2,beta1+beta2+beta12)
x.labels<-c("none","method1","method2","method12")
x.val <-x.labels[sample(1:4,100,replace=TRUE)]
x<-factor(x.val)
err<-rnorm(100,0,0.5)
y<-beta0 + beta.effect[x] + err
boxplot(y~as.factor(x))
beta.contrast<-rbind(c(0, 0,0),
C(1, 0,0),
C(0, 1,0),
C(1, 1,1))
colnames(beta.contrast)<-x.labels[-1]
p.lm<-lm(y~x,contrasts=list(x=beta.contrast))
Summary(p.lm)
Call:
lm(formula = y ~ x, contrasts = list(x = beta.contrast))
Residuals:
Min 1Q Median 3Q Max
-1.07775 -0.32030 -0.00897 0.28411 1.04254
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.05113 0.09267 54.509 < 2e-16 ***
xmethod1 0.99490 0.12983 7.663 1.46e-11 ***
xmethod2 1.67092 0.14054 11.890 < 2e-16 ***
xmethod12 -0.40585 0.19049 -2.130 0.0357 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4725 on 96 degrees of freedom
Multiple R-squared: 0.7738, Adjusted R-squared: 0.7668
F-statistic: 109.5 on 3 and 96 DF, p-value: < 2.2e-16
Описание набора данных RIKZ
| Sample | Номер образца |
| C1 | Chaetognatha species Щетинкочелюстны́е, или морски́е стре́лки |
| P1-P25 | Polychaeta species Многощетинковые черви, или полихеты |
| N1 | Nemertini species лентовидные черви |
| CR1-CR28 | Crustacea species Ракообразные |
| M1-M17 | Mollusca species Моллюски |
| I1-I5 | Insecta species |
| week | Неделя Июня (1-4), когда собирались данные |
| angle1 | Угол наклона станции |
| angle2 | Угол наклона для всего пляжа |
| exposure | Индекс, отражающий воздействие волн, длину береговой линии, уклон, мехсостав, и глубину анаэробного слоя (для пляжа) |
| salinity | Средняя соленость (для пляжа) |
| temperature | Температура (для пляжа) |
| NAP | Высота станции наблюдения по отношению к среднему уровню прилива |
| penetrability | Проницаемость |
| grainsize | Размер зерен песка |
| humus | Гумус (органика) |
| chalk | Мел |
| sorting1 | |
| Beach | Пляж |
p<-read.table("RIKZ.txt",header=TRUE,as.is=TRUE,dec=".")
# число видов
p$abundance<-rowSums(p[,2:76])
p1$fBeach<-as.factor(p1$Beach)
p1$lAbundance<-log(1+p1$abundance)
plot(p$abundance,pch=20,cex=2)
plot(p$lAbundance,pch=20,cex=2)
p.lm0<-lm(abundance~1,data=p)
Summary(p.lm0)
p.lm1<-lm(abundance~NAP,data=p)
Summary(p.lm1)
p.lm2<-lm(abundance~fBeach,data=p)
Summary(p.lm2)
p.lm3<-lm(abundance~NAP+fBeach,data=p)
Summary(p.lm3)
p.lm4<-lm(abundance~NAP*fBeach,data=p)
Summary(p.lm4)
p.lm5<-lm(abundance~fBeach/NAP,data=p)
Summary(p.lm5)