glmmPQL crashes on inclusion of corSpatial object

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]]

    [[alternative HTML version deleted]]