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

10

Фон

У меня есть данные полевого исследования, в котором есть четыре уровня обработки и шесть повторностей в каждом из двух блоков. (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
Графики 37 и 39 оба при x = 18, y = 10. Опечатка?
Аарон оставил переполнение стека
@ Аарон, спасибо, что указал на это. у = [0,10]. Исправлена.
Дэвид Лебауэр

Ответы:

10

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.

Аарон оставил переполнение стека
источник
спасибо за ваш полезный ответ; Непонятно, почему в части 2 вы обновились, modelа не m0, например. m2 <- update(m0, .~.+x*y)так что все три модели можно сравнить с помощью anova(m0,m1,m2); после этого m2большой неудачник (AIC = 64) кажется, что ваша роль
David LeBauer
ps последняя строка должна быть «после изменения значения y графика 37 на 5»; фактическое значение равно 0, но результаты эквивалентны.
Дэвид Лебауэр
Если вы сравниваете m0, m1и, m2как вы предлагаете, вы получите предупреждение: Fitted objects with different fixed effects. REML comparisons are not meaningful. Для сравнения фиксированных эффектов вы должны использовать обычную максимальную вероятность вместо REML. Смотрите редактировать.
Аарон оставил переполнение стека
спасибо за вашу помощь Я не уверен, почему, но я получаю ошибки, когда пытаюсь подогнать другие корреляционные структуры, например, используя corExp, как в примере с Пинейро и Бейтсом. Я открыл вопрос на SO об этой ошибке, но ваш вклад был бы оценен.
Дэвид Лебауэр
4

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 метров, я бы сказал, что совершенно уместно анализировать их отдельно.

sesqu
источник
Визуализация полезна, но если вы сравните нижний правый с верхним правым числами, мне кажется, что местоположение имеет более сильный эффект, чем уровень. (ps я думаю, что l [i] = ответ должен быть r [i] = ...)
Дэвид Лебауэр
Да. Эффект локации замечательный, и вам действительно нужна хорошая модель для этого, прежде чем вы начнете оценивать эффекты лечения. К сожалению, существует много пропущенных данных, поэтому трудно сказать, что это должно быть - лучшее, что я могу придумать, - это эффект моделирования местоположения как среднее значение ответа соседей + случайный компонент, а затем попробовать эту обработку , Это очень подозрительно, поэтому любые дополнительные знания предметной области были бы полезны. опечатка исправлена.
Sesqu
@sesqu ... нет пропущенных данных, есть данные со всех 24 случайно расположенных участков.
Дэвид Лебауэр
Данные отсутствуют в том смысле, что не каждая пара x, y содержит данные.
Аарон оставил переполнение стека