[R-sig-Geo] R2 values from SSErr fit.variogram attribute

Mon Mar 14 20:20:44 CET 2011

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