[R-sig-ME] glmmPQL crashes on inclusion of corSpatial object
Patrick Johann Schratz
patrick.schratz at gmail.com
Mon Jul 25 16:33:48 CEST 2016
Hi Paul,
indeed - thats the reason! What a simple error causing me headaches
since weeks.
Actually I thought I removed all NA´s in the "date" column but did not
check explicitly for it. Caused by a different date format using "/"
instead of "-" for some entries in the data set.
Kind of embarrassing...
Thank you all for your replies guys!!
Am 25.07.16 um 16:33 schrieb Patrick Johann Schratz:
> Hi Paul,
>
> indeed - thats the reason! What a simple error causing me headaches
> since weeks.
>
> Actually I thought I removed all NA´s in the "date" column but did not
> check explicitly for it. Caused by a different date format using "/"
> instead of "-" for some entries in the data set.
>
> Kind of embarrassing...
>
> Thank you all for your replies guys!!
>
>
> Am 25.07.16 um 15:50 schrieb Paul Debes:
>> Hi Patrick,
>>
>> Maybe the NA entries in the random "date" predictor are causing the
>> crash?
>> When I run the model on a subset without NA for "date" it runs.
>> Probably not what you want but it may help at least a bit.
>>
>> Hope this works for you too.
>> Paul
>>
>> str(data)
>> data$hail = as.integer(as.character(data$hail))
>>
>> # removing NA data in "date"
>> data.sub = data[!is.na(data$date),]
>>
>> library("MASS")
>> library("nlme")
>>
>> # changed 'ry' to 'rlat' and 'rx' to 'rlon'
>> correl = corSpatial(value = c(10000, 0.1), form = ~rlat + rlon,
>> nugget = TRUE, fixed = FALSE, type = "exponential")
>> correl = Initialize(correl, data = data.sub)
>>
>> model = glmmPQL(hail ~ 1 + prec_nov_apr + t_min_nov_apr +
>> srad_nov_apr + age,
>> random = ~ 1 | date,
>> data = data.sub,
>> correlation = correl,
>> family = binomial)
>>
>> summary(model)
>> Linear mixed-effects model fit by maximum likelihood
>> Data: data.sub
>> AIC BIC logLik
>> NA NA NA
>>
>> Random effects:
>> Formula: ~1 | date
>> (Intercept) Residual
>> StdDev: 0.7641223 1.042385
>>
>> Correlation Structure: Exponential spatial correlation
>> Formula: ~rlat + rlon | date
>> Parameter estimate(s):
>> range nugget
>> 159.37292936 0.05719671
>> Variance function:
>> Structure: fixed weights
>> Formula: ~invwt
>> Fixed effects: hail ~ 1 + prec_nov_apr + t_min_nov_apr + srad_nov_apr
>> + age
>> Value Std.Error DF t-value p-value
>> (Intercept) -8.686290 1.0712635 1056 -8.108453 0.0000
>> prec_nov_apr 0.043263 0.0057295 1056 7.550986 0.0000
>> t_min_nov_apr 0.197124 0.0772288 1056 2.552470 0.0108
>> srad_nov_apr 0.000183 0.0004113 1056 0.445830 0.6558
>> age 0.002953 0.0053051 1056 0.556672 0.5779
>> Correlation:
>> (Intr) prc_n_ t_mn__ srd_n_
>> prec_nov_apr -0.740
>> t_min_nov_apr 0.021 -0.308
>> srad_nov_apr -0.575 0.162 -0.207
>> age -0.012 -0.095 0.105 -0.074
>>
>> Standardized Within-Group Residuals:
>> Min Q1 Med Q3 Max
>> -1.1957073 -0.4895324 -0.3471985 -0.1436912 12.5749006
>>
>> Number of Observations: 1064
>> Number of Groups: 4
>>
>>
>>
>> On Mon, 25 Jul 2016 16:38:08 +0300, Patrick Johann Schratz
>> <patrick.schratz at gmail.com> wrote:
>>
>>> 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
>>
>>
>
More information about the R-sig-mixed-models
mailing list