Hi everyone,
I have fitted a poisson GLMM with the fallowing formula:
M1 <- glmer(ID ~ DE + PR +(1 | Plot/Site), data = DATA, family = poisson)
where ID = my dependent variable representing the count of a bird speciesDE =
bush height (independent variable)PR = number of preys (independent variable)(1
| Plot/Site) = my random structure
This is my output
Random effects: Groups Name Variance Std.Dev. Site:Plot (Intercept)
0.5663 0.7525 Plot (Intercept) 0.5778 0.7601 Number of obs: 491,
groups: Site:Plot, 41; Plot, 41
Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept)
0.40556 0.31162 1.301 0.19311 DE -0.14328 0.02647
-5.413 6.19e-08 ***PR 0.03933 0.01133 3.470 0.00052 ***
Since both the covariate (DE and PR) are significant, my idea is to plot
everything in a 3D graph (x = DE; y = PR, z = ID) in persp (i.e. a poisson
plane).
When I try to do that with some simulated data, I obtain the type of graph that
I need, with the poisson surface ad a nice grid on it.
The code that I used is:
x1 <- rnorm (100)x2 <- rnorm (100)abc <- function ( x1 , x2 ) { y <- exp
((1)*((+1 + (-0.2 * x1 + 0.2 * x2)))) } par(bg ="white")x1r <- range ( x1
) x1seq <- seq ( x1r [1], x1r[2], length=50)x2r <- range ( x2 ) x2seq <- seq (
x2r [1], x2r[2], length=50)z <- outer(x1seq, x2seq, abc)
persp (x = x1seq, y = x2seq, z= z, theta =-30, zlim = c(-0.2,10) )
Nevertheless, when I try to do it with my own data, I obtain the poisson
surface, but without the grid! The surface, instead, is colored in black and I
can't understand why.
The code is:
x1 <-DATA$DEx2 <-DATA$PRabc <- function ( x1 , x2 ) { y <- exp
((1)*((+0.4055 + (-0.1432 * x1 + 0.0393 * x2)))) } par(bg ="white")x1r <-
range ( DATA$DE ) x1seq <- seq ( x1r [1], x1r[2], length=491)
x2r <- range ( DATA2$DE) x2seq <- seq ( x2r [1], x2r[2], length=491)z <-
outer(x1seq, x2seq, abc) #I think that the problem is here, because probably I
have to do a prediction... but I can't
use the predict function because I have random factors...
persp (x = x1seq, y = x2seq, z= z, theta =-30, zlim = c(-0.2,10) )
Can someone help me? (I have attached the two graphs)I hope that I have
provided enough details, If not please ask.
Thanks a lot,
Marco
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
