I would appreciate your help with this very much:

I need to make a 3D graph for a model I created in E-SURGE. My model is:

- Code: Select all
`ad+juv.[i+xind+t*x(1)+t*x(2)+xind.t*x(2)]`

where xind is hatching date, and the two external covariates are year and year squared (standardized).

I want to make a 3D graph with them in R. I am using the command "expression" to create an object that contains the predicted values. But, in expression, do I need to use the standardized values I used to run the model, for the variables y (year) and y2 (year^2), or the original continuous values (y= 0 to 24, and y2= 0....576)?

In any case, the code runs, but the graph makes no sense after back-transforming the predicted values.

I also would like to confirm that I am looking at the right Betas, according to this http://www.phidot.org/forum/viewtopic.php?f=5&t=2346&p=7313&hilit=beta#p7313 the Betas inside transition are in order, please help me to confirm I am identifying the Betas in the right way:

Beta# 1# | +1.579353522 this is adult survival, no problem with this

Beta# 2# | -0.132888241 this should be juvenile survival

Beta# 3# | -0.192818063 t*x(1) or year ?

Beta# 4# | +0.453498834 t*x(1) or year2 ?

Beta# 5# | -0.082595948 marked as ind.cov in IVFV= xind

Beta# 6# | +0.105058405 marked as ind.cov in IVFV = xind.t*x(2) ?

These are the standardized values

- Code: Select all
`yval<-c(-1.6305,-1.4946,-1.3587,-1.2229,-1.087, -0.95111, -0.81524,`

-0.67937,-0.54349,-0.40762, -0.27175,-0.13587,0,0.13587,

0.27175,0.40762,0.54349,0.67937,0.81524,0.95111,1.087,

1.2229, 1.3587, 1.4946, 1.6305)

y2val<-c(-1.116,-1.1008,-1.0754,-1.0399,-0.99428,-0.93848,-0.87253,

-0.79644,-0.7102,-0.61382, -0.50729,-0.39061,-0.26379, -0.12682,0.020291,0.17755,0.34495,0.52251,0.7102,0.90804, 1.116,1.3342,1.5624,1.8009, 2.0494)

this is the code for the graph

- Code: Select all
`i <- 25`

int <- -0.132888241 #intercept juvenile

Bh <- -0.082595948 #beta for day

By <- -0.192818063 #beta for year

By2 <- 0.453498834 #beta for year squared

Bhy2<- 0.105058405 #beta for the interaction

xtemp <- seq(min(-10),max(10),length.out=i) #relative date

xrange <- rep(xtemp,times=i)

ytemp <- seq(min(0),max(24),length.out=i) #year

yrange <- rep(ytemp,each=i)

y2temp <- ytemp^2 #year squared term for nonlinearity

y2range <- rep(y2temp,each=i)

newdata <- data.frame(x=xrange,y=yrange,y2=y2range)

m<-expression(int+Bh*newdata$x+ By*yval + By2*y2val + Bhy2*(newdata$x*y2val))

m1<-eval(m)

m2<-1/(1+exp(-m1)) #transformation to the right scale

library(akima)

xyz<-interp(xrange,yrange,m2)

quartz()

persp(xyz,xlab="date",ylab="year",zlab="survival", cex.lab = 1.5,cex.axis = 1.3,

theta = 45, phi = 25,col="blue", border="grey40", ticktype = "detailed", zlim=c(0,1)

This is my alternative code for expression, using the scale for year 0-24, which isn't working either:

- Code: Select all
`newdata <- data.frame(x=xrange,y=yrange,y2=y2range)`

m<-expression(int+Bh*newdata$x+ By*newdata$y + By2*newdata$y2 + Bhy2*(newdata$x*newdata$y2))

Do you have any idea what am I doing wrong? Thanks a lot for your help!