[R-sig-ME] glmmPQL crashes on inclusion of corSpatial object
Ben Bolker
bbolker at gmail.com
Mon Jul 25 23:01:43 CEST 2016
that's amazing! thank you.
On 16-07-25 04:56 PM, Thierry Onkelinx wrote:
> Dear Patrick,
>
> I've put a quick and dirty example of your model fit with INLA on
> http://rpubs.com/INBOstats/spde
>
> Best regards,
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to say
> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
> 2016-07-25 15:38 GMT+02:00 Patrick Johann Schratz <patrick.schratz at gmail.com
>> :
>
>> Dear Thierry,
>>
>> thanks for your answer.
>>
>> I also recognized that the spatial correlation structure is causing the
>> problems as the random effects seem to work without any problems when
>> omitting the spatial correlations part.
>> However, trying the formula with a random effect of 2 levels
>> (=evaluation) instead of 4(=date) works:
>>
>> /fit <- glmmPQL(fo, random = ~1 | evaluation, data = d, //
>> // correlation = correl, family = binomial)/
>>
>> Which gives me troubles tracking down the problem. glmmPQL should be
>> able to handle 4 level random effects in combination with spatial
>> correlation structures?
>>
>> Setting "fixed = TRUE" in the corSpatial setup still causes the fitting
>> of glmmPQL() to break in my case (with "date" as random effect). Did you
>> also change other values of the corSpatial setup?
>>
>> Thanks for your link to the R-INLA package. Looks very promising.
>> However it will take me quite some time to adapt my code to R-INLA´s
>> syntax and I doubt I have the time right now.
>>
>> Best regards,
>> Patrick
>>
>>
>>
>> Am 25.07.16 um 13:36 schrieb Thierry Onkelinx:
>>> Dear Patrick,
>>>
>>> It seems like the correlation structure makes the model unstable. I
>>> get convergence when setting fixed = TRUE, but the estimate are very
>>> unstable.
>>>
>>> I'd suggest to have a look at the INLA package which allows to fit
>>> spatially correlated random intercepts. It's described in "Spatial and
>>> Spatio-Temporal Bayesian models with R-INLA" (Blangiardo & Cameletti,
>>> 2015)
>>>
>>> Best regards,
>>>
>>>
>>> ir. Thierry Onkelinx
>>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>>> and Forest
>>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>>> Kliniekstraat 25
>>> 1070 Anderlecht
>>> Belgium
>>>
>>> To call in the statistician after the experiment is done may be no
>>> more than asking him to perform a post-mortem examination: he may be
>>> able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
>>> The plural of anecdote is not data. ~ Roger Brinner
>>> The combination of some data and an aching desire for an answer does
>>> not ensure that a reasonable answer can be extracted from a given body
>>> of data. ~ John Tukey
>>>
>>> 2016-07-25 12:38 GMT+02:00 Patrick Johann Schratz
>>> <patrick.schratz at gmail.com <mailto:patrick.schratz at gmail.com>>:
>>>
>>> MacbookPro:
>>>
>>> > sessionInfo()
>>> R version 3.3.1 (2016-06-21)
>>> Platform: x86_64-apple-darwin15.5.0 (64-bit)
>>> Running under: OS X 10.11.6 (El Capitan)
>>>
>>> locale:
>>> [1] de_DE.UTF-8/de_DE.UTF-8/de_DE.UTF-8/C/de_DE.UTF-8/de_DE.UTF-8
>>>
>>> attached base packages:
>>> [1] stats graphics grDevices utils datasets methods base
>>>
>>> other attached packages:
>>> [1] lme4_1.1-12 Matrix_1.2-6 tibble_1.1 nlme_3.1-128
>>> MASS_7.3-45 gstat_1.1-3 sp_1.2-3
>>>
>>> loaded via a namespace (and not attached):
>>> [1] Rcpp_0.12.6 lattice_0.20-33 intervals_0.15.1 FNN_1.1
>>> spacetime_1.1-5 zoo_1.7-13 assertthat_0.1 grid_3.3.1
>>> pacman_0.4.1 minqa_1.2.4 nloptr_1.0.4
>>> [12] xts_0.9-7 splines_3.3.1 tools_3.3.1
>>>
>>>
>>>
>>> Am 25.07.16 um 12:34 schrieb Phillip Alday:
>>>
>>> Hi Patrick,
>>>
>>> can you send us your sessionInfo()? Knowing the R version is
>>> important, but a few other details matter for debugging this
>>> type of problem!
>>>
>>> Best,
>>> Phillip
>>>
>>>
>>> On 25 Jul 2016, at 19:55, Patrick Johann Schratz
>>> <patrick.schratz at gmail.com
>>> <mailto:patrick.schratz at gmail.com>> wrote:
>>>
>>> Link to data:
>>> <https://www.dropbox.com/s/yi3vf0bmqvydr8h/data.Rd?dl=0>
>>> (1170 obs, 9 variables, .Rd file) [plain link in case sth
>>> goes wrong
>>> with the hyperlink:
>>> https://www.dropbox.com/s/yi3vf0bmqvydr8h/data.Rd?dl=0]
>>>
>>> Simply read it in using `readRDS(file)`.
>>>
>>> I´m trying to setup a GLMM using the `glmmPQL` function
>>> from the `MASS`
>>> package including a random effects part and accounting for
>>> spatial
>>> autocorrelation. However, R (Version: 3.3.1) crashes upon
>>> execution.
>>>
>>> library(nlme)
>>>
>>> # setup model formula
>>> fo <- hail ~ prec_nov_apr + t_min_nov_apr +
>>> srad_nov_apr + age
>>>
>>> # setup corSpatial object
>>> correl = corSpatial(value = c(10000, 0.1), form = ~ry
>>> + rx, nugget
>>> = TRUE,
>>> fixed = FALSE, type = "exponential")
>>> correl = Initialize(correl, data = d)
>>>
>>> # fit model
>>> fit5 <- MASS::glmmPQL(fo, random = ~1 | date, data = d,
>>> correlation = correl, family = binomial)
>>>
>>> What I tried so far:
>>>
>>> - reduce number of observation
>>> - play with `corSpatial` parameters (range and nugget)
>>> - reduce number of fixed predictors
>>> - execute code on Windows, Linux (Debian) and Mac R
>>> installations
>>>
>>>
>>> While I get no error message on my local pc (RStudio just
>>> crashes),
>>> running the script on a server returns the following error
>>> message:
>>>
>>> `R: malloc.c:3540: _int_malloc: Assertion (fwd->size &
>>> 0x4) == 0'
>>> failed. Aborted`
>>>
>>> Debugging leads me to a "glibc" c++ library problem. I
>>> also run valgrind
>>> on it. If you need the output, just ask!
>>>
>>> Help ist highly appreciated!
>>> Cheers, Patrick
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org
>>> <mailto:R-sig-mixed-models at r-project.org> mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org
>>> <mailto:R-sig-mixed-models at r-project.org> mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list