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

Sébastien Durand v8extra at gmail.com
Mon Mar 14 18:07:17 CET 2011 

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)

More information about the R-sig-Geo mailing list