Error in splm random effects ML model
On Wed, 15 Feb 2023, Denys Dukhovnov wrote:
Dear Roger,?
Yes, I am running R version 4.2.2 and Windows 10 OS x64, packages used
updated to the latest version on CRAN. The panel data I am using is
complete for all the variables and spatial weights, which made me wonder
why the "within" specification runs perfectly fine, while the "random"
one breaks down.
Here is the error traceback:
11: .C64("aplsb1", SIGNATURE = c(SS$signature, SS$signature, "double",
? ? ? ? SS$signature, SS$signature, "double", "double", SS$signature,
? ? ? ? SS$signature, "double", SS$signature, SS$signature, SS$signature,
? ? ? ? SS$signature), nrow, ncol, A at entries, A at colindices, A at rowpointers,
? ? ? ? s, B at entries, B at colindices, B at rowpointers, entries = vector_dc("double",
? ? ? ? ? ? nzmax), colindices = vector_dc(SS$type, nzmax), rowpointers = vector_dc(SS$type,
? ? ? ? ? ? nrow + 1), nzmax + 1, ierr = vector_dc(SS$type, 1), INTENT = c("r",
? ? ? ? ? ? "r", "r", "r", "r", "r", "r", "r", "r", "w", "w", "w",
? ? ? ? ? ? "r", "w"), NAOK = getOption("spam.NAOK"), PACKAGE = SS$package)
10: spam_add(e1, e2, -1)
9: I - lambda * csrw
8: I - lambda * csrw
7: xprodB(lambda, w)
This is where the problem lies:
library(spam)
csrw <- as.spam(wgts)
n <- nrow(wgts)
I <- diag.spam(1, n, n)
lambda <- 0
B <- I - lambda * csrw
is OK. If lambda is NaN:
lambda <- NaN
B <- I - lambda * csrw
Error in .C64("aplsb1", SIGNATURE = c(SS$signature, SS$signature,
"double", :
NAs in argument 7 and 'NAOK = FALSE' (dotCall64)
If lambda is Inf:
lambda <- Inf
B <- I - lambda * csrw
Error in .C64("aplsb1", SIGNATURE = c(SS$signature, SS$signature,
"double", :
NAs in argument 7 and 'NAOK = FALSE' (dotCall64)
and so on. The default method="nlminb" can be changed to "BFGS":
system.time(reg.splm.random <- spml(formula = mx.standard ~ I.score,
data = analysis.subset, listw = listw.wgts, effect = "individual", model =
"random", lag = FALSE, spatial.error = "b", index = c("Group", "Year"),
method="BFGS")
+ )
user system elapsed
4191.258 391.390 4597.724
4597.724/60
[1] 76.62873 taking just over 75 minutes to complete
summary(reg.splm.random)
ML panel with , random effects, spatial error correlation
Call:
spreml(formula = formula, data = data, index = index, w =
listw2mat(listw),
w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl,
method = "BFGS")
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-5.75e-03 -1.08e-03 2.42e-05 -7.50e-06 1.07e-03 8.02e-03
Error variance parameters:
Estimate Std. Error t-value Pr(>|t|)
phi 0.045798 NaN NaN NaN
rho 0.725055 NaN NaN NaN
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.0284e-02 5.0441e-05 203.8892 < 2.2e-16 ***
I.score 5.2369e-03 1.0438e-03 5.0173 5.239e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The code involves inversion of the nxn matrix I - \lambda W in each call
to the log likelihood function, hence the long run time. I've added the
splm maintainer to CC in case this thread didn't yet reach him. I think
that checking that lambda is finite and returning an infinite ll value if
not may work, but the documentation of maxLik() and nlminb would need
checking to see how to signal an invalid lambda.
Hope this clarifies,
Roger
PS. With traceback, it is as simple as setting debug(splm:::semREmod) and
then setting more debug()s - reading the code also helps; patience is
however essential.
6: invSigma(philambda, n, t., w2)
5: objective(.par, ...)
4: nlminb(start = myparms0, objective = ll.c, gradient = NULL, hessian = NULL,
? ? ? ?y = y, X = X, n = n, t. = t., w = w, w2 = w2, scale = parscale,
? ? ? ?control = list(x.tol = x.tol, rel.tol = rel.tol, trace = trace),
? ? ? ?lower = lower.bounds, upper = upper.bounds)
3: est.fun(X, y, ind, tind, n, k, t, nT, w = w, w2 = w2, coef0 = coef0,
? ? ? ?hess = hess, trace = trace, x.tol = x.tol, rel.tol = rel.tol,
? ? ? ?...)
2: spreml(formula = formula, data = data, index = index, w = listw2mat(listw),
? ? ? ?w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl,
? ? ? ?...)
1: spml(formula = mx.standard ~ I.score, data = analysis.subset,
? ? ? ?listw = listw.wgts, effect = "individual", model = "random",
? ? ? ?lag = FALSE, spatial.error = "b", index = c("Group", "Year"))
My sessionInfo:
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8? LC_CTYPE=English_United States.utf8? ? LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C? ? ? ? ? ? ? ? ? ? ? ? ? ?LC_TIME=English_United States.utf8? ??
attached base packages:
[1] stats? ? ?graphics? grDevices utils? ? ?datasets? methods? ?base? ? ?
other attached packages:
?[1] splm_1.6-2? ? ? plm_2.6-2? ? ? ?spdep_1.2-7? ? ?spData_2.2.1? ? sp_1.6-0? ? ? ? forcats_1.0.0? ?stringr_1.5.0? ?dplyr_1.1.0? ? ?purrr_1.0.1? ??
[10] readr_2.1.4? ? ?tidyr_1.3.0? ? ?tibble_3.1.8? ? ggplot2_3.4.1? ?tidyverse_1.3.2 sf_1.0-9? ? ? ?
loaded via a namespace (and not attached):
?[1] nlme_3.1-160? ? ? ? fs_1.6.1? ? ? ? ? ? spatialreg_1.2-6? ? lubridate_1.9.2? ? ?httr_1.4.4? ? ? ? ? tools_4.2.2? ? ? ? ?backports_1.4.1? ??
?[8] utf8_1.2.3? ? ? ? ? R6_2.5.1? ? ? ? ? ? KernSmooth_2.23-20? DBI_1.1.3? ? ? ? ? ?colorspace_2.1-0? ? withr_2.5.0? ? ? ? ?tidyselect_1.2.0? ?
[15] compiler_4.2.2? ? ? cli_3.6.0? ? ? ? ? ?rvest_1.0.3? ? ? ? ?expm_0.999-7? ? ? ? xml2_1.3.3? ? ? ? ? sandwich_3.0-2? ? ? scales_1.2.1? ? ? ?
[22] lmtest_0.9-40? ? ? ?classInt_0.4-8? ? ? proxy_0.4-27? ? ? ? digest_0.6.31? ? ? ?pkgconfig_2.0.3? ? ?dbplyr_2.3.0? ? ? ? collapse_1.9.2? ? ?
[29] ibdreg_0.3.8? ? ? ? rlang_1.0.6? ? ? ? ?readxl_1.4.2? ? ? ? rstudioapi_0.14? ? ?generics_0.1.3? ? ? zoo_1.8-11? ? ? ? ? jsonlite_1.8.4? ? ?
[36] googlesheets4_1.0.1 magrittr_2.0.3? ? ? Formula_1.2-4? ? ? ?s2_1.1.2? ? ? ? ? ? dotCall64_1.0-2? ? ?Matrix_1.5-3? ? ? ? Rcpp_1.0.10? ? ? ??
[43] munsell_0.5.0? ? ? ?fansi_1.0.4? ? ? ? ?lifecycle_1.0.3? ? ?stringi_1.7.8? ? ? ?MASS_7.3-58.1? ? ? ?grid_4.2.2? ? ? ? ? parallel_4.2.2? ? ?
[50] bdsmatrix_1.3-6? ? ?crayon_1.5.2? ? ? ? deldir_1.0-6? ? ? ? lattice_0.20-45? ? ?haven_2.5.1? ? ? ? ?splines_4.2.2? ? ? ?hms_1.1.2? ? ? ? ??
[57] pillar_1.8.1? ? ? ? boot_1.3-28? ? ? ? ?LearnBayes_2.15.1? ?wk_0.7.1? ? ? ? ? ? reprex_2.0.2? ? ? ? glue_1.6.2? ? ? ? ? modelr_0.1.10? ? ??
[64] vctrs_0.5.2? ? ? ? ?spam_2.9-1? ? ? ? ? tzdb_0.3.0? ? ? ? ? Rdpack_2.4? ? ? ? ? miscTools_0.6-26? ? cellranger_1.1.0? ? gtable_0.3.1? ? ? ?
[71] assertthat_0.2.1? ? rbibutils_2.2.13? ? broom_1.0.3? ? ? ? ?e1071_1.7-13? ? ? ? coda_0.19-4? ? ? ? ?class_7.3-20? ? ? ? googledrive_2.0.0??
[78] gargle_1.3.0? ? ? ? units_0.8-1? ? ? ? ?maxLik_1.5-2? ? ? ? timechange_0.2.0? ? ellipsis_0.3.2? ? ?
I tried running the same code on Linux x64 platform as well, with much the same end result.
R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.6 LTS
Here is the link to the minimally reproducible example (code and the actual data):
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dropbox.com%2Fsh%2Fs680xpvfdxowuor%2FAAAH5yF4JtrzWUgdbvQIF6LBa%3Fdl%3D0&data=05%7C01%7Croger.bivand%40nhh.no%7C7b2114f83ed64bb06aad08db0f85b68f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638120839199952772%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=QX%2BWBZkxcnNypTVT%2FmJTZJC1aUmDluawNH3LGOVosTE%3D&reserved=0
Thank you for your help.
Denys Dukhovnov
On Wednesday, February 15, 2023 at 11:12:31 AM EST, Roger Bivand <roger.bivand at nhh.no> wrote:
On Wed, 15 Feb 2023, Denys Dukhovnov via R-sig-Geo wrote:
Dear community,
I am trying to run a spatial panel random effects SEM model using splm
package on a balanced spatial panel of N = 2472 and T = 4, and weights
are based on k = 6 nearest neighbors. I encounter this error, which I
can't find anything about in the documentation or online.
spml(formula = mx.standard ~ I.score,?
? ? ? ? ? ? ? ? ? ? ? ? ?data = analysis.subset,
? ? ? ? ? ? ? ? ? ? ? ? ?listw = listw.wgts,
? ? ? ? ? ? ? ? ? ? ? ? ?effect = "individual",
? ? ? ? ? ? ? ? ? ? ? ? ?model = "random",?
? ? ? ? ? ? ? ? ? ? ? ? ?lag = FALSE,
? ? ? ? ? ? ? ? ? ? ? ? ?spatial.error = "b",
? ? ? ? ? ? ? ? ? ? ? ? ?index = c("Group", "Year"))
Error in .C64("aplsb1", SIGNATURE = c(SS$signature, SS$signature,
"double", :?NAs in argument 7 and 'NAOK = FALSE' (dotCall64)
Please at least provide the output of traceback() after the error (it does not crash, it error-exits correctly). Also provide the output of sessionInfo() - versions of underlying packages may matter (I'm thinking of dotCall64 and packages it uses). .C64() is not called by splm functions directly. Does the error persist on a fully updated system on the same OS? What happens on a different platform if available? Ideally, provide a fully reproducible minimal example using built-in data, or provide your data and the code used to create analysis.subset and listw.wgts, so that the running the code creating the error can be reproduced exactly. The error message reports NAs, are the data complete? Roger
Unlike the same specification but with fixed effects (with model = "within" argument), the one above takes a long while to run until it crashes. As far as I can tell, there are no memory issues. An alternative spatial random effects using spgm() command runs OK, as do non-spatial fixed and random effects models using plm(). For the one above I could use spgm(), but I would prefer to keep all my models in the maximum likelihood estimation for consistency. Any clarifications as to the nature of the error are welcome.? Thank you. Denys Dukhovnov
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7Croger.bivand%40nhh.no%7C7b2114f83ed64bb06aad08db0f85b68f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638120839199952772%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=XEUJSidD5%2FFjLSMxlLy%2BX%2FOijPwhkBOyY0C2KnKl3DI%3D&reserved=0
Roger Bivand Emeritus Professor Department of Economics, Norwegian School of Economics, Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. e-mail: Roger.Bivand at nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en