[R] GSTAT Co-Kriging: How can i fix a non positive definte correlation matrix

Mon Jan 5 23:54:40 CET 2015

QUESTION: Can anyone give me advice on how to resolve an error message regarding a non positive definite correlation matrix so that the computed negative cross-variogram can be used to perform co-kriging?

I am performing co-kriging on negatively correlated parameters. 

data = dtrd 

x	y	n.ts	n.ae
3.2	6.4	-0.052514693	0.055104755
6.4	6.4	-0.375234492	-0.861381682
9.6	6.4	1.362810557	-0.23577909
12.8	6.4	-0.333087778	0.134741672
16	6.4	0.699916885	-2.348188189
22.4	6.4	0.605552409	-0.75492334
25.6	6.4	0.406868496	0.329847607
3.2	9.6	-0.529754386	0.251549504
6.4	9.6	-0.416509276	1.094807283
9.6	9.6	-0.450700276	0.988002955
12.8	9.6	-0.732937251	1.344413512
16	9.6	0.377457633	-0.542187169
19.2	9.6	-0.989029229	2.002433474
22.4	9.6	-1.259551904	-0.895859924
25.6	9.6	-0.470220669	0.370809323
3.2	12.8	1.046426438	0.168121785
6.4	12.8	-1.440639861	1.423433672
9.6	12.8	-0.899563555	0.83844681
12.8	12.8	-1.287567164	-0.086478673
16	12.8	-0.816931437	1.156633451
19.2	12.8	-0.782655308	0.299641538
22.4	12.8	0.048342822	-0.763091614
3.2	16	-0.5581308	-0.130509134
6.4	16	-0.217874377	0.815852429
9.6	16	-0.272702769	0.906621066
12.8	16	-0.996064	0.274459785
16	16	-1.97520507	2.017626766
19.2	16	-0.718163458	0.444173367
22.4	16	0.337266312	0.831151021
25.6	16	-0.565668608	0.513114809
3.2	19.2	0.651736313	-0.055114304
6.4	19.2	-0.715751462	0.205488257
9.6	19.2	0.650820529	0.558032798
12.8	19.2	-1.184380204	0.612006591
16	19.2	0.132388601	0.462059331
19.2	19.2	1.08502868	0.399240335
22.4	19.2	-0.85196341	0.113079949
25.6	19.2	0.099887288	-1.922971577
3.2	22.4	-0.233091874	-0.287556625
6.4	22.4	-0.989196908	0.10047403
9.6	22.4	-0.025054063	0.384754057
12.8	22.4	-0.394999955	-1.272051846
16	22.4	0.049839033	-2.568551091
19.2	22.4	1.327518037	-0.974203156
22.4	22.4	-0.421157849	-0.0101211
25.6	22.4	-0.178717505	0.408401438
3.2	25.6	-0.065413063	-0.195765074
6.4	25.6	1.167410569	-1.626236059
9.6	25.6	0.428968564	-1.287922066
12.8	25.6	1.550074254	-1.011990828
16	25.6	2.931334909	-1.872352342
19.2	25.6	4.126464345	-2.254471152
22.4	25.6	0.44561262	-0.709554064
25.6	25.6	0.543591219	1.03999109
3.2	28.8	0.716228164	-0.671315106
6.4	28.8	1.856183532	-1.403466298
9.6	28.8	0.774644882	-0.952601604
12.8	28.8	0.425336383	-0.776073605
16	28.8	0.294908064	-0.00019281
19.2	28.8	2.38260471	-0.778811406
22.4	28.8	-0.315391214	-0.223383771
25.6	28.8	1.136916243	-0.238938091
3.2	32	-0.233091874	-0.520269718
6.4	32	-0.352017426	-0.556041646
9.6	32	-0.933099317	-0.441354858
16	32	-0.865945243	-1.011088256
19.2	32	-0.594339105	-0.399550741
22.4	32	-0.653328511	1.213119402
25.6	32	0.076670222	0.819568017
3.2	35.2	-0.14022361	0.166527242
6.4	35.2	-0.545492977	-0.613535469
9.6	35.2	-0.084386566	-0.423754708
12.8	35.2	2.107283841	-0.891753223
16	35.2	-1.645006796	0.566968259
22.4	35.2	0.512684145	0.515250895
25.6	35.2	0.517794479	-1.144578624
3.2	38.4	-0.359495901	0.191964723
6.4	38.4	0.548288806	-1.011825356
9.6	38.4	-0.597741695	0.046560408
12.8	38.4	-1.003803022	0.319904274
16	38.4	0.060157729	0.824035747
19.2	38.4	0.251793973	0.632946256
22.4	38.4	-1.148625922	1.257360463
25.6	38.4	-0.957779058	0.30201831
3.2	41.6	0.267364885	1.139304073
6.4	41.6	0.171656399	0.623288738
9.6	41.6	-0.912461925	0.820049388
12.8	41.6	-0.645228334	1.041375034
16	41.6	0.83148026	2.914000894
19.2	41.6	-1.187664129	2.756035789
22.4	41.6	-0.467591981	0.467880919
25.6	41.6	-0.121964677	0.878460826
3.2	44.8	-0.158281328	-0.266451488
6.4	44.8	-0.235932095	1.440612621
12.8	44.8	-0.436274739	-0.298237058
16	44.8	-1.136811015	-0.688794913
19.2	44.8	0.87349541	0.449964869
22.4	44.8	-0.661067533	-0.08882536
25.6	44.8	0.301101862	-0.255259598
3.2	48	1.265698729	0.948786211
6.4	48	-0.806040053	0.352201302
9.6	48	1.006815543	1.12999254
12.8	48	-0.0132082	0.581935908
16	48	-1.214201235	1.431842632
19.2	48	-0.323473333	-0.083379843
22.4	48	-0.003250658	-0.248324838
25.6	48	-0.168398809	-0.916859761
3.2	51.2	-0.854793312	-0.231928117
6.4	51.2	0.507014021	-0.49107152
9.6	51.2	3.104090518	-2.654897126
12.8	51.2	-0.193785381	0.067244345
16	51.2	-0.71632415	-1.46047875
19.2	51.2	2.480632323	-0.336852087
22.4	51.2	-0.534663506	-1.093553232
25.6	51.2	-0.124544351	-0.603005465

Using the following code I compute a negative cross-variogram:

dtrd = as.data.frame(dtrd)
# create a grid onto which i will interpolate
# first get a range in data
x.range <- as.integer(range(dtrd[,1]))
y.range <- as.integer(range(dtrd[,2]))
# now expand to a grid with desired spacing
grid <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.4), 
y=seq(from=y.range[1], to=y.range[2], by=1.6))
plot(grid)

# Convert to spatial data frame
coordinates(grid) = c("x", "y")
coordinates(dtrd) = c("x", "y")

# Compute cross-variogram 
g <- gstat(id = "n.ae", formula = n.ae ~ + 1, data = dtrd, nmax=10)
g <- gstat(g, "n.ts", n.ts ~ 1, dtrd, nmax=10)
g1 <- gstat(g, id=c("n.ae", "n.ts"), model = vgm(psill=-0.65, model="Sph", range=15, nugget=-0.05), fill.all = T)
v <- variogram(g)
g.fit<-fit.lmc(v, g1, correct.diagonal = 1.01)

# Plot cross-variogram with fitted model
plot(v[v$id == "n.ae.n.ts",], g.fit$model[[2]], main="Cross-variogram: 
Air Entry and Saturated Water Content", xlab="Lag Distance (mm)", 
ylab="Semivariance (-)", ylim=c(0,-0.8))

When attempting to use the negative cross-variogram in the co-kriging plan I receive the following error message:

# Perform co-kriging predictions
ck <- predict.gstat(g1, grid)

non-positive definite coefficient matrix in structure 1non-positive definite coefficient matrix in structure 2Warning: No Intrinsic Correlation or Linear Model of Coregionalization found
Reason: coefficient matrix not positive definite
Warning: [add `set nocheck = 1;' to the command file to ignore the following error]

Error in predict.gstat(g1, grid) : 
  gstat: value not allowed for: variograms do not satisfy a legal model

I have successfully produced kriged maps of each of the parameters individually, but do not know exactly how to deal with the negative cross-variogram and resolve the issue of a non positive definite correlation matrix. I think the problem arises from using a negative sill and nugget to model the sample cross-variogram with the fit.lmc function. If I use the same code but make the sill and nugget positive I get the same fit of the model to the sample cross-variogram and a message reading : Intrinsic Correlation found. Good. [using ordinary cokriging], which I initial thought that this was a good sign, but once again I can not resolve the following error message:

# Compute cross-variogram 
g <- gstat(id = "n.ae", formula = n.ae ~ + 1, data = dtrd, nmax=10)
g <- gstat(g, "n.ts", n.ts ~ 1, dtrd, nmax=10)
g1 <- gstat(g, id=c("n.ae", "n.ts"), model = vgm(psill=0.65, model="Sph", range=15, nugget=0.05), fill.all = T)
v <- variogram(g)
g.fit<-fit.lmc(v, g1, correct.diagonal = 1.01)

# Perform co-kriging predictions
ck <- predict.gstat(g1, grid)

Intrinsic Correlation found. Good. 
[using ordinary cokriging]

ERROR: "chfactor.c", line 131: singular matrix in function LDLfactor()
Error in predict.gstat(g1, grid) : LDLfactor

The gstat manual shows that this problem can be fixed by adding correct.diagonal = 1.01 to the fit.lmc function, but this does not eliminate the error message

I am very confused about all of this and would appreciate any help to fix the non positive correlation matrix so that the negative cross-variogram can be used to perform co-kriging.

Thanks in advance for any help,

Blake Roberts