Closed PatWright closed 5 years ago
Sounds like a problem with spaMM. So before I dive into it, did you make sure that you have the latest version of spaMM installed and that your model have been fitted using the same version of spaMM than the one you are using to run calibfit? In case of doubt I stored the spaMM version in the object created by isofit and to know which spaMM version you are actually using, just type library(spaMM) and it should be displayed.
Yes, I think it is an issue with spaMM. I'm currently using spaMM 2.5.0.
and did you fit the models with the same version?
Yes.
[models fitted with spaMM version 2.5.0]
I see... The current version of spaMM has a bug for predictions with offsets, which we use here, so it may come from that. Here is the development version of spaMM (the patched version will be on CRAN soon): spaMM_2.5.4.tar.gz. Please install it, refit the model and try to fit the calibration. Let me know if that solves the problem.
Hi,
as so often, without debugging information, it is a waste of time to speculate... could you show us the traceback that appears when the error occurs (affter setting options(error=recover) or using the Rstudio traceback feature) ? .get_cP_stuff() is internal to spaMM and I currently do not understand in which context it is called here. I do not see how the newer spaMM version would change anything.
F.
@PatWright, @f-rousset is the developer of spaMM. Follow his advice (and forget mine). I will assist you if you don't understand. ++
Just to be sure, @PatWright are you using IsoriX 0.8.1?
Is this enough information? Indeed, the latest version didn't help. I started again from scratch with 2.5.0 and it's working again.
> CalibFA<- calibfit(data = TotFurAut, isofit = FitAut)#furaut
[1] "predicting the isoscape value in each calibration site..."
Error in .get_cP_stuff(ranFix, "nu", which = char_rd) :
could not find function ".get_cP_stuff"
Timing stopped at: 0.03 0 0.03
> traceback()
10: MaternCorr.default(nu = .get_cP_stuff(ranFix, "nu", which = char_rd),
Nugget = .get_cP_stuff(ranFix, "Nugget", which = char_rd),
d = distmat)
9: MaternCorr(nu = .get_cP_stuff(ranFix, "nu", which = char_rd),
Nugget = .get_cP_stuff(ranFix, "Nugget", which = char_rd),
d = distmat)
8: .get_sub_corr_info(object)$corr_families[[old_rd]]$calc_corr_from_dist(ranFix = object$ranFix,
char_rd = old_char_rd, distmat = uuCnewold)
7: structure(.get_sub_corr_info(object)$corr_families[[old_rd]]$calc_corr_from_dist(ranFix = object$ranFix,
char_rd = old_char_rd, distmat = uuCnewold), ranefs = ranefs[[new_rd]])
6: .make_new_corr_lists(object = object, locdata = locdata, which_mats = which_mats,
newZAlist = newZAlist, newinold = newinold)
5: .calc_new_X_ZAC(object = object, newdata = newdata, re.form = re.form,
variances = variances)
4: .predict_body(object = object, newdata = newdata, re.form = re.form,
variances = variances, binding = binding, intervals = intervals,
level = level, blockSize = blockSize)
3: spaMM::predict.HLfit(isofit[["mean_fit"]], newdata = data, variances = list(predVar = TRUE,
cov = TRUE))
2: system.time({
if (any(class(isofit) %in% "MULTIISOFIT")) {
stop("object 'isofit' of class MULTIISOFIT; calibration have not yet been implemented for this situation.")
}
species_rand <- FALSE
species_info <- any(colnames(data) %in% "species_ID")
if (!is.null(species_rand)) {
if (!species_info & species_rand) {
stop("The random effect for species cannot be fit if data does not contain a column called species_ID")
}
}
data <- .prepare_data_calib(data)
if (!species_info) {
species_rand <- FALSE
}
else {
nb_species <- length(unique(data$species_ID))
if (is.null(species_rand) & nb_species > 4) {
species_rand <- TRUE
}
else {
species_rand <- FALSE
}
}
if (verbose) {
print("predicting the isoscape value in each calibration site...")
}
mean_calib <- spaMM::predict.HLfit(isofit[["mean_fit"]],
newdata = data, variances = list(predVar = TRUE, cov = TRUE))
data$mean_source_value <- c(mean_calib)
predcov_matrix_isofit_full <- attr(mean_calib, "predVar")
data$mean_predVar_source <- diag(predcov_matrix_isofit_full)
firstoccurences <- match(levels(data$site_ID), data$site_ID)
predcov_isofit <- predcov_matrix_isofit_full[firstoccurences,
firstoccurences]
rownames(predcov_isofit) <- levels(data$site_ID)
colnames(predcov_isofit) <- levels(data$site_ID)
if (verbose) {
print("fitting the calibration function...")
}
objective_fn_calib <- function(param, data, predcov, species_rand,
return_fit = FALSE, lik_method = "REML") {
data$intercept <- param[1]
data$slope <- param[2]
lambda_list <- list(lambda = c(1e-06 + unique(data$slope)^2,
NA))
calib_formula <- "sample_value ~ 0 + offset(intercept+slope*mean_source_value) +\n corrMatrix(1|site_ID) + (1|site_ID)"
if (species_rand) {
lambda_list$lambda <- c(lambda_list$lambda, NA)
calib_formula <- paste(calib_formula, "+ (1|species_ID)")
}
calib_fit <- spaMM::HLCor(formula = stats::formula(calib_formula),
corrMatrix = predcov, ranPars = lambda_list, data = data,
method = lik_method)
if (return_fit)
return(calib_fit)
return(calib_fit$APHLs$p_v)
}
opt_res <- stats::optim(par = c(0, 1), fn = objective_fn_calib,
control = c(list(fnscale = -1), control_optim), data = data,
predcov = predcov_isofit, species_rand = species_rand,
lik_method = "REML")
param_calibfit <- opt_res$par
names(param_calibfit) <- c("intercept", "slope")
calib_fit <- objective_fn_calib(param = param_calibfit, data = data,
predcov = predcov_isofit, species_rand = species_rand,
lik_method = "REML", return_fit = TRUE)
if (verbose) {
print("computing the covariance matrix of fixed effects...")
}
fixefCov_calibfit <- solve(-numDeriv::hessian(objective_fn_calib,
param_calibfit, data = data, predcov = predcov_isofit,
species_rand = species_rand, lik_method = "ML"))
rownames(fixefCov_calibfit) <- names(param_calibfit)
colnames(fixefCov_calibfit) <- names(param_calibfit)
})
1: calibfit(data = TotFurAut, isofit = FitAut)
@courtiol Yes, IsoriX 0.8.1
still puzzling. So the next step is
options(error = quote({dump.frames(to.file = TRUE)}))
and send me the resulting dump file that is created upon error.
Please also paste here the result of
head(TotFurAut)
and the call you used to create FitAut
.
One more quick question (but I will be able to seethe answer in the dump). With which spaMM version was created object FitAut
? I.e., what is the value of FitAut$[["mean_fit"]]$spaMM.version
?
... and what is the sessionInfo()
of the session where the error occurs...
I hope this info can help. It seems to work fine now if I start from the beginnig.
@courtiol
FitAut <- isofit(data = GNIPDataAutumn,
mean_model_fix = list(elev = TRUE, lat_abs = TRUE))
head(TotFurAut)
site_ID long lat elev sample_ID species sample_value
1 Alvaiazere -8.41042 39.83198 595.2 ALV101 Miniopterus -23.21
2 Alvaiazere -8.41042 39.83198 595.2 ALV102 Miniopterus -15.86
@f-rousset
FitAut$mean_fit$spaMM.version
[1] ‘2.5.0’
sessionInfo()
R version 3.4.3 (2017-11-30)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
Matrix products: default
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spaMM_2.5.0 raster_2.6-7 sp_1.3-1 IsoriX_0.8.1 gridExtra_2.3 digest_0.6.15
[7] ggplot2_3.0.0
loaded via a namespace (and not attached):
[1] Rcpp_0.12.18 pillar_1.3.0 compiler_3.4.3 RColorBrewer_1.1-2
[5] plyr_1.8.4 bindr_0.1.1 tools_3.4.3 rasterVis_0.45
[9] nlme_3.1-131 tibble_1.4.2 gtable_0.2.0 lattice_0.20-35
[13] viridisLite_0.3.0 rlang_0.2.2 Matrix_1.2-12 rstudioapi_0.7
[17] parallel_3.4.3 hexbin_1.27.2 withr_2.1.2 grid_3.4.3
[21] R6_2.2.2 pbapply_1.3-4 latticeExtra_0.6-28 scales_1.0.0
[25] MASS_7.3-47 assertthat_0.2.0 colorspace_1.3-2 labeling_0.3
[29] proxy_0.4-22 lazyeval_0.2.1 munsell_0.5.0 crayon_1.3.4
[33] zoo_1.8-3
$`calibfit(data = TotFurAut, isofit = FitAut)`
<environment: 0x000000003330d158>
$`system.time({\n if (any(class(isofit) %in% "MULTIISOFIT")) {\n stop("object 'isofi`
<environment: 0x000000003330d2e0>
$`spaMM::predict.HLfit(isofit[["mean_fit"]], newdata = data, variances = list(predVar = TRUE`
<environment: 0x000000002bbce390>
$`.predict_body(object = object, newdata = newdata, re.form = re.form, variances = variances`
<environment: 0x0000000031a948b8>
$`.calc_new_X_ZAC(object = object, newdata = newdata, re.form = re.form, variances = varianc`
<environment: 0x0000000031a94b90>
$`.make_new_corr_lists(object = object, locdata = locdata, which_mats = which_mats, newZAlis`
<environment: 0x0000000033712b20>
$`structure(.get_sub_corr_info(object)$corr_families[[old_rd]]$calc_corr_from_dist(ranFix = `
<environment: 0x000000003366e3b0>
$`.get_sub_corr_info(object)$corr_families[[old_rd]]$calc_corr_from_dist(ranFix = object$ran`
<environment: 0x000000003366c828>
$`MaternCorr(nu = .get_cP_stuff(ranFix, "nu", which = char_rd), Nugget = .get_cP_stuff(ranFi`
<environment: 0x000000003366cac8>
$`MaternCorr.default(nu = .get_cP_stuff(ranFix, "nu", which = char_rd), Nugget = .get_cP_stu`
<environment: 0x000000003366c150>
attr(,"error.message")
[1] "Error in .get_cP_stuff(ranFix, \"nu\", which = char_rd) : \n could not find function \".get_cP_stuff\"\n"
attr(,"class")
[1] "dump.frames"
@PatWright
Thanks for this. Is the dump file somewhere that I missed ?
@PatWright we are not sure of the current state of your situation. Is everything actually now working?
Sorry, I've been away. Yes, it works. I just need to redo my isoscapes when I reopen my environment.
This is presumably a tricky issue with respect to the namespaces attached to a saved fit object and/or to member functions defined within the fit object. Hence I guess it depends on the environment of the sesion where a fit object was created and how it was saved to file. I am trying to figure a way to avoid problems but I am unable to recreate them.
As an illustration of what I am concerned about, in a virgin session with only IsoriX attached, I run the isofit
example to create object GermanFit
, I saved GermanFit
to a file, then started another virgin session where I loaded this object. Simply loading this object has the effect that all packages loaded as a namespace in the first session are loaded in the second session (this includes spaMM since IsoriX imports spaMM). Without attaching any package, I then can use the predict
function from spaMM to predict from GermanFit$mean_fit
in new positions, which calls .get_cP_stuff()
. No bug occurs. I fail to see why it would differ on your side.
Hi there,
Just started to have this issue with calibfit. I was wondering if anyone knew how to fix it? It's a code that worked fine with the same dataset last week.
Thanks, Patrick