[R] plotting a 3D poisson surface with persp package

stè snafa snafa at hotmail.it
Fri Apr 18 16:54:18 CEST 2014


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



 
 		 	   		  


More information about the R-help mailing list