Wie kann ich die räumliche Kovarianz in einem linearen Modell berücksichtigen?

10

Hintergrund

Ich habe Daten aus einer Feldstudie, in der es vier Behandlungsstufen und sechs Wiederholungen in jedem von zwei Blöcken gibt. (4x6x2 = 48 Beobachtungen)

Die Blöcke sind ungefähr 1 Meile voneinander entfernt, und innerhalb der Blöcke gibt es ein Raster von 42, 2 mx 4 m großen Parzellen und einen 1 m breiten Gehweg; Meine Studie verwendete nur 24 Parzellen in jedem Block.

Ich möchte die räumliche Kovarianz bewerten.

Hier ist eine Beispielanalyse unter Verwendung der Daten aus einem einzelnen Block ohne Berücksichtigung der räumlichen Kovarianz. Im Datensatz plotist die Plot-ID xdie x-Position und ydie y-Position jedes Plots, wobei Plot 1 auf 0, 0 zentriert ist. Dies levelist die Behandlungsstufe und responsedie Antwortvariable.

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)

Fragen

  1. Wie kann ich eine Kovarianzmatrix berechnen und in meine Regression einbeziehen?
  2. Die Blöcke sind sehr unterschiedlich und es gibt starke Wechselwirkungen zwischen Behandlung und Block. Ist es angebracht, sie separat zu analysieren?
David LeBauer
quelle
1
Die Diagramme 37 und 39 liegen beide bei x = 18, y = 10. Tippfehler?
Aaron verließ Stack Overflow
@ Aaron danke für den Hinweis. y = [0,10]. Fest.
David LeBauer

Antworten:

10

1) Sie können die räumliche Korrelation mit der nlmeBibliothek modellieren . Es gibt mehrere mögliche Modelle, die Sie auswählen können. Siehe Seiten 260-266 von Pinheiro / Bates.

Ein guter erster Schritt besteht darin, ein Variogramm zu erstellen, um zu sehen, wie die Korrelation von der Entfernung abhängt.

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

Hier nimmt das Probensemivariogramm mit der Entfernung zu, was darauf hinweist, dass die Beobachtungen tatsächlich räumlich korreliert sind.

Eine Option für die Korrelationsstruktur ist eine sphärische Struktur; das könnte folgendermaßen modelliert werden.

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

Dieses Modell scheint besser zu passen als das Modell ohne Korrelationsstruktur, obwohl es durchaus möglich ist, dass es auch mit einer der anderen möglichen Korrelationsstrukturen verbessert werden könnte.

> 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) Sie können auch versuchen, xund ydirekt in das Modell aufzunehmen. Dies könnte angemessen sein, wenn das Korrelationsmuster von mehr als nur der Entfernung abhängt. In Ihrem Fall (wenn Sie sich die Bilder von sesqu ansehen) scheint es, dass Sie für diesen Block ohnehin ein diagonales Muster haben.

Hier aktualisiere ich das Originalmodell anstelle von m0, da ich nur die festen Effekte ändere, sodass beide Modelle mit maximaler Wahrscheinlichkeit angepasst werden sollten.

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

Um alle drei Modelle zu vergleichen, müssen Sie sie alle mit glsund der Maximum-Likelihood-Methode anstelle der Standardmethode von REML anpassen.

> 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

Denken Sie daran, dass Sie insbesondere mit Ihrem Wissen über die Studie möglicherweise ein Modell entwickeln können, das besser ist als jedes andere. Das heißt, das Modell m2bsollte noch nicht unbedingt als das beste angesehen werden.

Hinweis: Diese Berechnungen wurden durchgeführt, nachdem der x-Wert von Diagramm 37 auf 0 geändert wurde.

Aaron verließ Stack Overflow
quelle
Vielen Dank für Ihre hilfreiche Antwort. Es ist nicht klar, warum Sie in Teil 2 modelanstelle von m0z. m2 <- update(m0, .~.+x*y)so dass alle drei Modelle mit verglichen werden können anova(m0,m1,m2); nachdem dies getan wurde, m2ist ein großer Verlierer (AIC = 64) es scheint, dass Ihr Teil
David LeBauer
ps die letzte Zeile sollte 'nach Änderung des y-Wertes von Diagramm 37 auf 5' sein; Der tatsächliche Wert ist 0, aber die Ergebnisse sind äquivalent.
David LeBauer
Wenn Sie vergleichen m0, m1und m2wie Sie vorschlagen , dass Sie die Warnung erhalten: Fitted objects with different fixed effects. REML comparisons are not meaningful. Um feste Effekte vergleichen Sie haben Wahrscheinlichkeit regelmäßig maximal zu nutzen , anstatt REML. Siehe Bearbeiten.
Aaron verließ Stack Overflow
danke für all deine Hilfe. Ich bin mir nicht sicher warum, aber ich erhalte Fehler, wenn ich versuche, andere Korrelationsstrukturen anzupassen, z. B. mit corExp wie im Beispiel von Pinheiro und Bates. Ich habe auf SO eine Frage zu diesem Fehler gestellt, aber Ihre Eingabe wäre willkommen.
David LeBauer
4

1) Was ist Ihre räumliche Erklärungsvariable? Es sieht so aus, als wäre die x * y-Ebene ein schlechtes Modell für den räumlichen Effekt.

Darstellung der Behandlungen und Reaktionen

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) Angesichts der Tatsache, dass die Blöcke 1 Meile voneinander entfernt sind und Sie Effekte für nur 30 Meter erwarten, würde ich sagen, dass es völlig angemessen ist, sie separat zu analysieren.

sesqu
quelle
Die Visualisierung ist hilfreich, aber wenn Sie die untere rechte mit der oberen rechten Seite der Figuren vergleichen, scheint es mir, dass die Position einen stärkeren Effekt als die Ebene hat. (ps Ich denke, l [i] = Antwort sollte r [i] = ... sein)
David LeBauer
Ja. Der Standorteffekt ist bemerkenswert, und deshalb möchten Sie wirklich ein gutes Modell dafür, bevor Sie mit der Abschätzung der Behandlungseffekte beginnen. Leider fehlen viele Daten, so dass es schwierig ist zu sagen, was es sein soll. Das Beste, was ich mir einfallen lassen kann, ist ein Modellierung des Standorteffekts als Durchschnitt der Antwort des Nachbarn + der zufälligen Komponente, und dann wird die Behandlung darauf versucht . Das ist sehr verdächtig, daher wäre jedes zusätzliche Domänenwissen wertvoll. Tippfehler behoben.
Sesqu
@sesqu ... es fehlen keine Daten, Daten aus allen 24 zufällig angeordneten Plots sind vorhanden.
David LeBauer
Es fehlen Daten in dem Sinne, dass nicht jedes x, y-Paar Daten hat.
Aaron verließ Stack Overflow