Nest Survival distribution and covariates help

[mikad]

Hello all,

I am fairly new to Program R, especially RMark, so I apologize if these questions are very basic. Currently, I am working on a nest survival model in RMark following the Mallard example closely. However, I have encountered some issues that I am hoping someone may be able to clarify for me.

First of all, my observed nest predation distribution is not normal or linear, it looks like it would fit a poisson or negative binomial distribution over time. Is there a way, or is it statistically valid, to fit my data to a different distribution (such as the poisson or negative binomial) and re-analyze my models to see if these distributions provide a better fit? I have not been able to find a way to fit my data to these distributions when following the mallard example.

Second, Once my models run, I would like to try and plot them. However, when I use the find.covariates function, 0 covariates are found. When I look at the output to my model, it also states that there are 0 covariates. In my processed data set and ddl, I have type, cam, and road, and time should be one as well. Why does it not recognize these covariates? Does that mean my current model estimates are wrong if the covariates are not being recognized?

Unfortunately, I can not share my code. I can only say that I am very closely following the mallard example. I recognize that this may make assisting me difficult, but can anyone help me out using the mallard example as sample code?

I appreciate any help I can get. Thank you all.

[jlaake]

The nest survival model is a longitudinal analysis with 0/1 outcomes. It is like a std survival analysis except the exact time of failure is not known. It is not fitting a statistical distribution to count data like you are thinking.

Without sharing your code and a snippet of your data there is little I can help with in regards to RMark. I don't understand why you can't share your code and a couple of lines of your data.

--jeff

[mikad]

Hello again,

I don't understand the best way to upload my script and sample data, so here's hoping I did it correctly. I had some trouble understanding the help pages. In any case, here should be my script as well as a snip of the data that I am using.

Thank you for any help you may be able to provide.

Code: Select all

   
id FirstFound LastPresent LastChecked Fate Road_Dist   Type Cam
1    1          1          22          22    0  77.64347   Open   Y
2    8          1          22          22    0  98.71662   Open   N
3   12          1          22          22    0 125.18906   Open   N
4  190          2          21          22    1 101.11220 Active   N
5  237          2          16          17    1  86.69883 Active   Y
6  743          1          17          18    1 111.36745 Active   Y
7  744          1          22          22    0  95.32873 Active   Y
8  758          2          23          23    0  97.84255    Old   Y
9  763          2          22          23    1 181.74645    Old   Y
10 784          2           6           7    1 109.46130    Old   N
11 788          2          10          11    1  60.66460    Old   Y


Code: Select all

nests$id=factor(nests$id)
nests$Cam=factor(nests$Cam);nests
str(nests)

nests.processed <- process.data(nests,model="Nest",nocc=23,groups=c("Type","Cam"),)
nests.ddl <- make.design.data(nests.processed)
nests.ddl$S$Time

run.nests=function()
{
  Dot=mark(nests.processed,nests.ddl)   #the default model is dot
  Dot$results$real
 
  Type=mark(nests.processed,nests.ddl,
            model.parameters=list(S=list(formula=~Type)),) 
 
  TimeType=mark(nests.processed,nests.ddl,
                model.parameters=list(S=list(formula=~Time*Type)),)
 
  TimeTypeInt=mark(nests.processed, nests.ddl,
                   model.parameters=list(S=list(formula=~Time:Type)),)
 
  TimeandType=mark(nests.processed, nests.ddl,
                   model.parameters=list(S=list(formula=~Time+Type)),)
 
  Cam=mark(nests.processed,nests.ddl,
           model.parameters=list(S=list(formula=~Cam)),)
 
  road=mark(nests.processed, nests.ddl,
            model.parameters=list(S=list(formula=~Road_Dist)),)
 
  TimeTrend2=mark(nests.processed, nests.ddl,
                  model.parameters=list(S=list(formula=~Time)),)
 
  TypeRoad=mark(nests.processed, nests.ddl,
                model.parameters=list(S=list(formula=~Road_Dist+Type)),)
 
  Global=mark(nests.processed, nests.ddl,
              model.parameters=list(S=list(formula=~Type+Time+Cam+Road_Dist)),)
 
  return(collect.models())
}

nests.results<-run.nests()
nests.results
#my.table <- collect.models(, adjust = FALSE, table = TRUE)
#model.table(my.table, model.name = FALSE, adjust=FALSE)

avg <- model.average(nests.results,"S",vcv=T)$estimates
names(avg)
avg$estimate

nests.results$Type
nests.results$Type$design.matrix
nests.results$Type$results$beta
nests.results$Type$results$beta.vcv
nests.results$Type$results$real

nests.results$TimeType
nests.results$TimeType$results$beta
nests.results$TimeType$results$beta.vcv

TimeType <- nests.results$TimeType
fc <- find.covariates(TimeType, nests)


[jlaake]

It is important to recognize that RMark has evolved in the nine years since it was first released and some of the examples use an older approach to fitting models and computing predictions.

Also, you only have a single individual covariate which is Road_Dist. I don't know why it couldn't find the covariates but the example you sent me was a model that was not using Road_Dist. Type and Cam are grouping variables. A group variable is a factor variable and for each level a different PIM is created so the estimates for each group have a separate real parameter that is computed automatically.

One final note, that if you were to specify a c-hat, it is not a good idea to extract the results directly from the model because the std errors etc do not reflect the chat. It is safest to use summary function to extract parameters, std errors etc.

What I've done below is to rewrite your code using mark.wrapper which is a cleaner approach and then show how you can use covariate.predictions to compute estimates for a range of individual covariate values which are model averaged across the set of models and then plot values.

Code: Select all

nests.processed <-
process.data(nests,model="Nest",nocc=23,groups=c("Type","Cam"))
nests.ddl <- make.design.data(nests.processed)

run.nests=function()
{
  S.Dot=list(formula=~1)
  S.Type=list(formula=~Type)
  S.TimeType=list(formula=~Time*Type)
  S.TimeTypeInt=list(formula=~Time:Type)
  S.TimeandType=list(formula=~Time+Type)
  S.Cam=list(formula=~Cam)
  S.road=list(formula=~Road_Dist)
  S.TimeTrend2=list(formula=~Time)
  S.TypeRoad=list(formula=~Road_Dist+Type)
  S.Global=list(formula=~Type+Time+Cam+Road_Dist)
  cml=create.model.list("Nest")
  return(mark.wrapper(cml,data=nests.processed,ddl=nests.ddl))
}

nests.results<-run.nests()
nests.results
# compute predictions for different distances from the road for all of the 132 different survival estimates which
# are for 11 values of Road_Dist, for 6 groups and for 22 survival intervals (23 occasions is 22 intervals)
pred <- covariate.predictions(nests.results,data=data.frame(Road_Dist=rep(10*(0:10),each=132),index=rep(1:132,times=11)))$estimates
predwithdata <- merge(pred,nests.ddl$S,by="par.index")
predwithdata <- predwithdata[order(predwithdata$par.index,predwithdata$Road_Dist),]
# you can plot all values as points
plot(predwithdata$Road_Dist,predwithdata$estimate)
# or you can show a line plot with confidence limits, for a single group and time; recognize that depending on the
# other covariates in the models there are potentially 132 different plots.
with(predwithdata[predwithdata$Time==0&predwithdata$group=="ActiveN",],plot(Road_Dist,estimate,ylim=c(0.0,1)))
with(predwithdata[predwithdata$Time==0&predwithdata$group=="ActiveN",],lines(Road_Dist,lcl,lty=2))
with(predwithdata[predwithdata$Time==0&predwithdata$group=="ActiveN",],lines(Road_Dist,ucl,lty=2))




[mikad]

Hello,

Thank you for the advice and more intuitive code. This is certainly not the way that I learned to write code, but it does seem more straight-forward.

[jlaake]

That is simpler question because you get the dsr value directly from summary or from model.average. Your original post suggested that you were interested in covariates.

Here is code that uses the best model and produces a plot over time for each group and then the same thing using the model averaged estimates which use the mean value of road_dist. If the models including road_dist have little weight then the value of road_dist will not matter.

It is often difficult to learn something new, but if you are new to R and RMark you should spend some time with the documentation so you understand the code. There are examples in the RMark documentation that do something similar.

--jeff

Code: Select all

# get  best model for list of results
best=nests.results[[as.numeric(rownames(nests.results$model.table))[1]]]
# extract estimates of S from this model with std errors and conf intervals
s.estimates=summary(best,se=TRUE)$reals$S
# plot values for each of the 6 groups
par(mfrow=c(2,3))
for(type in levels(s.estimates$Type))
  for(cam in levels(s.estimates$Cam))
  {
    with(s.estimates[s.estimates$Type==type & s.estimates$Cam==cam,],
{
  plot(Time,estimate,ylim=c(0,1),main=paste("Type:",type," Cam:",cam),ylab="DSR")
  lines(Time,lcl,lty=2)
  lines(Time,ucl,lty=2)
})
  }


# model average estimates using mean value of road_dist values
s.estimates=model.average(nests.results,parameter="S",vcv=T)$estimates
# plot values for each group
par(mfrow=c(2,3))
for(type in levels(s.estimates$Type))
for(cam in levels(s.estimates$Cam))
{
  with(s.estimates[s.estimates$Type==type & s.estimates$Cam==cam,],
       {
         plot(Time,estimate,ylim=c(0,1),main=paste("Type:",type," Cam:",cam),ylab="DSR")
         lines(Time,lcl,lty=2)
         lines(Time,ucl,lty=2)
       })
}


[mikad]

Hello Jeff,

I see where (some of) my confusion was. I think I just haven't gotten a grasp on the language yet and was using the term "covariate" incorrectly. I was thinking that Type and Cam should be covariates, and that I needed find.covariates to recognize them as such in order to plot my estimates.

I will be spending a lot of time working in R and RMark in the future, so I will certainly take your advice and spend time with the documentation.

Thank you very much for all your help.

[jlaake]

No worries. I've been loose with some of my terminology. Type and Cam are covariates but because they are grouping covariates they are handled differently in MARK than numeric individual covariates which are specified in the design matrix as a variable name like "Road_Dist". find.covariates looks for variable names in the design matrix but it is outdated now and covariate.predictions is the better way to handle predictions for individual numeric covariates. Grouping covariates are in the design data (ddl) because they are attached to a PIM value; whereas individual numeric covariates are in the data but not the design data because their value can differ for a specific PIM value. Your confusion is not just your own as this is not the typical way data are handled in most modelling software but is necessary with the way the MARK interface was written because you have to specify the design matrix by hand. If you were to have to create the design matrix by hand for each animal/nest to handle its unique individual covariate value it would get very cumbersome because the design matrix would be very large. Thus the shortcut by specifying variable names. However, if you use the facilities in R, namely model.matrix function, to create the design matrix via formula then the artificial split between these types of covariates is unnecessary. That is the approach that the mra and marked packages use but neither includes a nest survival model. However, if you look through this phidot list, the web and in the literature you will see that you can fit nest survival models with standard modelling software in R and other stat packages. The nest survival model is like the known fate model where you don't have p - capture probability and there is lots of other software for std survival analysis where p=1. MARK's niche is for problems where p<1 and it must be estimated as well.

--jeff

[mikad]

One more question.

Is there a way using the RMark interface to conduct model averaged parameter estimates for only the top models? In my case, I am looking to get averaged parameter estimates for only my top 3 models.

Thank you for any help you can provide.

[cooch]

One more question.

Is there a way using the RMark interface to conduct model averaged parameter estimates for only the top models? In my case, I am looking to get averaged parameter estimates for only my top 3 models.

Thank you for any help you can provide.

Sure, with a fair bit of work, but it is entirely inappropriate to do so, unless there are clear and defensible structural reasons why models 4 -> N are 'invalid' in some fashion. The only place I've seen this happen with any regularity is with some closed abundance estimators.

Otherwise, you model average over the entire candidate model set. Models 4 -> N will probably have much less weight, and thus will contribute very little to the overall average.