MATLAB: “Tilde” model formula of “R” applied to Matlab

MATLAB and Simulink Student Suiteregressorresidualsvariogram

Hi everybody,
I'm facing a problem which bothering me a lot at the moment… I'm trying to get the residuals variogram with Matlab and "BMElib2.0" toolbox from famous "Meuse" data (package "sp" within "R"). The problem is that "lscov" function returns something wrong (I think so) and my residuals are overestimated by Matlab (see attachment "Variogram_residuals"). It goes over 0.5 on y-axis!
However, I should obtain this (see attachment "Minasny"). It's a caption from 2007 Minasny & McBratney article "Spatial prediction of soil properties using EBLUP with the Matérn covariance function" and as you see the residuals don't go over 0.3.
I've tried to get results from "R" and it's the same as Minasny & McBratney (see attachment "Residuals_R"). The blue circles (residuals data) peak at 0.3 on y-axis.
I've discovered that "R" uses a special command named "model formula" which allows to specify a regressor to another variable. It's the "tilde" in below code. I need to make a dependency between log(zinc) and sqrt(dist).
My question is simple: is there a "model formula"-ish command for Matlab in order to get the same residuals as "R" and Minasny & McBratney article? Or maybe my "lscov" function is not suited for this one?
Thanks in advance for any help you'll provide,
Julien Kaba.
#################################
My "R" code (the "tilde" I'm interested in is visible in "v1" command):
library(lattice)
library(gstat)
library(sp)
load(system.file("data", "meuse.rda", package = "sp"))
v1 = variogram(log(zinc) ~ sqrt(dist), locations = ~x + y, data = meuse, width=50, cutoff=2000)
v2 = variogram(log(zinc) ~ 1, locations = ~x + y, data = meuse, width=50, cutoff=2000)
m1 = fit.variogram(v1, vgm(psill = 0.1089, "Mat", range = 40, nugget = 0.084, kappa = 8)) 
m2 = fit.variogram(v2, vgm(psill = 4.9335, "Exp", range = 1412, nugget = 0.086, kappa = 1)) 
plot(gamma~dist, v2, ylim = c(0, 1.05*max(v2$gamma)),col='red', ylab = 
       'semivariance', xlab = 'distance') 
lines(variogramLine(m2, 2000), col='red',ylab ='',xlab='') 
points(gamma~dist, v1, col='blue') 
lines(variogramLine(m1, 2000), col='blue')
My Matlab code for residuals ("BMElib2.0" is used to get variograms):
% "coord" (coordinates of the 155 points) and "z" (log of zinc concentration for this 155 points) 
% are available in attachments (.mat files)
X_Zn = [ones(size(coord(:,2))) coord(:,2) coord(:,1)];
b_derive_Zn = lscov(X_Zn,z);
Zn_derive = X_Zn*b_derive_Zn;
Zn_vec = z - Zn_derive;
figure
plZn_vec = (0:40:2000);
[dZn_vec,vZn_vec,oZn_vec] = vario(coord,Zn_vec,plZn_vec,'kron');
plot(dZn_vec,vZn_vec,'r*');
xlabel('Distance (m)','FontSize',11); ylabel('Variogram','FontSize',11);

Best Answer

So, I'm answering myself: "lscov" is perfectly suited but my Matlab code was wrong. I tried:
X_Zn = [ones(size(coord(:,2))) coord(:,2) coord(:,1)];
but it's not correct. Coordinates don't matter here, it's "sqrt(ndist)" which is important (a linear relation is specified in Minasny & McBratney article, p. 327) ! So, the right code is:
X_Zn = [ones(size(coord(:,2))) sqrt(ndist)];
Running the script gives the correct residuals variogram :-)... Hope it might help other users!
Related Question