Modelo de ajuste de variograma compatibilidad entre los paquetes geoR, gstat y nlme en R

Question

Estoy tratando de especificar los parámetros de estructura de covarianza en un modelo mixto lineal (usando las facilidades de estructura de correlación en nlme). Los parámetros se estiman a partir de un modelo de variograma ajustado a partir de un modelo gstat o geoR ajustado al semivariograma empírico.

Mi justificación para especificar el modelo de covarianza derivado de gstat en nlme es porque nlme utiliza solo el estimador de semivarianza clásico ajustado al conjunto de datos completo. gstat y geoR permiten el ajuste de modelos de variograma al semivariograma robusto (Cressie) y el establecimiento de una distancia máxima de influencia. Por ejemplo, vea mi pregunta buscando orientación sobre la selección entre el semivariograma clásico y el robusto.

Puedo producir semivariogramas empíricos casi idénticos a los residuos de un modelo solo con intercepto (primera figura a continuación). Asimismo, los ajustes del modelo de variograma exponencial de gstat y geoR son esencialmente idénticos. Sin embargo, ajustar el modelo de variograma exponencial en nlme produce estimaciones de parámetros extrañas.

nlme utiliza una especificación ligeramente diferente de los modelos de variograma y estructuras de correlación en comparación con gstat y geoR: específicamente, nlme estandariza los errores dentro del grupo a una varianza unitaria y utiliza un efecto de pepita multiplicativo en lugar de aditivo (consulte, p. ej., págs. 230-232 en Pinheiro y Bates, 2000, Modelos de efectos mixtos en S y S-Plus). Sin embargo, en ausencia de un efecto de pepita (el ejemplo a continuación), los ajustes exponenciales deberían ser fácilmente traducibles entre los programas.

Aquí hay una ilustración simple de lo que me desconcierta. Este ejemplo utiliza el conjunto de datos de concentraciones de metales pesados a lo largo del río Mosa en los Países Bajos (?meuse en el paquete gstat); la autocorrelación es espacial (las concentraciones de metales pesados más cercanas entre sí en el espacio son más similares en promedio que las concentraciones más alejadas entre sí).

Los variogramas empíricos de los residuos del modelo usando nlme, gstat y geoR son esencialmente idénticos.

# Cargar paquetes necesarios
library(gstat) # ajuste del modelo de variograma
library(sp) # clases y métodos para datos espaciales
library(geoR) # ajuste del modelo de variograma
library(nlme) # modelos lineales mixtos con autocorrelación
library(ggplot2) #gráficos

# Utilizar datos de metales pesados de `gstat`
# Enfoque en concentraciones de zinc
data(meuse)
meuse$dGroup <- 1 # variable de agrupación ficticia para poder utilizar la función `lme`
m1 <- lme(log(zinc) ~ 1 , random = ~ 1|dGroup, data=meuse)

# Variograma en `nlme`
lVar <- Variogram(m1, form = ~ x + y, maxDist=1700, breaks=seq(100, 1700, 100))
#lVar; plot(lVar)

# Variograma en `gstat`
meuse2 <- data.frame(x = meuse$x, y = meuse$y, resid = resid(m1, type='n'))
gmeuse <- meuse2
coordinates(gmeuse) = ~x+y
v <- variogram(resid~1, gmeuse, cutoff=1700, width=100)
# Convertir a clase `variogram` para graficar simultáneamente con el variograma de `nlme`
vdf <- with(v, data.frame(variog = gamma, dist = dist, n.pairs = np)) 
class(vdf) <- c("Variogram", "data.frame")
#vdf; plot(vdf)

# Variograma en `geoR`
geodat <- as.geodata(meuse2)
geodat.v1 <- variog(geodat, breaks = seq(100, 1700, 100), max.dist=1700, option='bin')
# Convertir a clase `variogram` para graficar simultáneamente con el variograma de `nlme`
gVar <- with(geodat.v1, data.frame(variog = v, dist = u, n.pairs = n))
class(gVar) <- c("Variogram", "data.frame")
#gVar; plot(gVar)

# Consolidar variogramas empíricos y graficar
# Variogramas empíricos esencialmente idénticos
# El ajuste es un suavizado LOESS, no un modelo de variograma
empirVar <- data.frame(rbind(gVar, vdf, lVar), 
                       model = rep(c('geoR', 'gstat', 'nlme'), each=17))
theme_update(legend.position = 'right')
empirVarPlot <- ggplot(empirVar, aes(x=dist, y=variog, colour = model)) + 
    geom_point() + geom_smooth(fill='white', alpha=0) + 
    xlab('Distancia') + ylab('Semivariograma') + 
    scale_colour_manual(values=c('red', 'blue', 'green')) + theme_bw()
empirVarPlot

Estimar el ajuste del modelo (por ejemplo, una estructura de covarianza exponencial) también es trivial en todos estos paquetes. geoR y gstat producen ajustes de modelos exponenciales prácticamente idénticos, pero el ajuste del variograma de nlme es completamente diferente tanto en el ajuste como en los valores estimados del semivariograma.

# Estimar modelo de variograma en `geoR`
geoExp <- variofit(geodat.v1, ini=c(1, 300), nugget=0.1, cov.model='exponential')
geoExp

> variofit: parámetros del modelo estimados por WLS (mínimos cuadrados ponderados):
>covariance model is: exponential
>estimaciones de parámetros:
>tausq sigmasq phi
>0.0000 1.2457 338.3446 

# Estimar modelo de variograma en `gstat`
gstExp <- fit.variogram(v, vgm(1, "Exp", 300, 1), fit.method=1)
gstExp

> model psill range
> 1 Nug 0.000000 0.0000
> 2 Exp 1.246994 340.6168

# Estimar modelo de variograma en `nlme`
# Nota: el valor de pepita (segundo valor en el vector) se especifica de manera diferente en `nlme`
# En particular, es 'uno menos la correlación entre dos observaciones tomadas 
# arbitrariamente cercanas' (consulte ?corExp y Pinheiro y Bates 2000)
m2 <- update( m1, corr = corExp(c(300, 0.7), form = ~ x + y, nugget = T) )
m2$modelStruct$corStruct

> Estructura de correlación de clase corExp que representa
> gama pepita
> 1.051978e+08 4.192519e-07

# Gráficarlos todos
plot(geodat.v1, main = 'geoR'); lines(geoExp, col='red')
plot(v, model = gstExp, main = 'gstat')
plot(Variogram(m2, maxDist = 1700, breaks=seq(100, 1700, 100)), main = 'nlme')

Corregir los valores de la correlación de la estructura de nlme para que coincidan con el modelo de gstat parece ser una mejora, pero el ajuste no parece coincidir con las estimaciones del semivariograma. El ajuste exponencial está asintotando a 1, lo cual tiene sentido dado la especificación de las estructuras de correlación en nlme, pero no estoy seguro de cómo modificar este comportamiento o 'traducir' los parámetros de gstat o geoR a términos de nlme.

# Corregir la estructura de correlación para que coincida con el ajuste del modelo de variograma de `gstat`
m3 <- update(m1, corr = corExp(gstExp[2, 3], form = ~ x + y, 
                               nugget=FALSE, fixed=TRUE))
m3$modelStruct$corStruct

> Estructura de correlación de clase corExp que representa
> rango
> 340.6168

plot(Variogram(m3, maxDist = 1700, breaks=seq(100, 1700, 100)), main = 'nlme')

Mis preguntas, para resumir:

1. ¿Por qué existe una gran discrepancia entre los ajustes del modelo de variograma exponencial entre gstat/geoR y nlme?
2. ¿Cómo traducir un ajuste de variograma de gstat/geoR a nlme dadas las diferentes especificaciones de la estructura de correlación y los efectos de pepita?

Answer 1

Veo una salida diferente, a saber:

> m2 <- update( m1, corr = corExp(c(300, 0.7), form = ~ x + y, nugget = T) )
Error in lme.formula(fixed = log(zinc) ~ 1, data = meuse, random = ~1 |  :   
  problema nlminb, código de error de convergencia = 1
  mensaje = falsa convergencia (8)

esto es nlme_3.1-119 en R versión 3.1.2 (2014-10-31) Plataforma: x86_64-pc-linux-gnu (64-bit). ¿Qué arroja tu sessionInfo()?

Parece que nlme asume un correlograma (tope = 1), donde geoR/gstat modelan semivariogramas, donde el tope no está limitado a ser 1. Eso parece explicar el "ajuste insuficiente" de tu última figura.

Answer 2

Comprobé esto y encontré que usando el método ML en lugar del REML predeterminado solucionará el problema de subestimación.