[R-sig-ME] glmmPQL crashes on inclusion of corSpatial object
Thierry Onkelinx
thierry.onkelinx at inbo.be
Mon Jul 25 22:56:53 CEST 2016
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 op 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 op gmail.com <mailto:patrick.schratz op 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 op gmail.com
> > <mailto:patrick.schratz op 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 op r-project.org
> > <mailto:R-sig-mixed-models op r-project.org> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> >
> > _______________________________________________
> > R-sig-mixed-models op r-project.org
> > <mailto:R-sig-mixed-models op r-project.org> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> >
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list