R crushes when I run anim.fit function in windows 10

[Mangeni-Sande]

R crushes when I run anim.fit () function below.
Below is part of the script I am using map dagaa densities on Lake Victoria. The data is also attached. Please help me to identify where the problem culd be. Thanks

Code: Select all

## Load required packages
library(Rcpp)
library(RGeostats) # RGeostats_11.2.1
library(spatstat)  # spatstat_1.55-1
library(colorspace) # colorspace_1.4-0
library(dplyr)
library(fields)
library(ggplot2)
library(mapdata)
areaPoly <- read.csv("lakeVictoriaRGeostatPolygon.csv", header=T)
dat <- read.csv("dagaa_NASC_by_EDSU3.csv", header=T)# Fish density data
densDb <- db.create(dat, flag.grid=F)#flag.grid=F indicates do not overlay the coastline
densDb = db.locate(densDb, "long", "x",1)#attribute long corresponds to the first (x1) coordinate
densDb = db.locate(densDb, "lat", "x",2)#
densDb = db.locate(densDb, "NASC","z")#attribute of interest NASC is desiginated locator "Z"
poly.data <- polygon.create(areaPoly)
projec.query()# the projection is turned off
plot(densDb, zmin=0.001, pch.low=3, cex.low=0.25, las=1,
     pch=21, col=1, inches=5, title="dagaa_2015", asp=1)#asp=1, distances btn pts
densDb <- infl(densDb, nodes=c(400,400), origin=c(min(areaPoly$X),min(areaPoly$Y)),
              extend=c(abs(max(areaPoly$X)-min(areaPoly$X)), abs(max(areaPoly$Y))+abs(min(areaPoly$Y))),
               polygon=poly.data, plot=T, asp=1)
projec.define(projection="mean", db=densDb ,flag.update=T)
plot(densDb, zmin=0.001, pch.low=3, cex.low=0.25, las=1,
     pch=21, col=1, inches=5, title="dagaa_2015 (Mean Projection)", asp=1)
pvals <- projec.operate(densDb@items$lon,
                        densDb@items$lat)
llag <- mean(nndist(cbind(pvals$x,pvals$y)))
nlag <- round((max(pairdist(cbind(pvals$x,pvals$y)))/2)/llag)
max(pairdist(cbind(pvals$x,pvals$y)))
hist(densDb@items$NASC, main="",  xlab="",  breaks=50,  col="#7f7f7f", border="#7f7f7f")
mtext(expression(atop(s[a]~m^2~nmi^-2)), side=1, line=4.5)
mean(densDb@items$NASC)
var(densDb@items$NASC)
# NOTE: skewed distribution with many zeros and low high values
model.anam <- anam.fit(densDb,
                       type="emp", # Empirical anamorphosis with logarithmic dispersion
                       draw=T, title="Model of anamorphosis", ndisc=densDb$nech, sigma2e=40)

[Mangeni-Sande]

HI everone,

I resubmit that the anim.fit function of the code for performing conditional simulation of the backscatter data with many zero values seem not to work when using r in windows 10 . I adopted the code from RGeostats.herring.sa.scot.R demo script and Petitgas et al. 2017-Handbook of Geostatistics in R for Fisheries and Marine Ecology-page 171)

I need urgent help, please.
Thanks

[Fabien Ors]

Hi Richard

We are currently working on this issue.
To prevent anam.fit from crashing, I think that you can deactivate the draw option (draw=F).

But whatever the above workaround, I assume that anamorphosis calculations are currently broken in RGeostats.
We will come back to you as soon as the problem will be fixed (next week I hope).

Sorry for inconvenience.

[Fabien Ors]

Dear Richard,
Could you provide us the two data files lakeVictoriaRGeostatPolygon.csv and dagaa_NASC_by_EDSU3.csv please ?
Best regards,
Fabien

[Mangeni-Sande]

Dear Fabien,

Find attached the Dagaa_NASC File. Sending both files was aproblem..I will send it seperately. Thanks

[Mangeni-Sande]

Dear Richard,
Could you provide us the two data files lakeVictoriaRGeostatPolygon.csv and dagaa_NASC_by_EDSU3.csv please ?
Best regards,
Fabien

lakeVictoriaRGeostatPolygon data is here attached (changed from .csv to *.xlsb to reduce the size)

Thanks

[Mangeni-Sande]

Hi Fabien,
Surprisingly, the anamorphosis calculations work perfectly in RGeostats version 12.0.1 . The calclations are broken in other versions of RGoestatsv

Thanks

[Fabien Ors]

Nice to hear that !
Best regards

[Mangeni-Sande]

Nice to hear that !
Best regards

Hi Fabien,
Like I noted, the anim.fit works well in RGoestats 12.0.1 BUT the Gibbs sampling is not giving good results as the histogram of the truncated Gaussian NASC still contain a substantial proportion of zero values (attachement).
Could you please find out what the problem could be (I have also attached the code I am using)

Code: Select all

##############################################################################################
areaPoly <- read.csv("lakeVictoriaRGeostatPolygon.csv", header=T)
areaPoly
dat <- read.csv("dagaa_NASC_by_EDSU3.csv", header=T)# Fish density data

densDb <- db.create(dat, flag.grid=F)#flag.grid=F indicates do not overlay the coastline

densDb = db.locate(densDb, "long", "x",1)#attribute long corresponds to the first (x1) coordinate
densDb = db.locate(densDb, "lat", "x",2)#attribute lat corresponds to the second (x2) coordinate
densDb = db.locate(densDb, "NASC","z")#attribute of interest NASC is desiginated locator "Z"

poly.data <- polygon.create(areaPoly)

projec.query()# the projection is turned off

plot(densDb, zmin=0.001, pch.low=3, cex.low=0.25, las=1,
     pch=21, col=1, inches=5, title="dagaa_2015", asp=1)

densDb <- infl(densDb,
               nodes=c(400,400),
               origin=c(min(areaPoly$X),min(areaPoly$Y)),
               extend=c(abs(max(areaPoly$X)-min(areaPoly$X)),
                        abs(max(areaPoly$Y))+abs(min(areaPoly$Y))),
               polygon=poly.data,
               plot=T,
               asp=1)
projec.define(projection="mean", db=densDb ,flag.update=T)
plot(densDb, zmin=0.001, pch.low=3, cex.low=0.25, las=1,
     pch=21, col=1, inches=5,
     title="dagaa_2015 (Mean Projection)", asp=1)

pvals <- projec.operate(densDb@items$lon,
                        densDb@items$lat)
llag <- mean(nndist(cbind(pvals$x,pvals$y)))
nlag <- round((max(pairdist(cbind(pvals$x,pvals$y)))/2)/llag)

max(pairdist(cbind(pvals$x,pvals$y)))
hist(densDb@items$NASC,
     main="",
     xlab="", # Leave x-label blank
     breaks=50,
     col="#7f7f7f",
     border="#7f7f7f")
mtext(expression(atop(s[a]~m^2~nmi^-2)),
      side=1,
      line=4.5)

mean(densDb@items$NASC)
var(densDb@items$NASC)

model.anam <- anam.fit(densDb,
                       type="emp", # Empirical anamorphosis with logarithmic dispersion
                       draw=T,
                       title="Model of anamorphosis",
                       ndisc=densDb$nech,
                       sigma2e=40)

db.data.trans <- anam.z2y(densDb, anam=model.anam)

db.data.trans <- db.rename(db.data.trans,
                           name="Gaussian.NASC",
                           newname="Yp")

ycut <- qnorm(sum(db.extract(db.data.trans,"NASC") == 0) /
                length(db.extract(db.data.trans, name="NASC",
                                  flag.compress=T)))
Y <- db.extract(db.data.trans, "Yp", flag.compress=T)

hist(Y, breaks = 50, main = "",
     col="#b3b3b3", border="#b3b3b3",
     xlab="Anamorphosis Transformed NASC")

nasc <- db.extract(db.data.trans, "NASC")
Y[nasc == 0] <- ycut
db.data.trans <- db.replace(db.data.trans, "Yp", Y)

vario.Yp <- vario.calc(db.data.trans, lag=llag, nlag=nlag)
n.H <- 50
vario.Y  <- vario.trans.cut(vario.Yp, ycut, n.H)#
plot(vario.Y, xlab="Distance (n.mi)",col="red",lwd=2)
model.vario.Y <- model.auto(vario.Y, struc=melem.name(c(1,2)),col="red", draw=T)

plot(vario.Y ,npairdw=T, npairpt=F, inches=0.06, col="red", xlab="Distance (n.mi.)")
plot(vario.Yp, npairdw=T, npairpt=F, inches=0.06, col="#000000", add=T)
plot(model.vario.Y, add=T, col="blue", lwd=3)

legend("bottomright", cex=0.75, y.intersp=1.2,
       legend=c("Gaussian Transformed Variogram",
                "Truncated Variable Variogram",
                "Model Variogram"),
       col=c("red","black","blue"), lty=1, pch=c(19,19,NA), lwd=c(1,1,2))# pch controls

## a Gibbs Sampler
Ymax <- db.extract(db.data.trans, name="Yp", flag.compress=F)
Ymin <- db.extract(db.data.trans, name="Yp", flag.compress=F)
Ymin[Ymin <= ycut] <- -10
db.data.trans <- db.add(db.data.trans,Ymax)
db.data.trans <- db.locate(db.data.trans, db.data.trans$natt,"upper")
db.data.trans <- db.add(db.data.trans, Ymin)
db.data.trans <- db.locate(db.data.trans,db.data.trans$natt,"lower")
db.data.trans <- db.locate(db.data.trans, "Yp", "z")

db.data.trans <- gibbs(db = db.data.trans,
                       model = model.vario.Y,
                       seed = 1337, # Seed for reproducibility
                       nboot = 1000, # Number of iterations performed for bootstrap
                       niter = 1000, # Maximum iterations for calculating the
                       # conditional expectation
                       flag.norm = FALSE, # For normalizing the model
                       percent = 0, # Nugget effect added to model (0 = model unchanged)
                       toleps = 1, # Relative difference between the previous and the
                       # new value at a sample location below which the sample
                       # is considered as immobile
                       radix = "Gibbs", # For naming results variable
                       modify.target = TRUE) # Define locator in new variable?

db.data.trans <- db.rename(db.data.trans, "Gibbs.G1", "Y2")

# Truncated variable
histYp <- hist(db.data.trans[,"Yp"],plot=F,breaks=seq(-4,4,.1))
hist(db.data.trans[,"Yp"], proba=T, breaks=seq(-4,4,.1), 
     col="#b3b3b3", border="#b3b3b3",
     main="", xlab="Lower-cut Gaussian Variable",
     xlim=c(-4,4), ylim=c(0,ceiling(max(histYp$density))))

# Truncated variable after Gibbs sampling
hist(db.extract(db.data.trans,"Y2",flag.compress=T), proba=T, breaks=seq(-4,4,.1),
     xlab="Gaussian Transformed Variable", main="",
     col="#b3b3b3", border="#b3b3b3",
     xlim=c(-4,4), ylim=c(0,ceiling(max(histYp$density))))

[Fabien Ors]

I agree with you. We found some issues in the herring.sa.scot whole workflow.
We are currently working hard in order to fix it.
Please, wait the next RGeostats release.