Plotting Nest Survival Models

[kdillon]

Hi all,

I'm hoping someone might have some tips for me on modeling prediction surfaces for nest survival models. Using the code from the mallard example, I am able to plot the prediction surface for a model with 2 explanatory variables. However, I am struggling to extend it to the case of a model with 3 variables. The code below plots the model: ~fail+elevation, where fail is a 4 level factor and elevation is continuous. I'm trying to plot the model: ~fail+elevation+init, where init is also continuous, plotting the first two variables and holding init constant at its mean value. No adjustments I've made to the code have accomplished this successfully. If anyone has any tips, or other code that has worked for them, it would be much appreciated!

Thank you!

Kristen

Code: Select all

failelev=rfwa2.results$failstage.elevation
fc<-find.covariates(failelev, rfwa2)
seq.elevation <- seq(1799, 2779, by=10)
seq.fail<-c(1,2,3,4)
point <- matrix(nrow=4, ncol=length(seq.elevation))
lower <- matrix(nrow=4, ncol=length(seq.elevation))
upper <- matrix(nrow=4, ncol=length(seq.elevation))
colnum<-0
for (ielevation in seq.elevation) {
  fc$value[1:4]=ielevation
  colnum <- colnum + 1
  fc$value[fc$var=="stagefailed"]=seq.fail
  design=fill.covariates(failelev,fc)
  point[,colnum] <- compute.real(failelev,design=design)[,"estimate"]
  lower[,colnum] <- compute.real(failelev,design=design)[,"lcl"]
  upper[,colnum] <- compute.real(failelev,design=design)[,"ucl"]
}

library(rgl)
open3d()
bg3d("white")
material3d(col="black")
persp3d(seq.fail, seq.elevation, point, aspect=c(1, 1, 0.5), col = "lightblue",
        xlab = "stage", ylab = "elevation", zlab = "DSR", zlim=range(c(upper,lower)),
        main="Daily survival rate")


[bacollier]

Hi all,

I'm hoping someone might have some tips for me on modeling prediction surfaces for nest survival models. Using the code from the mallard example, I am able to plot the prediction surface for a model with 2 explanatory variables. However, I am struggling to extend it to the case of a model with 3 variables. The code below plots the model: ~fail+elevation, where fail is a 4 level factor and elevation is continuous. I'm trying to plot the model: ~fail+elevation+init, where init is also continuous, plotting the first two variables and holding init constant at its mean value. No adjustments I've made to the code have accomplished this successfully. If anyone has any tips, or other code that has worked for them, it would be much appreciated!

Thank you!

Kristen

Code: Select all

failelev=rfwa2.results$failstage.elevation
fc<-find.covariates(failelev, rfwa2)
seq.elevation <- seq(1799, 2779, by=10)
seq.fail<-c(1,2,3,4)
point <- matrix(nrow=4, ncol=length(seq.elevation))
lower <- matrix(nrow=4, ncol=length(seq.elevation))
upper <- matrix(nrow=4, ncol=length(seq.elevation))
colnum<-0
for (ielevation in seq.elevation) {
  fc$value[1:4]=ielevation
  colnum <- colnum + 1
  fc$value[fc$var=="stagefailed"]=seq.fail
  design=fill.covariates(failelev,fc)
  point[,colnum] <- compute.real(failelev,design=design)[,"estimate"]
  lower[,colnum] <- compute.real(failelev,design=design)[,"lcl"]
  upper[,colnum] <- compute.real(failelev,design=design)[,"ucl"]
}

library(rgl)
open3d()
bg3d("white")
material3d(col="black")
persp3d(seq.fail, seq.elevation, point, aspect=c(1, 1, 0.5), col = "lightblue",
        xlab = "stage", ylab = "elevation", zlab = "DSR", zlim=range(c(upper,lower)),
        main="Daily survival rate")


Kristen,
If you don't mind sharing your data/RMark code, there are some different ways to create the plot I think you want, so the trick is finding what works best. If you want, shoot me an email offline <bret@tamu.edu> and I will have a look at it and see what I can come up with.

\bret

[kdillon]

Thanks to Bret for suggesting two simple solutions to this problem. The first solution treats my factor (stagefailed) as a continuous variable for the sake of plotting, as I do in my original code, and plots DSR by stagefailed and elevation, while holding init at it's mean. The second solution uses lattice to plot DSR across all values of elevation and init, creating a separate plot for each level of the factor stagefailed. The code for both solutions is below.

Solution #1:

Code: Select all

elevmodel=rfwa2.results$failinit.elev
fc<-find.covariates(elevmodel,rfwa2)
seq.elevation <- seq(1799, 2779, by=10)
seq.fail<-c(1,2,3,4)
point <- matrix(nrow=4, ncol=length(seq.elevation))
lower <- matrix(nrow=4, ncol=length(seq.elevation))
upper <- matrix(nrow=4, ncol=length(seq.elevation))
colnum<-0
for (ielevation in seq.elevation) {
fc$value=c(rep(143,4), rep(ielevation, 4))
colnum <- colnum + 1
fc$value[fc$var=="stagefailed"]=seq.fail
design=fill.covariates(elevmodel,fc)
point[,colnum] <- compute.real(elevmodel,design=design)[,"estimate"]
lower[,colnum] <- compute.real(elevmodel,design=design)[,"lcl"]
upper[,colnum] <- compute.real(elevmodel,design=design)[,"ucl"]
}

library(rgl)
open3d()
bg3d("white")
material3d(col="black")
persp3d(seq.fail, seq.elevation, point, aspect=c(1, 1, 0.5), col = "lightblue",
xlab = "stage", ylab = "elevation", zlab = "DSR", zlim=range(c(upper,lower)),
main="Daily survival rate")
grid3d(c("x","y+","z"))

Solution #2:

Code: Select all

rfwa2.results$failinit.elev   #get parameter estimates for each variable in model
x=seq(from=1799, to=2779,length=60)     #range of elevation values
y=seq(from=122, to=182, length=60)     #range of init values
xx=expand.grid(x=x,y=y)   
z1=plogis(2.8311847+0.00003389*x-0.0089609*y+0.7481460*0+ 1.2397910*0 +22.44213*0)   #stagefailed1
z2=plogis(2.8311847+0.00003389*x-0.0089609*y+0.7481460*1+ 1.2397910*0 +22.44213*0)   #stagefailed2
z3=plogis(2.8311847+0.00003389*x-0.0089609*y+0.7481460*0+ 1.2397910*1 +22.44213*0)    #stagefailed3
z4=plogis(2.8311847+0.00003389*x-0.0089609*y+0.7481460*0+ 1.2397910*0 +22.44213*1)    #stagefailed4
df1=data.frame(xx,z=z1)
df2=data.frame(xx,z=z2)
df3=data.frame(xx,z=z3)
df4=data.frame(xx,z=z4)
wireframe(z~x*y, data=df1, drape=TRUE)
wireframe(z~x*y, data=df2, drape=TRUE)
wireframe(z~x*y, data=df3, drape=TRUE)
wireframe(z~x*y, data=df4, drape=TRUE)