Roberts, Blake
2015-Jan-05 22:54 UTC
[R] GSTAT Co-Kriging: How can i fix a non positive definte correlation matrix
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