Problem with missing interactions?

questions concerning analysis/theory using the R package 'marked'

Problem with missing interactions?

Postby lrmont » Tue May 30, 2017 11:19 am

Hi
I am analysing capture recapture data from a study with bats, with a main purpose to determine the influence of individual asymmetry (Left-Right forearm length) in survival probability. We have 904 marked individuals over a period of two years and 10 capture sessions.

I have been using the package marked for a while with smaller similar data sets and find it excellent. Lately, I have stumbled on some errors (certainly due to my inexperience with the package) and would like to ask for some help. In my study, I found that the model with by far the best fit has the following parameters = Phi=list(formula=~cohort+time+AF+Sex),p=list(formula=~Sex+time). I am having a difficulty predicting values to estimate Phi and p with the following commands:
Code: Select all
> newaf<-expand.grid(Sexo=c("F","M"),AF=seq(0,1,by=0.01))
> mr.asy.pred<-predict(mrasy.crm,newdata=newaf,se=TRUE)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
  contrasts can be applied only to factors with 2 or more levels


I was wondering if the error above has something to do with the fact that there are missing interactions, such as no females in cohort 6, or it could be caused by something else. I get no error if cohort is left out of the model.

A second related question is that it might be better to create bins and group cohorts according to season. I could not find an example in the documentation on how to define the bins and would greatly appreciate any tips that might help.

Thanks a lot for the attention.
Leandro
lrmont
 
Posts: 7
Joined: Tue May 30, 2017 10:00 am
Location: UENF, Rio de Janeiro, Brazil

Re: Problem with missing interactions?

Postby jlaake » Tue May 30, 2017 4:46 pm

Has nothing to do with interactions. The newdata argument you specified is not complete for the model you specified. It should contain all of the variables. Note that you used sexo instead of sex. See example with help for predict.crm.
jlaake
 
Posts: 1413
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Problem with missing interactions?

Postby lrmont » Tue May 30, 2017 10:01 pm

Thanks, I managed to get the prediction to work with the design data, but not with the new data set, even specifying all variables. I was confused because before including cohort, it predicted the values although I did not have to specify time as a variable in the new data.

This is the new error

Code: Select all
>mrasy.crm<-crm(mr.asym.proc,ddl=mr.asym.des,model.parameters=list(Phi=list(formula=~AF+cohort+time+Sexo),p=list(formula=~time+Sexo)),hessian=T)
> newaf<-expand.grid(time=factor(c(1:9)),cohort=factor(c(1:9)),AF=seq(0,1,by=0.1), Sexo=c("F","M"))
> mr.asy.pred<-predict(mrasy.crm,newdata = newaf)
Error in `$<-.data.frame`(`*tmp*`, "Time", value = numeric(0)) :
  replacement has 0 rows, data has 16038
In addition: Warning message:
In min(df$time) : no non-missing arguments to min; returning Inf


Thanks again for the help and sorry for the inconvenience.
lrmont
 
Posts: 7
Joined: Tue May 30, 2017 10:00 am
Location: UENF, Rio de Janeiro, Brazil

Re: Problem with missing interactions?

Postby jlaake » Wed Jun 14, 2017 12:46 pm

It sounded like you had a way forward so haven't jumped on this. It is hard for me to debug this without the data etc that you used. But you can work it out. predict.crm is a fairly short function. Here is the piece of code that is run if you specify a non-null value for newdata.
Code: Select all
      
if(is.data.frame(newdata))
{
      newdata$ch=paste(rep("1",object$data$nocc),collapse="")
      newdata.proc=process.data(newdata,model=object$model,begin.time=object$data$begin.time,
                                groups=names(object$data$group.covariates),accumulate=FALSE)
      ddl=make.design.data(newdata.proc,parameters=object$design.parameters)
          dml=create.dml(ddl,model.parameters=object$model.parameters,
               design.parameters=object$design.parameters)
      object$results$model_data$Phi.dm=dml$Phi$fe
      object$results$model_data$p.dm=dml$p$fe      
}else
   stop("Invalid newdata")
jlaake
 
Posts: 1413
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to analysis help

Who is online

Users browsing this forum: No registered users and 1 guest