Как я могу объяснить пространственную ковариацию в линейной модели?

Фон

У меня есть данные полевого исследования, в котором есть четыре уровня обработки и шесть повторностей в каждом из двух блоков. (4x6x2 = 48 наблюдений)

Блоки находятся примерно в 1 миле друг от друга, а внутри блоков есть сетка из 42, 2 х 4 м участков и дорожка шириной 1 м; мое исследование использовало только 24 участка в каждом блоке.

Я хотел бы оценить оценку пространственной ковариации.

Вот пример анализа с использованием данных из одного блока без учета пространственной ковариации. В наборе данных plotэто идентификатор графика, xэто местоположение x и местоположение yy каждого графика с графиком 1, центрированным на 0, 0. levelэто уровень обработки и responseпеременная отклика.

layout <- structure(list(plot = c(1L, 3L, 5L, 7L, 8L, 11L, 12L, 15L, 16L, 
17L, 18L, 22L, 23L, 26L, 28L, 30L, 31L, 32L, 35L, 36L, 37L, 39L, 
40L, 42L), level = c(0L, 10L, 1L, 4L, 10L, 0L, 4L, 10L, 0L, 4L, 
0L, 1L, 0L, 10L, 1L, 10L, 4L, 4L, 1L, 1L, 1L, 0L, 10L, 4L), response = c(5.93, 
5.16, 5.42, 5.11, 5.46, 5.44, 5.78, 5.44, 5.15, 5.16, 5.17, 5.82, 
5.75, 4.48, 5.25, 5.49, 4.74, 4.09, 5.93, 5.91, 5.15, 4.5, 4.82, 
5.84), x = c(0, 0, 0, 3, 3, 3, 3, 6, 6, 6, 6, 9, 9, 12, 12, 12, 
15, 15, 15, 15, 18, 18, 18, 18), y = c(0, 10, 20, 0, 5, 20, 25, 
10, 15, 20, 25, 15, 20, 0, 15, 25, 0, 5, 20, 25, 0, 10, 20, 
25)), .Names = c("plot", "level", "response", "x", "y"), row.names = c(NA, 
-24L), class = "data.frame")

model <- lm(response ~ level, data = layout)      
summary(model)

Вопросов

1. Как я могу вычислить ковариационную матрицу и включить ее в свою регрессию?
2. Блоки сильно отличаются друг от друга, и есть сильные взаимодействия с блоком лечения. Уместно ли их анализировать отдельно?

1) Вы можете моделировать пространственную корреляцию с nlmeбиблиотекой; Есть несколько возможных моделей, которые вы можете выбрать. Смотрите страницы 260-266 Пиньейру / Бейтса.

Хороший первый шаг - сделать вариограмму, чтобы увидеть, как корреляция зависит от расстояния.

library(nlme)
m0 <- gls(response ~ level, data = layout)  
plot(Variogram(m0, form=~x+y))

Здесь выборочная вариограмма увеличивается с расстоянием, указывая, что наблюдения действительно пространственно коррелированы.

Одним из вариантов корреляционной структуры является сферическая структура; это можно смоделировать следующим образом.

m1 <- update(m0, corr=corSpher(c(15, 0.25), form=~x+y, nugget=TRUE))

Эта модель, кажется, подходит лучше, чем модель без корреляционной структуры, хотя вполне возможно, что ее также можно улучшить с помощью одной из других возможных корреляционных структур.

> anova(m0, m1)
   Model df     AIC      BIC    logLik   Test  L.Ratio p-value
m0     1  3 46.5297 49.80283 -20.26485                        
m1     2  5 43.3244 48.77961 -16.66220 1 vs 2 7.205301  0.0273

2) Вы также можете попробовать включить xи yнепосредственно в модель; это может быть целесообразно, если картина корреляции зависит не только от расстояния. В вашем случае (глядя на картинки sesqu) кажется, что для этого блока в любом случае вы можете иметь диагональный рисунок.

Здесь я обновляю исходную модель вместо m0, потому что я изменяю только фиксированные эффекты, поэтому модели должны подходить с максимальной вероятностью.

> model2 <- update(model, .~.+x*y)
> anova(model, model2)
Analysis of Variance Table

Model 1: response ~ level
Model 2: response ~ level + x + y + x:y
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1     22 5.3809                                
2     19 2.7268  3    2.6541 6.1646 0.004168 **

Чтобы сравнить все три модели, вам нужно согласовать их с glsметодом максимального правдоподобия вместо метода по умолчанию REML.

> m0b <- update(m0, method="ML")
> m1b <- update(m1, method="ML")
> m2b <- update(m0b, .~x*y)
> anova(m0b, m1b, m2b, test=FALSE)
    Model df      AIC      BIC     logLik
m0b     1  3 38.22422 41.75838 -16.112112
m1b     2  5 35.88922 41.77949 -12.944610
m2b     3  5 29.09821 34.98847  -9.549103

Помните, что, особенно учитывая ваши знания об исследовании, вы сможете найти модель, которая лучше, чем любая из них. То есть модель еще m2bне обязательно должна считаться лучшей.

Примечание. Эти расчеты были выполнены после изменения значения x графика 37 на 0.

1) Какова ваша пространственная объясняющая переменная? Похоже, плоскость x * y была бы плохой моделью для пространственного эффекта.

i=c(1,3,5,7,8,11,14,15,16,17,18,22,23,25,28,30,31,32,35,36,39,39,41,42)
l=rep(NA,42)[i];l[i]=level
r=rep(NA,42)[i];r[i]=response
image(t(matrix(-l,6)));title("treatment")
image(t(matrix(-r,6)));title("response")

2) Учитывая, что блоки находятся на расстоянии 1 мили друг от друга, и вы ожидаете увидеть эффекты всего лишь для 30 метров, я бы сказал, что совершенно уместно анализировать их отдельно.