Задание структуры корреляции для линейной смешанной модели с использованием пакета рампы в R

Question

Я пытаюсь создать линейную смешанную модель (lmm), которая позволяет пространственную корреляцию между точками (имеет длину lat/long для каждой точки). Я хотел бы, чтобы пространственная корреляция основывалась на большом круговом расстоянии между точками.

В комплект пакета ramps входит корреляционная структура, которая вычисляет расстояние «haversine» - хотя у меня возникают проблемы с ее реализацией. Я ранее использовал другие корреляционные структуры (corGaus, corExp) и не испытывал никаких затруднений. Я предполагаю, что corRGaus с метрикой 'haversine' может быть реализован точно так же.

Я могу успешно создать lmm с пространственной корреляцией, рассчитанной на плоском расстоянии, используя функцию lme.

Я также могу создать линейную модель (не смешанную) с пространственной корреляцией, рассчитанную с использованием большого кругового расстояния, хотя есть ошибки в структуре корреляции с использованием команды gls.

При попытке использования команды gls для линейной модели с большой круговой дистанции у меня есть следующие ошибки:

x = runif(20, 1,50) 
y = runif(20, 1,50) 
gls(x ~ y, cor = corRGaus(form = ~ x + y)) 

Generalized least squares fit by REML 
Model: x ~ y 
Data: NULL 
Log-restricted-likelihood: -78.44925 

Coefficients: 
(Intercept)   y 
24.762656602 0.007822469 

Correlation Structure: corRGaus 
Formula: ~x + y 
Parameter estimate(s): 
Error in attr(object, "fixed") && unconstrained : 
invalid 'x' type in 'x && y' 

Когда я увеличить размер данных Есть ошибки распределения памяти (по-прежнему очень небольшой набор данных):

x = runif(100, 1, 50) 
y = runif(100, 1, 50) 
lat = runif(100, -90, 90) 
long = runif(100, -180, 180) 
gls(x ~ y, cor = corRGaus(form = ~ x + y)) 

Error in glsEstimate(glsSt, control = glsEstControl) : 
'Calloc' could not allocate memory (18446744073709551616 of 8 bytes) 

При попытке запуска смешанной модели с использованием lme команды и corRGaus из ramps пакета следующие результаты:

x = runif(100, 1, 50) 
y = runif(100, 1, 50) 
LC = c(rep(1, 50) , rep(2, 50)) 
lat = runif(100, -90, 90) 
long = runif(100, -180, 180) 

lme(x ~ y,random = ~ y|LC, cor = corRGaus(form = ~ long + lat)) 

Error in `coef<-.corSpatial`(`*tmp*`, value = value[parMap[, i]]) : 
    NA/NaN/Inf in foreign function call (arg 1) 
In addition: Warning messages: 
1: In nlminb(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt, lmePars), : 
    NA/NaN function evaluation 
2: In nlminb(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt, lmePars), : 
    NA/NaN function evaluation 

Я не уверен, как это сделать. Функция «haversine» - это то, что я хочу использовать для завершения моих моделей, но у меня возникли проблемы с их реализацией. Вопросов о пакете ramps очень мало, и я видел очень мало реализаций. Любая помощь будет принята с благодарностью.

Я ранее пытался изменить пакет nlme и не смог это сделать. Я задал вопрос о this, где мне рекомендовали использовать пакет ramps.

Я использую R 3.0.0 на компьютере под управлением Windows 8.

Answer 1

ОК, здесь есть опция, которая реализует различные пространственные корреляционные структуры в gls/nlme с расстоянием до Хаверсина.

Различные классы классов corSpatial уже имеют механизмы для создания корреляционной матрицы из пространственных ковариаций с учетом метрики расстояния.К сожалению, dist не реализует расстояние Хаверсина, а dist - это функция, называемая corSpatial для вычисления матрицы расстояния от пространственных ковариаций.

Расчет расстояний матрицы выполняется в getCovariate.corSpatial. Измененная форма этого метода пройдет соответствующее расстояние до других методов, и большинство методов не нужно будет изменять.

Здесь я создаю новый класс corStruct, corHaversine и изменять только getCovariate и один другой метод (Dim), который определяет, какая функция корреляции используется. Те методы, которые не нуждаются в модификации, копируются из эквивалентных методов corSpatial. Аргумент (новый) mimic в corHaversine берет имя пространственного класса с интересующей функцией корреляции: по умолчанию он установлен на «corSpher».

Предостережение: за исключением того, что этот код работает для сферических и гауссовых корреляционных функций, я на самом деле не проверил много проверок.

#### corHaversine - spatial correlation with haversine distance 

# Calculates the geodesic distance between two points specified by radian latitude/longitude using Haversine formula. 
# output in km 
haversine <- function(x0, x1, y0, y1) { 
    a <- sin((y1 - y0)/2)^2 + cos(y0) * cos(y1) * sin((x1 - x0)/2)^2 
    v <- 2 * asin(min(1, sqrt(a))) 
    6371 * v 
} 

# function to compute geodesic haversine distance given two-column matrix of longitude/latitude 
# input is assumed in form decimal degrees if radians = F 
# note fields::rdist.earth is more efficient 
haversineDist <- function(xy, radians = F) { 
    if (ncol(xy) > 2) stop("Input must have two columns (longitude and latitude)") 
    if (radians == F) xy <- xy * pi/180 
    hMat <- matrix(NA, ncol = nrow(xy), nrow = nrow(xy)) 
    for (i in 1:nrow(xy)) { 
     for (j in i:nrow(xy)) { 
      hMat[j,i] <- haversine(xy[i,1], xy[j,1], xy[i,2], xy[j,2]) 
      } 
     } 
    as.dist(hMat) 
} 

## for most methods, machinery from corSpatial will work without modification 
Initialize.corHaversine <- nlme:::Initialize.corSpatial 
recalc.corHaversine <- nlme:::recalc.corSpatial 
Variogram.corHaversine <- nlme:::Variogram.corSpatial 
corFactor.corHaversine <- nlme:::corFactor.corSpatial 
corMatrix.corHaversine <- nlme:::corMatrix.corSpatial 
coef.corHaversine <- nlme:::coef.corSpatial 
"coef<-.corHaversine" <- nlme:::"coef<-.corSpatial" 

## Constructor for the corHaversine class 
corHaversine <- function(value = numeric(0), form = ~ 1, mimic = "corSpher", nugget = FALSE, fixed = FALSE) { 
    spClass <- "corHaversine" 
    attr(value, "formula") <- form 
    attr(value, "nugget") <- nugget 
    attr(value, "fixed") <- fixed 
    attr(value, "function") <- mimic 
    class(value) <- c(spClass, "corStruct") 
    value 
} # end corHaversine class 
environment(corHaversine) <- asNamespace("nlme") 

Dim.corHaversine <- function(object, groups, ...) { 
    if (missing(groups)) return(attr(object, "Dim")) 
    val <- Dim.corStruct(object, groups) 
    val[["start"]] <- c(0, cumsum(val[["len"]] * (val[["len"]] - 1)/2)[-val[["M"]]]) 
    ## will use third component of Dim list for spClass 
    names(val)[3] <- "spClass" 
    val[[3]] <- match(attr(object, "function"), c("corSpher", "corExp", "corGaus", "corLin", "corRatio"), 0) 
    val 
} 
environment(Dim.corHaversine) <- asNamespace("nlme") 


## getCovariate method for corHaversine class 
getCovariate.corHaversine <- function(object, form = formula(object), data) { 
    if (is.null(covar <- attr(object, "covariate"))) {   # if object lacks covariate attribute 
     if (missing(data)) {         # if object lacks data 
      stop("need data to calculate covariate") 
      } 
     covForm <- getCovariateFormula(form) 
     if (length(all.vars(covForm)) > 0) {     # if covariate present 
      if (attr(terms(covForm), "intercept") == 1) {  # if formula includes intercept 
       covForm <- eval(parse(text = paste("~", deparse(covForm[[2]]),"-1",sep=""))) # remove intercept 
       } 
      # can only take covariates with correct names 
      if (length(all.vars(covForm)) > 2) stop("corHaversine can only take two covariates, 'lon' and 'lat'") 
      if (!all(all.vars(covForm) %in% c("lon", "lat"))) stop("covariates must be named 'lon' and 'lat'") 
      covar <- as.data.frame(unclass(model.matrix(covForm, model.frame(covForm, data, drop.unused.levels = TRUE)))) 
      covar <- covar[,order(colnames(covar), decreasing = T)] # order as lon ... lat 
      } 
     else { 
      covar <- NULL 
      } 

     if (!is.null(getGroupsFormula(form))) {     # if groups in formula extract covar by groups 
      grps <- getGroups(object, data = data) 
      if (is.null(covar)) { 
       covar <- lapply(split(grps, grps), function(x) as.vector(dist(1:length(x))))  # filler? 
       } 
      else { 
       giveDist <- function(el) { 
        el <- as.matrix(el) 
        if (nrow(el) > 1) as.vector(haversineDist(el))      
        else numeric(0) 
        } 
       covar <- lapply(split(covar, grps), giveDist) 
       } 
      covar <- covar[sapply(covar, length) > 0] # no 1-obs groups 
      } 
     else {         # if no groups in formula extract distance 
      if (is.null(covar)) { 
       covar <- as.vector(dist(1:nrow(data))) 
       } 
      else { 
       covar <- as.vector(haversineDist(as.matrix(covar))) 
       } 
      } 
     if (any(unlist(covar) == 0)) {   # check that no distances are zero 
      stop("cannot have zero distances in \"corHaversine\"") 
      } 
     } 
    covar 
    } # end method getCovariate 
environment(getCovariate.corHaversine) <- asNamespace("nlme") 

Чтобы проверить, что это работает, данный параметр диапазон 1000:

## test that corHaversine runs with spherical correlation (not testing that it WORKS ...) 
library(MASS) 
set.seed(1001) 
sample_data <- data.frame(lon = -121:-22, lat = -50:49) 
ran <- 1000 # 'range' parameter for spherical correlation 
dist_matrix <- as.matrix(haversineDist(sample_data)) # haversine distance matrix 
# set up correlation matrix of response 
corr_matrix <- 1-1.5*(dist_matrix/ran)+0.5*(dist_matrix/ran)^3 
corr_matrix[dist_matrix > ran] = 0 
diag(corr_matrix) <- 1 
# set up covariance matrix of response 
sigma <- 2 # residual standard deviation 
cov_matrix <- (diag(100)*sigma) %*% corr_matrix %*% (diag(100)*sigma) # correlated response 
# generate response 
sample_data$y <- mvrnorm(1, mu = rep(0, 100), Sigma = cov_matrix) 

# fit model 
gls_haversine <- gls(y ~ 1, correlation = corHaversine(form=~lon+lat, mimic="corSpher"), data = sample_data) 
summary(gls_haversine) 

# Correlation Structure: corHaversine 
# Formula: ~lon + lat 
# Parameter estimate(s): 
# range 
# 1426.818 
# 
# Coefficients: 
#     Value Std.Error t-value p-value 
# (Intercept) 0.9397666 0.7471089 1.257871 0.2114 
# 
# Standardized residuals: 
#  Min   Q1  Med   Q3  Max 
# -2.1467696 -0.4140958 0.1376988 0.5484481 1.9240042 
# 
# Residual standard error: 2.735971 
# Degrees of freedom: 100 total; 99 residual 

тестирование, что он работает с гауссовой корреляции с параметром диапазон = 100:

## test that corHaversine runs with Gaussian correlation 
ran = 100 # parameter for Gaussian correlation 
corr_matrix_gauss <- exp(-(dist_matrix/ran)^2) 
diag(corr_matrix_gauss) <- 1 
# set up covariance matrix of response 
cov_matrix_gauss <- (diag(100)*sigma) %*% corr_matrix_gauss %*% (diag(100)*sigma) # correlated response 
# generate response 
sample_data$y_gauss <- mvrnorm(1, mu = rep(0, 100), Sigma = cov_matrix_gauss) 

# fit model 
gls_haversine_gauss <- gls(y_gauss ~ 1, correlation = corHaversine(form=~lon+lat, mimic = "corGaus"), data = sample_data) 
summary(gls_haversine_gauss) 

С lme:

## runs with lme 
# set up data with group effects 
group_y <- as.vector(sapply(1:5, function(.) mvrnorm(1, mu = rep(0, 100), Sigma = cov_matrix_gauss))) 
group_effect <- rep(-2:2, each = 100) 
group_y = group_y + group_effect 
group_name <- factor(group_effect) 
lme_dat <- data.frame(y = group_y, group = group_name, lon = sample_data$lon, lat = sample_data$lat) 
# fit model 
lme_haversine <- lme(y ~ 1, random = ~ 1|group, correlation = corHaversine(form=~lon+lat, mimic = "corGaus"), data = lme_dat, control=lmeControl(opt = "optim")) 
summary(lme_haversine) 

# Correlation Structure: corHaversine 
# Formula: ~lon + lat | group 
# Parameter estimate(s): 
# range 
# 106.3482 
# Fixed effects: y ~ 1 
#     Value Std.Error DF  t-value p-value 
# (Intercept) -0.0161861 0.6861328 495 -0.02359033 0.9812 
# 
# Standardized Within-Group Residuals: 
#  Min   Q1  Med   Q3  Max 
# -3.0393708 -0.6469423 0.0348155 0.7132133 2.5921573 
# 
# Number of Observations: 500 
# Number of Groups: 5 

Answer 2

Смотрите, если этот ответ на R-Help полезно: http://markmail.org/search/?q=list%3Aorg.r-project.r-help+winsemius+haversine#query:list%3Aorg.r-project.r-help%20winsemius%20haversine+page:1+mid:ugecbw3jjwphu2pb+state:results

Я только что проверил и и не появляется, что ramps или nlme пакеты были изменены, чтобы включить эти изменения, предложенные Malcolm Fairbrother, так что вы будете необходимо сделать некоторые взлома. Я не хочу, чтобы меня рассматривали за щедрость, так как я не отправляю протестированное решение, и я тоже не мечтал об этом.