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