Hi,
I tried to account for the spatial autocorrelation in my data, but got error message. Basically, the correspond variable is species richness in each grid and predictor variables are bioclimatic variables, soil fertility and speciation rate of species in each grid. All the predictor variables have been z-transformed. I searched relating posts in this discussion web and tried few ways to incorporate the spatial weight matrix. The example data set has been uploaded, which I sampled 100 rows out from the whole data set. And here is the code I used.
library(dplyr)
example_dt = read.csv("100rows_example_20200829.csv", sep = ",", header = T)
example_dt$grid_no = as.factor(example_dt$grid_no)
# calculate a distance matrix using the coordinates
library(geosphere)
dist = distm(example_dt[,14:15], example_dt[,14:15], fun = distVincentyEllipsoid)
head(dist)
# to name the columns and rows
rownames(dist) = example_dt$grid_no
colnames(dist) = example_dt$grid_no
head(dist)
# using z-transformation to scale the predictors
example_dt$bio1_scaled = (example_dt$bio1-mean(example_dt$bio1))/sd(example_dt$bio1)
example_dt$bio4_scaled = (example_dt$bio4-mean(example_dt$bio4))/sd(example_dt$bio4)
example_dt$bio11_scaled = (example_dt$bio11-mean(example_dt$bio11))/sd(example_dt$bio11)
example_dt$bio12_scaled = (example_dt$bio12-mean(example_dt$bio12))/sd(example_dt$bio12)
example_dt$bio15_scaled = (example_dt$bio15-mean(example_dt$bio15))/sd(example_dt$bio15)
example_dt$cec_scaled = (example_dt$cec-mean(example_dt$cec))/sd(example_dt$cec)
example_dt$pet_scaled = (example_dt$pet-mean(example_dt$pet))/sd(example_dt$pet)
example_dt$spra_scaled = (example_dt$speci_rate-mean(example_dt$speci_rate))/sd(example_dt$speci_rate)
head(example_dt)
library(brms)
brm0 = brm(richness ~ bio1_scaled + bio4_scaled + bio11_scaled + bio12_scaled + bio15_scaled +
cec_scaled + pet_scaled + spra_scaled + (1|grid_no),
data = example_dt, family = poisson(),
cov_ranef = list(grid_no = dist),
chains = 4, iter = 2000, warmup = 500, cores = 8)
##### Error: Within-group covariance matrices must be positive definite.
brm1 = brm(richness ~ bio1_scaled + bio4_scaled + bio11_scaled + bio12_scaled + bio15_scaled +
cec_scaled + pet_scaled + spra_scaled + (1|grid_no),
data = example_dt, family = poisson(),
autocor = cor_car(dist),
chains = 4, iter = 2000, warmup = 500, cores = 8)
##### it ran but warning message said "Converting all non-zero values in âMâ to 1."
brm2 = brm(richness ~ bio1_scaled + bio4_scaled + bio11_scaled + bio12_scaled + bio15_scaled +
cec_scaled + pet_scaled + spra_scaled + (1|gr(grid_no, cov = dist)),
data = example_dt, family = poisson(),
data2 = list(dist = dist),
chains = 4, iter = 2000, warmup = 500, cores = 8)
##### Error: Within-group covariance matrices must be positive definite.
brm3 = brm(richness ~ bio1_scaled + bio4_scaled + bio11_scaled + bio12_scaled + bio15_scaled +
cec_scaled + pet_scaled + spra_scaled + (1|grid_no),
data = example_dt, family = poisson(),
autocor = cor_car(dist),
chains = 4, iter = 2000, warmup = 500, cores = 8)
##### it ran but warning message said "Converting all non-zero values in âMâ to 1."
- Operating System: Windows 10 and the R version is 4.0.1
- brms Version: 2.13.0
I checked that no negative value exists in the âdistâ matrix. The minimum value was 0. And I tried to add an extremely small number (0.00000001) to the âdistâ matrix, but the error message was still the same.
So, how should I correctly produce the spatial weight matrix and incorporate the spatial autocorrelation?
100rows_example_20200829.csv (13.8 KB)