Como posso explicar a covariância espacial em um modelo linear?

10

fundo

Eu tenho dados de um estudo de campo em que existem quatro níveis de tratamento e seis repetições em cada um dos dois blocos. (4x6x2 = 48 observações)

Os blocos estão a cerca de 1,6 km de distância e, dentro dos blocos, há uma grade de parcelas de 42, 2m x 4m e uma passarela de 1m de largura; meu estudo utilizou apenas 24 parcelas em cada bloco.

Eu gostaria de avaliar avaliar covariância espacial.

Aqui está um exemplo de análise usando os dados de um único bloco, sem considerar a covariância espacial. No conjunto de dados, ploté o ID do gráfico, xé o local x, e yo local y de cada gráfico com o gráfico 1 centrado em 0, 0. levelé o nível de tratamento e responsea variável de resposta.

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)

Questões

  1. Como posso calcular uma matriz de covariância e incluí-la na minha regressão?
  2. Os blocos são muito diferentes e existem fortes interações de tratamento *. É apropriado analisá-los separadamente?
David LeBauer
fonte
11
Os gráficos 37 e 39 estão ambos em x = 18, y = 10. Erro de digitação?
Aaron deixou o Stack Overflow em
@ Aaron obrigado por apontar isso. y = [0,10]. Fixo.
David LeBauer

Respostas:

10

1) Você pode modelar correlação espacial com a nlmebiblioteca; existem vários modelos possíveis que você pode escolher. Veja as páginas 260-266 de Pinheiro / Bates.

Um bom primeiro passo é fazer um variograma para ver como a correlação depende da distância.

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

Aqui, o semivariograma amostral aumenta com a distância, indicando que as observações são de fato espacialmente correlacionadas.

Uma opção para a estrutura de correlação é uma estrutura esférica; que podem ser modelados da seguinte maneira.

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

Esse modelo parece se encaixar melhor que o modelo sem estrutura de correlação, embora seja inteiramente possível que ele também possa ser aprimorado com uma das outras estruturas de correlação possíveis.

> 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) Você também pode tentar incluir xe ydiretamente no modelo; isso pode ser apropriado se o padrão de correlação depender mais do que apenas a distância. No seu caso (olhando as fotos de sesqu), parece que, para esse bloco, você pode ter um padrão diagonal.

Aqui, estou atualizando o modelo original em vez de m0, porque estou alterando apenas os efeitos fixos; portanto, os modelos devem ser ajustados usando a máxima probabilidade.

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

Para comparar todos os três modelos, você precisará ajustá-los a todos glse ao método de máxima verossimilhança, em vez do método padrão de 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

Lembre-se de que, especialmente com o seu conhecimento do estudo, você poderá criar um modelo melhor do que qualquer um deles. Ou seja, o modelo m2bnão deve necessariamente ser considerado o melhor ainda.

Nota: Esses cálculos foram realizados após a alteração do valor x do gráfico 37 para 0.

Aaron deixou Stack Overflow
fonte
obrigado pela sua resposta útil; não está claro por que, na parte 2, você atualizou em modelvez de m0, por exemplo. m2 <- update(m0, .~.+x*y)para que todos os três modelos possam ser comparados usando anova(m0,m1,m2); Após fazer isso, m2é um grande perdedor (AIC = 64), parece que a sua parte
David LeBauer
ps a última linha deve ser 'depois de alterar o valor y do gráfico 37 para 5'; o valor real é 0, mas os resultados são equivalentes.
David LeBauer
Se você comparar m0, m1e m2conforme sugerir, receberá o aviso: Fitted objects with different fixed effects. REML comparisons are not meaningful. Para comparar efeitos fixos, use a probabilidade máxima regular em vez de REML. Veja editar.
Aaron deixou o Stack Overflow
Obrigado por toda sua ajuda. Não sei por que, mas estou recebendo erros quando tento ajustar outras estruturas de correlação, por exemplo, usando corExp como no exemplo de Pinheiro e Bates. Eu abriu uma pergunta sobre SO sobre este erro, mas sua entrada seria apreciada.
David LeBauer
4

1) Qual é a sua variável de explicação espacial? Parece que o plano x * y seria um modelo ruim para o efeito espacial.

trama de tratamentos e respostas

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) Visto que os blocos estão separados por uma milha e você espera ver efeitos por meros 30 metros, eu diria que é inteiramente apropriado analisá-los separadamente.

sesqu
fonte
A visualização é útil, mas se você comparar o canto inferior direito ao canto superior direito das figuras, parece-me que o local tem um efeito mais forte que o nível. (ps eu acho l [i] = resposta deve ser r [i] = ...)
David LeBauer
Sim. O efeito da localização é notável e, portanto, você realmente deseja um bom modelo para isso antes de começar a estimar os efeitos do tratamento. Infelizmente, há muitos dados ausentes, por isso é difícil dizer o que deveria ser - o melhor que posso apresentar seria um efeito de localização de modelagem como uma média da resposta dos vizinhos + componente aleatório e, em seguida, tentar o tratamento com base nisso. . Isso é muito suspeito, portanto, qualquer conhecimento adicional sobre o domínio seria valioso. erro de digitação corrigido.
sesqu 15/05
@sesqu ... não há dados ausentes, os dados de todas as 24 parcelas localizadas aleatoriamente estão lá.
David LeBauer
Faltam dados no sentido de que nem todo par x, y possui dados.
Aaron deixou Stack Overflow