Page 1 of 1

Plotting Nest Survival Models

PostPosted: Thu Jul 11, 2013 9:38 pm
by 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")

Re: Plotting Nest Survival Models

PostPosted: Mon Jul 15, 2013 11:42 am
by bacollier
kdillon wrote: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

Re: Plotting Nest Survival Models

PostPosted: Fri Aug 02, 2013 1:04 pm
by 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)