使用R中的渐变包为线性混合模型指定相关结构

距离矩阵的计算是在 getCovariate.corSpatial 中执行的.此方法的修改形式将与其他方法保持适当的距离，并且大多数方法都不需要修改.在这里，我创建一个新的 corStruct 类 corHaversine ，并仅修改 getCovariate 和其他方法( Dim )，确定使用哪个相关函数.这些不需要修改的方法是从等效的 corSpatial 方法中复制的. corHaversine 中的(新) mimic 参数采用具有相关功能的空间类的名称:默认情况下，它设置为" corSpher ".注意:除了确保该代码针对球形和高斯相关函数运行之外，我还没有做过很多检查. #### corHaversine-与Haversine距离的空间相关性#使用Haversine公式计算由弧度纬度/经度指定的两点之间的测地距离.#输出(以公里为单位)Haversine<-function(x0，x1，y0，y1){a<-sin((y1-y0)/2)^ 2 + cos(y0)* cos(y1)* sin((x1-x0)/2)^ 2v<-2 * asin(min(1，sqrt(a)))6371 * v}#函数在给定经度/纬度的两列矩阵的情况下计算测地线的hasversine距离如果弧度= F，则以十进制形式假定#输入#注意fields :: rdist.earth效率更高haversineDist<-函数(xy，弧度= F){如果(ncol(xy)> 2)stop(输入必须具有两列(经度和纬度)")如果(弧度== F)xy 1)as.vector(haversineDist(el))其他数值(0)}covar<-lapply(split(covar，grps)，GiveDist)}covar<-covar [sapply(covar，length)>0]#没有1-obs组}否则{#如果公式中没有组提取距离如果(is.null(covar)){covar<-as.vector(dist(1:nrow(data)))}别的 {covar<-as.vector(haversineDist(as.matrix(covar)))}}if(any(unlist(covar)== 0)){#检查距离是否为零停止(在\"corHaversine \"中不能有零距离)}}科瓦#end方法getCovariate环境(getCovariate.corHaversine)<-asNamespace("nlme") 要测试它是否运行，请给定范围参数1000: ##测试corHaversine以球形相关性运行(不测试它是否起作用...)图书馆(MASS)set.seed(1001)sample_data<-data.frame(lon = -121:-22，lat = -50:49)ran<-1000#'范围'参数用于球形相关dist_matrix<-as.matrix(haversineDist(sample_data))#Haversine距离矩阵#建立响应的相关矩阵corr_matrix<-1-1.5 *(dist_matrix/ran)+ 0.5 *(dist_matrix/ran)^ 3corr_matrix [dist_matrix>跑] = 0diag(corr_matrix)<-1#建立响应协方差矩阵sigma<-2#残留标准偏差cov_matrix<-(diag(100)* sigma)％*％corr_matrix％*％(diag(100)* sigma)#相关响应#产生回应sample_data $ y<-mvrnorm(1，mu = rep(0，100)，Sigma = cov_matrix)#适合型号gls_haversine<-gls(y〜1，相关= corHaversine(form =〜lon + lat，mimic ="corSpher")，data = sample_data)摘要(gls_haversine)#关联结构:corHaversine#公式:〜lon + lat#参数估算值:# 范围#1426.818##系数:#值标准错误t值p值#(拦截)0.9397666 0.7471089 1.257871 0.2114##标准化残差:#最小Q1中Q3最大#-2.1467696 -0.4140958 0.1376988 0.5484481 1.9240042##残留标准误:2.739971#自由度:总共100；残留99 使用范围参数= 100的高斯相关性进行测试: ##测试corHaversine以高斯相关性运行ran = 100#高斯相关参数corr_matrix_gauss<-exp(-(dist_matrix/ran)^ 2)diag(corr_matrix_gauss)<-1#建立响应协方差矩阵cov_matrix_gauss<--(diag(100)* sigma)％*％corr_matrix_gauss％*％(diag(100)* sigma)#相关响应#产生回应sample_data $ y_gauss<-mvrnorm(1，mu = rep(0，100)，Sigma = cov_matrix_gauss)#适合型号gls_haversine_gauss<-gls(y_gauss〜1，相关= corHaversine(form =〜lon + lat，mimic ="corGaus")，data = sample_data)摘要(gls_haversine_gauss) 使用 lme : ##与lme一起运行#设置具有组效果的数据group_y<-as.vector(sapply(1:5，function(.)mvrnorm(1，mu = rep(0，100)，Sigma = cov_matrix_gauss)))group_effect<-rep(-2:2，每个= 100)group_y = group_y + group_effectgroup_name<-factor(group_effect)lme_dat<-data.frame(y = group_y，group = group_name，lon = sample_data $ lon，lat = sample_data $ lat)#适合型号lme_haversine<-lme(y〜1，随机=〜1 | group，相关= corHaversine(form =〜lon + lat，mimic ="corGaus")，data = lme_dat，control = lmeControl(opt ="optim")))摘要(lme_haversine)#关联结构:corHaversine#公式:〜lon + lat |团体#参数估算值:# 范围#106.3482#固定效果:y〜1#值标准错误DF t值p值#(拦截)-0.0161861 0.6861328 495 -0.02359033 0.9812##标准化组内残差:#最小Q1中Q3最大#-3.0393708 -0.6469423 0.0348155 0.7132133 2.5921573##观察数:500#团体人数:5 I am trying to create a linear mixed model (lmm) that allows for a spatial correlation between points (have lat/long for each point). I would like the spatial correlation to be based upon the great circular distance between points.The package ramps includes a correlation structure that computes the ‘haversine’ distance – although I am having trouble implementing it. I have previously used other correlation structures (corGaus, corExp) and not had any difficulties. I am assuming the corRGaus with the 'haversine' metric can be implemented in the same way.I am able to successfully create an lmm with spatial correlation calculated on a planar distance using the lme function.I am also able to create a linear model (not mixed) with spatial correlation calculated using great circular distance although there are errors with the correlation structure using the gls command.When trying to the use the gls command for a linear model with the great circular distance I have the following errors: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: NULLLog-restricted-likelihood: -78.44925Coefficients: (Intercept) y24.762656602 0.007822469Correlation Structure: corRGaus Formula: ~x + y Parameter estimate(s):Error in attr(object, "fixed") && unconstrained : invalid 'x' type in 'x && y'When I increase the size of the data there are memory allocation errors (still a very small dataset):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)When trying to run a mixed model using the lme command and the corRGaus from the ramps package the following results: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 evaluation2: In nlminb(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt, lmePars), : NA/NaN function evaluationI am unsure about how to proceed with this method. The "haversine" function is what I want to use to complete my models, but I am having trouble implementing them. There are very few questions anywhere about the ramps package, and I have seen very few implementations. Any helps would be greatly appreciated.I have previously attempted to modify the nlme package and was unable to do so. I posted a question about this, where I was recommended to use the ramps package.I am using R 3.0.0 on a Windows 8 computer. 解决方案 OK, here is an option that implements various spatial correlation structures in gls/nlme with haversine distance.The various corSpatial-type classes already have machinery in place to construct a correlation matrix from spatial covariates, given a distance metric. Unfortunately, dist does not implement haversine distance, and dist is the function called by corSpatial to compute a distance matrix from the spatial covariates.The distance matrix computations are performed in getCovariate.corSpatial. A modified form of this method will pass the appropriate distance to other methods, and the majority of methods will not need to be modified.Here, I create a new corStruct class, corHaversine, and modify only getCovariate and one other method (Dim) that determines which correlation function is used. Those methods which do not need modification, are copied from equivalent corSpatial methods. The (new) mimic argument in corHaversine takes the name of the spatial class with the correlation function of interest: by default, it is set to "corSpher".Caveat: beyond ensuring that this code runs for spherical and Gaussian correlation functions, I haven't really done a lot of checking.#### corHaversine - spatial correlation with haversine distance# Calculates the geodesic distance between two points specified by radian latitude/longitude using Haversine formula.# output in kmhaversine <- 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 efficienthaversineDist <- 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 modificationInitialize.corHaversine <- nlme:::Initialize.corSpatialrecalc.corHaversine <- nlme:::recalc.corSpatialVariogram.corHaversine <- nlme:::Variogram.corSpatialcorFactor.corHaversine <- nlme:::corFactor.corSpatialcorMatrix.corHaversine <- nlme:::corMatrix.corSpatialcoef.corHaversine <- nlme:::coef.corSpatial"coef<-.corHaversine" <- nlme:::"coef<-.corSpatial"## Constructor for the corHaversine classcorHaversine <- 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 classenvironment(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 classgetCovariate.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 getCovariateenvironment(getCovariate.corHaversine) <- asNamespace("nlme")To test that this runs, given range parameter of 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 correlationdist_matrix <- as.matrix(haversineDist(sample_data)) # haversine distance matrix# set up correlation matrix of responsecorr_matrix <- 1-1.5*(dist_matrix/ran)+0.5*(dist_matrix/ran)^3corr_matrix[dist_matrix > ran] = 0diag(corr_matrix) <- 1# set up covariance matrix of responsesigma <- 2 # residual standard deviationcov_matrix <- (diag(100)*sigma) %*% corr_matrix %*% (diag(100)*sigma) # correlated response# generate responsesample_data$y <- mvrnorm(1, mu = rep(0, 100), Sigma = cov_matrix)# fit modelgls_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 residualTesting that it runs with Gaussian correlation, with range parameter = 100:## test that corHaversine runs with Gaussian correlationran = 100 # parameter for Gaussian correlationcorr_matrix_gauss <- exp(-(dist_matrix/ran)^2)diag(corr_matrix_gauss) <- 1# set up covariance matrix of responsecov_matrix_gauss <- (diag(100)*sigma) %*% corr_matrix_gauss %*% (diag(100)*sigma) # correlated response# generate responsesample_data$y_gauss <- mvrnorm(1, mu = rep(0, 100), Sigma = cov_matrix_gauss)# fit modelgls_haversine_gauss <- gls(y_gauss ~ 1, correlation = corHaversine(form=~lon+lat, mimic = "corGaus"), data = sample_data)summary(gls_haversine_gauss)With lme:## runs with lme# set up data with group effectsgroup_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_effectgroup_name <- factor(group_effect)lme_dat <- data.frame(y = group_y, group = group_name, lon = sample_data$lon, lat = sample_data$lat)# fit modellme_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 这篇关于使用R中的渐变包为线性混合模型指定相关结构的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！ 上岸，阿里云！