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, plot
jest identyfikatorem wykresu, x
jest położeniem xi położeniem y
y każdego wykresu z wykresem 1 wyśrodkowanym na 0, 0. level
jest poziomem leczenia i response
jest 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
- Jak obliczyć macierz kowariancji i uwzględnić ją w regresji?
- Bloki są bardzo różne i występują silne interakcje bloków leczenia *. Czy należy je analizować osobno?
r
spatial
linear-model
covariance
David LeBauer
źródło
źródło
Odpowiedzi:
1) Możesz modelować korelację przestrzenną z
nlme
biblioteką; 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.
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.
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.
2) Możesz także spróbować włączyć
x
iy
bezpoś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.
Aby porównać wszystkie trzy modele, musisz dopasować je wszystkie
gls
i metodę największej wiarygodności zamiast domyślnej metody REML.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
m2b
nie musi być uważany za najlepszy jak dotąd.Uwaga: Obliczenia te wykonano po zmianie wartości x wykresu 37 na 0.
źródło
model
zamiastm0
, 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,m2
jest dużym przegranym (AIC = 64) wydaje się, że twoja częśćm0
,m1
im2
jak 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ę.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.
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.
źródło