StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
primarykey
Id
19013386
data
AcceptedAnswerId
0
AnswerCount
0
ClosedDate
CommentCount
4
CommunityOwnedDate
CreationDate
2013-09-25T19:18:53.590
FavoriteCount
0
LastActivityDate
2013-09-26T17:21:16.127
LastEditDate
2013-09-26T17:21:16.127
LastEditorUserId
2711730
OwnerUserId
2711730
ParentId
18857443
PostTypeId
2
Score
4
ViewCount
0
LastEditorDisplayName
text
Body
OK, here is an option that implements various spatial correlation structures in <code>gls</code>/<code>nlme</code> with haversine distance. The various <code>corSpatial</code>-type classes already have machinery in place to construct a correlation matrix from spatial covariates, given a distance metric. Unfortunately, <code>dist</code> does not implement haversine distance, and <code>dist</code> is the function called by <code>corSpatial</code> to compute a distance matrix from the spatial covariates. The distance matrix computations are performed in <code>getCovariate.corSpatial</code>. 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 <code>corStruct</code> class, <code>corHaversine</code>, and modify only <code>getCovariate</code> and one other method (<code>Dim</code>) that determines which correlation function is used. Those methods which do not need modification, are copied from equivalent <code>corSpatial</code> methods. The (new) <code>mimic</code> argument in <code>corHaversine</code> takes the name of the spatial class with the correlation function of interest: by default, it is set to "<code>corSpher</code>". Caveat: beyond ensuring that this code runs for spherical and Gaussian correlation functions, I haven't really done a lot of checking. <pre><code>#### 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") </code></pre> To test that this runs, given range parameter of 1000: <pre><code>## 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 </code></pre> Testing that it runs with Gaussian correlation, with range parameter = 100: <pre><code>## 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) </code></pre> With <code>lme</code>: <pre><code>## 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 </code></pre>
Tags
Title
singulars
PostAcceptedAnswerId
1. This table or related slice is empty.
PostParentId
1. POSpecifying a correlation structure for a linear mixed model using the ramps package in R
 singulars
 PostTypePostTypeId
 PTQuestion
PostTypePostTypeId
1. PTAnswer
UserLastEditorUserId
1. USNate Pope
UserOwnerUserId
1. USNate Pope
plurals
PostLinksPostIdRelatedPostId
1. This table or related slice is empty.
PostLinksRelatedPostIdPostId
1. This table or related slice is empty.
PostsAcceptedAnswerId
1. POSpecifying a correlation structure for a linear mixed model using the ramps package in R
 singulars
 PostTypePostTypeId
 PTQuestion
PostsParentIdCreationDate
1. This table or related slice is empty.
VotesPostIdCreationDate
1. VO
 singulars
 PostPostId
 PO
 UserUserId
 This table or related slice is empty.
 VoteTypeVoteTypeId
 VTUpMod
2. VO
 singulars
 PostPostId
 PO
 UserUserId
 This table or related slice is empty.
 VoteTypeVoteTypeId
 VTUpMod
3. VO
 singulars
 PostPostId
 PO
 UserUserId
 This table or related slice is empty.
 VoteTypeVoteTypeId
 VTAcceptedByOriginator
CommentsPostId

Querying!

Guidance

A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.