Dear Dr. Pebesma,
I appreciate your help very much and I agree with you!
My head is clearer now but I still have some questions.
My first question is if I set fit.variogram to fit.sills=FALSE, fit.ranges=FALSE, to I get computation of SSErr... (I guess so).
That been said lets move forward :
Based on Geostatistic Modeling Spatial uncertainty (Chiles) page 109, in
the weighted least square method, the weight are also squared :
$sum{w^2 * [?(h_j)-?(h_j:b)]^2}$
So this agrees with what you have suggested me (and which I have applied). THANK GOD
Now, based on the gstat user's manual page 42 table 4.2 values for fit.
# 1 = $N_j$
# 2 = $N_j/{?(h_j)}^2$
# 6 = no weights (OLS) Ordinary Least Square
# 7 = $N_j/h_j^2$
I have made this code... that computes the R2 based of the fitted method used.
I have also include an example of the data... and my code.
If you tell me the code is correct, I will move forward, but I still do not comprehend why I cannot obtain similar values.
I must be missing something... any clue...???
I how the thing about comparing results with other software is a risky business but it is statistics, so I would at least be happy with simial result but look what I get... and tell me this is not mind puzzling?
Anyway I hope this code is good:
findR2 <- function(vario, fit, fit.method=c(1,2,6,7)){
# fit = fit.variogram output
# vario = variogram output
if(!is(vario, "gstatVariogram")) stop("Fit must be a gstatVariogram!\n")
if(!is(fit, "variogramModel")) stop("Fit must be a variogramModel!\n")
if(length(fit.method)!=1) stop("One fit.method must be selected!\n")
if(!any(fit.method==c(1,2,6,7))) stop("The selected fit.method is not managed!\n")
SSErr<-attr(fit,"SSErr")
if(fit.method==1){
# weights : $N_j$
# with $N_j$ the number of point pairs.
weig<-vario$np
SSTot<- sum( (weig* (vario$gamma-mean(vario$gamma)))^2 )
}
if(fit.method==2){
# weights : $N_j/{?(h_j)}^2$
# with $N_j$ the number of point pairs and $?(h_j)$ the variance of the group of point-pairs.
weig<-vario$np/vario$gamma^2
SSTot<- sum( (weig * (vario$gamma-mean(vario$gamma)))^2 )
}
if(fit.method==6){
# this method uses no weights
SSTot<- sum((vario$gamma-mean(vario$gamma))^2 )
}
if(fit.method==7){
# The default method used by fit.variogram
# weights : $N_j/h_j^2$
# with $N_j$ the number of point pairs and $h_j$ the distance.
weig<-vario$np/vario$dist^2
SSTot<- sum( (weig * (vario$gamma-mean(vario$gamma)))^2 )
}
SSReg <- SSTot-SSErr
R2 <- 1-SSErr/SSTot
cat("\nFind2R values -> SSErr: ", SSErr, " SSReg: ", SSReg, " SSTot: ", SSTot, " R2: ", R2, "\n")
return(R2)
}
x<-c(728.42462,730.28539,716.32417,709.38916,723.78518,711.48527,710.63429,702.78617,711.34134,727.15081,733.52192,747.83184,762.91814,772.14823,754.44948,756.01782,782.16277,805.00559,791.57293,778.95638,751.28784,747.50682,729.00907,737.52659,752.25245,765.02277,775.26606,758.40180,761.04191,786.26503,801.47612,805.93667,784.75853,811.48931,808.96591,810.76875,833.36320,852.31984,848.53256,862.31869,864.78162,874.39381,860.98220,863.13351,846.52117,839.26401,823.99134,824.17076,808.31054,796.94449,788.72534,783.46914,760.01992,748.42461,704.66752,705.03418,706.80017,735.23510,760.99192,748.97461,743.58937,757.76718,745.37049,735.69955,713.69304,711.76608,696.13015,698.68165,703.38914,697.03829,674.89737,691.64862,686.96580,670.46693,660.36596,652.17736,655.94281,631.23998,625.32977,608.40446,600.01346,585.69355,565.28714,556.27343,536.13812,515.69045,501.85434,468.78282,452.12268,462.46920,477.41710,502.99157,526.49228,533.27104,565.42016,580.46385,593.30132,597.04050,612.55727,616.51646,632.00233,639.82386,618.77567,632.42059,699.40611,698.19514,688.99576,683.27393,664.85813,672.86755,673.97141,660.74718,651.24705,618.13733,630.26338,611.28937,610.41164,599.03826,576.78329,561.36347,556.46450,553.23343,546.01103,554.50072,542.04608,521.94745,510.27728,512.65755,488.60461,405.01785,361.92731,346.21866,340.41320,380.65148,383.12679,370.70358,351.20737,300.50567,247.95438,262.63726,268.59934,270.24494,281.41292,297.47110,320.07809,313.60984,357.62665,350.65600,384.75226,401.94217,408.48164,422.53176,429.00488,425.17268,437.22270,456.20930,461.29643,474.43676,458.18572,439.57278,421.30108,422.40514,405.29465,392.81110,374.30654,351.43988,317.06783,333.12881,309.35842,265.56919,255.88517,245.52172,229.21670,233.48720,243.23252,259.24005,253.14360,239.99223,214.56455,183.59160,168.31185,152.99117,138.34226,137.52623,115.98087,107.06663,99.69656,85.75688,78.31777,70.17954,79.94090);
y<-c(124.7121,155.7824,162.0367,174.9860,182.7883,193.3682,204.0938,209.6457,226.2255,230.6839,216.3763,206.8596,211.5037,220.2471,224.9453,232.8240,237.5390,247.7069,248.0330,252.4336,247.5189,255.2388,260.4527,271.2947,280.9344,296.3364,320.4039,321.0515,335.8206,348.7338,362.3612,373.6949,380.5346,408.6241,418.1795,440.9424,435.0106,458.6167,476.4112,472.4167,481.8280,497.4812,511.2659,524.5774,505.7546,518.3096,505.7870,475.0471,474.8810,474.4909,456.7841,449.8332,427.9383,425.1974,420.6297,402.9187,397.7280,392.1504,403.7409,393.0887,374.5793,362.1227,340.7810,329.5665,332.4625,314.0868,319.9376,303.6359,284.9916,261.5924,254.7595,243.9278,223.1862,209.6194,218.8728,233.3998,248.0205,228.0175,237.4241,238.7267,227.7562,234.9984,230.4800,246.1515,245.6719,259.3798,271.8863,276.0670,281.5663,304.5483,288.5204,297.2949,298.9224,285.7951,273.1466,292.0257,282.3901,276.2681,254.0934,269.8184,271.7436,287.5611,288.3422,296.2242,339.5946,351.0233,387.9683,380.3363,371.8982,363.1685,347.2458,357.5835,366.7113,391.5945,387.2011,369.5583,380.6504,380.5077,367.6701,392.5244,381.4005,375.1102,390.4030,403.0003,417.9508,396.0499,386.9360,416.5784,416.0306,406.7233,395.8752,394.1087,375.7258,383.5885,408.1062,407.8005,402.2309,380.0658,393.0108,383.8640,384.5255,377.9253,359.1873,365.2347,371.8960,354.6252,349.7077,360.4012,349.1288,345.7168,360.3246,356.7448,346.6512,360.1093,356.0386,344.5684,325.8377,324.5357,310.5146,314.3928,306.2765,314.1672,313.2553,310.6806,321.2840,321.6360,321.7791,337.2096,350.8569,340.5812,359.0478,375.2724,371.9662,361.9425,353.8320,323.5796,305.8561,311.9475,297.7328,288.8309,298.5439,289.7562,276.5800,255.2413,242.7040,259.3922,254.0665,236.4128,232.4118,214.5125,222.0566);
z<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,0,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1);
# The variogram used in the other soft.
v.m=vgm(0.311, "Sph", 482, 0.049)
dat=data.frame(x,y,z)
coordinates(dat) = ~x+y
vario=variogram(z~1, dat, width=68.28, cutoff=800);
# The problem I face is that the other software result is
# R2: 0.982 and Residual Sum of Square: 1.65.1E-03
# So for fit.method=7
f.m=7
v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m);
r2=findR2(vario, v.fit, f.m)
# So for fit.method=6
f.m=6
v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m)
r2=findR2(vario, v.fit, f.m)
f.m=2
v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m)
r2=findR2(vario, v.fit, f.m)
f.m=1
v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m)
r2=findR2(vario, v.fit, f.m)
R2 values from SSErr fit.variogram attribute
2 messages · Sébastien Durand, Edzer Pebesma
Sebastien, let me first express my happiness that your head is clearer now; I hope that I can speak for the other 1990 subscribers of this list too. With which software do you compare? And how does that software exactly document how it computes R2 values? In your message of Mar 7 to this list you quote a reply from me to kili at grf.rs. In the computations you've showed so far you're still ignoring the comment I made there.
On 03/14/2011 06:07 PM, S?bastien Durand wrote:
Dear Dr. Pebesma,
I appreciate your help very much and I agree with you!
My head is clearer now but I still have some questions.
My first question is if I set fit.variogram to fit.sills=FALSE, fit.ranges=FALSE, to I get computation of SSErr... (I guess so).
That been said lets move forward :
Based on Geostatistic Modeling Spatial uncertainty (Chiles) page 109, in
the weighted least square method, the weight are also squared :
$sum{w^2 * [?(h_j)-?(h_j:b)]^2}$
So this agrees with what you have suggested me (and which I have applied). THANK GOD
Now, based on the gstat user's manual page 42 table 4.2 values for fit.
# 1 = $N_j$
# 2 = $N_j/{?(h_j)}^2$
# 6 = no weights (OLS) Ordinary Least Square
# 7 = $N_j/h_j^2$
I have made this code... that computes the R2 based of the fitted method used.
I have also include an example of the data... and my code.
If you tell me the code is correct, I will move forward, but I still do not comprehend why I cannot obtain similar values.
I must be missing something... any clue...???
I how the thing about comparing results with other software is a risky business but it is statistics, so I would at least be happy with simial result but look what I get... and tell me this is not mind puzzling?
Anyway I hope this code is good:
findR2 <- function(vario, fit, fit.method=c(1,2,6,7)){
# fit = fit.variogram output
# vario = variogram output
if(!is(vario, "gstatVariogram")) stop("Fit must be a gstatVariogram!\n")
if(!is(fit, "variogramModel")) stop("Fit must be a variogramModel!\n")
if(length(fit.method)!=1) stop("One fit.method must be selected!\n")
if(!any(fit.method==c(1,2,6,7))) stop("The selected fit.method is not managed!\n")
SSErr<-attr(fit,"SSErr")
if(fit.method==1){
# weights : $N_j$
# with $N_j$ the number of point pairs.
weig<-vario$np
SSTot<- sum( (weig* (vario$gamma-mean(vario$gamma)))^2 )
}
if(fit.method==2){
# weights : $N_j/{?(h_j)}^2$
# with $N_j$ the number of point pairs and $?(h_j)$ the variance of the group of point-pairs.
weig<-vario$np/vario$gamma^2
SSTot<- sum( (weig * (vario$gamma-mean(vario$gamma)))^2 )
}
if(fit.method==6){
# this method uses no weights
SSTot<- sum((vario$gamma-mean(vario$gamma))^2 )
}
if(fit.method==7){
# The default method used by fit.variogram
# weights : $N_j/h_j^2$
# with $N_j$ the number of point pairs and $h_j$ the distance.
weig<-vario$np/vario$dist^2
SSTot<- sum( (weig * (vario$gamma-mean(vario$gamma)))^2 )
}
SSReg <- SSTot-SSErr
R2 <- 1-SSErr/SSTot
cat("\nFind2R values -> SSErr: ", SSErr, " SSReg: ", SSReg, " SSTot: ", SSTot, " R2: ", R2, "\n")
return(R2)
}
x<-c(728.42462,730.28539,716.32417,709.38916,723.78518,711.48527,710.63429,702.78617,711.34134,727.15081,733.52192,747.83184,762.91814,772.14823,754.44948,756.01782,782.16277,805.00559,791.57293,778.95638,751.28784,747.50682,729.00907,737.52659,752.25245,765.02277,775.26606,758.40180,761.04191,786.26503,801.47612,805.93667,784.75853,811.48931,808.96591,810.76875,833.36320,852.31984,848.53256,862.31869,864.78162,874.39381,860.98220,863.13351,846.52117,839.26401,823.99134,824.17076,808.31054,796.94449,788.72534,783.46914,760.01992,748.42461,704.66752,705.03418,706.80017,735.23510,760.99192,748.97461,743.58937,757.76718,745.37049,735.69955,713.69304,711.76608,696.13015,698.68165,703.38914,697.03829,674.89737,691.64862,686.96580,670.46693,660.36596,652.17736,655.94281,631.23998,625.32977,608.40446,600.01346,585.69355,565.28714,556.27343,536.13812,515.69045,501.85434,468.78282,452.12268,462.46920,477.41710,502.99157,526.49228,533.27104,565.42016,580.46385,593.30132,597.04050,612
.55727,616.51646,632.00233,639.82386,618.77567,632.42059,699.40611,698.19514,688.99576,683.27393,664.85813,672.86755,673.97141,660.74718,651.24705,618.13733,630.26338,611.28937,610.41164,599.03826,576.78329,561.36347,556.46450,553.23343,546.01103,554.50072,542.04608,521.94745,510.27728,512.65755,488.60461,405.01785,361.92731,346.21866,340.41320,380.65148,383.12679,370.70358,351.20737,300.50567,247.95438,262.63726,268.59934,270.24494,281.41292,297.47110,320.07809,313.60984,357.62665,350.65600,384.75226,401.94217,408.48164,422.53176,429.00488,425.17268,437.22270,456.20930,461.29643,474.43676,458.18572,439.57278,421.30108,422.40514,405.29465,392.81110,374.30654,351.43988,317.06783,333.12881,309.35842,265.56919,255.88517,245.52172,229.21670,233.48720,243.23252,259.24005,253.14360,239.99223,214.56455,183.59160,168.31185,152.99117,138.34226,137.52623,115.98087,107.06663,99.69656,85.75688,78.31777,70.17954,79.94090);
y<-c(124.7121,155.7824,162.0367,174.9860,182.7883,193.3682,204.0938,209.6457,226.2255,230.6839,216.3763,206.8596,211.5037,220.2471,224.9453,232.8240,237.5390,247.7069,248.0330,252.4336,247.5189,255.2388,260.4527,271.2947,280.9344,296.3364,320.4039,321.0515,335.8206,348.7338,362.3612,373.6949,380.5346,408.6241,418.1795,440.9424,435.0106,458.6167,476.4112,472.4167,481.8280,497.4812,511.2659,524.5774,505.7546,518.3096,505.7870,475.0471,474.8810,474.4909,456.7841,449.8332,427.9383,425.1974,420.6297,402.9187,397.7280,392.1504,403.7409,393.0887,374.5793,362.1227,340.7810,329.5665,332.4625,314.0868,319.9376,303.6359,284.9916,261.5924,254.7595,243.9278,223.1862,209.6194,218.8728,233.3998,248.0205,228.0175,237.4241,238.7267,227.7562,234.9984,230.4800,246.1515,245.6719,259.3798,271.8863,276.0670,281.5663,304.5483,288.5204,297.2949,298.9224,285.7951,273.1466,292.0257,282.3901,276.2681,254.0934,269.8184,271.7436,287.5611,288.3422,296.2242,339.5946,351.0233,387.9683,380.3363,371.8982,36
3.1685,347.2458,357.5835,366.7113,391.5945,387.2011,369.5583,380.6504,380.5077,367.6701,392.5244,381.4005,375.1102,390.4030,403.0003,417.9508,396.0499,386.9360,416.5784,416.0306,406.7233,395.8752,394.1087,375.7258,383.5885,408.1062,407.8005,402.2309,380.0658,393.0108,383.8640,384.5255,377.9253,359.1873,365.2347,371.8960,354.6252,349.7077,360.4012,349.1288,345.7168,360.3246,356.7448,346.6512,360.1093,356.0386,344.5684,325.8377,324.5357,310.5146,314.3928,306.2765,314.1672,313.2553,310.6806,321.2840,321.6360,321.7791,337.2096,350.8569,340.5812,359.0478,375.2724,371.9662,361.9425,353.8320,323.5796,305.8561,311.9475,297.7328,288.8309,298.5439,289.7562,276.5800,255.2413,242.7040,259.3922,254.0665,236.4128,232.4118,214.5125,222.0566);
z<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,0,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1); # The variogram used in the other soft. v.m=vgm(0.311, "Sph", 482, 0.049) dat=data.frame(x,y,z) coordinates(dat) = ~x+y vario=variogram(z~1, dat, width=68.28, cutoff=800); # The problem I face is that the other software result is # R2: 0.982 and Residual Sum of Square: 1.65.1E-03 # So for fit.method=7 f.m=7 v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m); r2=findR2(vario, v.fit, f.m) # So for fit.method=6 f.m=6 v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m) r2=findR2(vario, v.fit, f.m) f.m=2 v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m) r2=findR2(vario, v.fit, f.m) f.m=1 v.fit=fit.variogram(vario, v.m, fit.sills=FALSE, fit.ranges=FALSE, fit.method=f.m) r2=findR2(vario, v.fit, f.m)
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Edzer Pebesma Institute for Geoinformatics (ifgi), University of M?nster Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics e.pebesma at wwu.de