Jak mogę uwzględnić kowariancję przestrzenną w modelu liniowym?

10

tło

Mam dane z badań terenowych, w których istnieją cztery poziomy leczenia i sześć powtórzeń w każdym z dwóch bloków. (4x6x2 = 48 obserwacji)

Bloki są oddalone od siebie o około 1 milę, aw obrębie bloków znajduje się siatka o powierzchni 42, 2m x 4m i chodnik o szerokości 1m; moje badanie wykorzystało tylko 24 wykresy w każdym bloku.

Chciałbym ocenić ocenę kowariancji przestrzennej.

Oto przykładowa analiza wykorzystująca dane z jednego bloku, bez uwzględnienia kowariancji przestrzennej. W zestawie danych, plotjest identyfikatorem wykresu, xjest położeniem xi położeniem yy każdego wykresu z wykresem 1 wyśrodkowanym na 0, 0. leveljest poziomem leczenia i responsejest zmienną odpowiedzi.

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)

pytania

  1. Jak obliczyć macierz kowariancji i uwzględnić ją w regresji?
  2. Bloki są bardzo różne i występują silne interakcje bloków leczenia *. Czy należy je analizować osobno?
David LeBauer
źródło
1
Wykresy 37 i 39 mają wartości x = 18, y = 10. Literówka?
Aaron opuścił Stack Overflow
@Aaron dziękuję za zwrócenie na to uwagi. y = [0,10]. Naprawiony.
David LeBauer

Odpowiedzi:

10

1) Możesz modelować korelację przestrzenną z nlmebiblioteką; istnieje kilka możliwych modeli, które możesz wybrać. Zobacz strony 260-266 Pinheiro / Bates.

Dobrym pierwszym krokiem jest wykonanie wariogramu, aby zobaczyć, jak korelacja zależy od odległości.

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

W tym przypadku przykładowy semiwariogram rośnie wraz z odległością, wskazując, że obserwacje są rzeczywiście skorelowane przestrzennie.

Jedną opcją dla struktury korelacji jest struktura sferyczna; które można modelować w następujący sposób.

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

Ten model wydaje się pasować lepiej niż model bez struktury korelacji, chociaż jest całkowicie możliwe, że można go również ulepszyć za pomocą jednej z innych możliwych struktur korelacji.

> 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) Możesz także spróbować włączyć xi ybezpośrednio do modelu; może to być właściwe, jeśli wzór korelacji zależy od więcej niż odległości. W twoim przypadku (patrząc na zdjęcia sesqu) wydaje się, że i tak dla tego bloku możesz mieć ukośny wzór.

Tutaj aktualizuję oryginalny model zamiast m0, ponieważ zmieniam tylko stałe efekty, więc oba modele powinny być dopasowane z maksymalnym prawdopodobieństwem.

> 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 **

Aby porównać wszystkie trzy modele, musisz dopasować je wszystkie glsi metodę największej wiarygodności zamiast domyślnej metody 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

Pamiętaj, że zwłaszcza dzięki swojej wiedzy na temat badania możesz być w stanie wymyślić model lepszy niż którykolwiek z nich. Oznacza to, że model m2bnie musi być uważany za najlepszy jak dotąd.

Uwaga: Obliczenia te wykonano po zmianie wartości x wykresu 37 na 0.

Aaron opuścił Stack Overflow
źródło
dziękuję za pomocną odpowiedź; nie jest jasne, dlaczego w części 2 zaktualizowałeś modelzamiast m0, np. m2 <- update(m0, .~.+x*y)dzięki czemu można porównać wszystkie trzy modele za pomocą anova(m0,m1,m2); po zrobieniu tego, m2jest dużym przegranym (AIC = 64) wydaje się, że twoja część
David LeBauer
ps ostatni wiersz powinien być „po zmianie wartości y wykresu 37 na 5”; rzeczywista wartość wynosi 0, ale wyniki są równoważne.
David LeBauer
Jeśli porównasz m0, m1i m2jak sugerujesz otrzymasz ostrzeżenie: Fitted objects with different fixed effects. REML comparisons are not meaningful. Aby porównać ustalone efekty, musisz użyć regularnego maksymalnego prawdopodobieństwa zamiast REML. Zobacz edycję.
Aaron opuścił Stack Overflow
Dziękuję Ci za całą Twoją pomoc. Nie jestem pewien dlaczego, ale dostaję błędy, gdy próbuję dopasować inne struktury korelacji, np. Używając corExp jak w przykładzie Pinheiro i Batesa. Mam otwarte pytanie na SO o tym błędzie, ale Twój wkład będzie mile widziane.
David LeBauer
4

1) Jaka jest twoja zmienna wyjaśniająca przestrzennie? Wygląda na to, że płaszczyzna x * y byłaby złym modelem dla efektu przestrzennego.

wykres leczenia i odpowiedzi

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) Biorąc pod uwagę, że bloki są oddalone o 1 milę od siebie, a spodziewasz się zobaczyć efekty na zaledwie 30 metrach, powiedziałbym, że całkowicie właściwe jest analizowanie ich osobno.

sesqu
źródło
Wizualizacja jest pomocna, ale jeśli porównasz prawy dolny róg z prawym górnym rysunkiem, wydaje mi się, że lokalizacja ma silniejszy efekt niż poziom. (ps Myślę, że l [i] = odpowiedź powinna być r [i] = ...)
David LeBauer
Tak. Efekt lokalizacji jest niezwykły, dlatego naprawdę potrzebujesz dobrego modelu, zanim zaczniesz szacować efekty leczenia. Niestety, brakuje wielu danych, więc trudno powiedzieć, co powinno być - najlepszym, co mogę wymyślić, byłby efekt lokalizacji modelowania jako średnia odpowiedź sąsiadów + losowy składnik, a następnie wypróbowanie leczenia na tym . To bardzo podejrzane, więc każda dodatkowa wiedza na temat domen byłaby cenna. literówka naprawiona.
sesqu
@sesqu ... nie brakuje danych, są dane ze wszystkich 24 losowo rozmieszczonych wykresów.
David LeBauer
Brakuje danych w tym sensie, że nie każda para x, y ma dane.
Aaron opuścił Stack Overflow