RDOccupEG problem

[Buco]

Hi!

I try to run this function from the RDOccupancy help file (the original help file use the RDOccupPG model, so there is no Epsilon):

data(RDOccupancy)
run.RDExample=function()
{
Psi.area=list(formula=~samplearea)
Psi.cover=list(formula=~cover)
Psi.areabycover=list(formula=~samplearea*cover)
Psi.dot=list(formula=~1)
Psi.time=list(formula=~time)
p.dot=list(formula=~1)
p.occ=list(formula=~occ)
p.area=list(formula=~sample.area)
p.coverbyocc=list(formula=~occ*cover)
gam.area=list(formula=~samplearea)
epsilon.area=list(formula=~samplearea)
gam.dot=list(formula=~1)
epsilon.dot=list(formula=~1)
time_intervals=c(0,0,1,0,0,1,0,0)
RD_process=process.data(RDOccupancy, model="RDOccupEG", time.intervals=time_intervals)
RD_ddl=make.design.data(RD_process)
model.p.dot.Psi.dot.gam.dot.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
model.p.dot.Psi.time.gam.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.time, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
model.p.dot.Psi.area.gam.area<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.area, Gamma=gam.area), invisible=TRUE)
model.p.dot.Psi.cover.gam.area<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.cover, Gamma=gam.area), invisible=TRUE)
model.p.occ.Psi.dot.gam.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.occ, Psi=Psi.dot, Gamma=gam.dot), invisible=TRUE)
model.p.occ.Psi.area.gam.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.occ, Psi=Psi.area, Gamma=gam.dot), invisible=TRUE)
return(collect.models())
}
robustexample<-run.RDExample()

but when i hit enter, i get this message:

Error in `contrasts<-`(`*tmp*`, value = "contr.treatment") :
contrasts can be applied only to factors with 2 or more levels

Sorry but that, I do not understand what does it mean, or how can i fix it! Any idea?

Thanks in advance for the answers!
Buco

[bacollier]

Buco,
You have a few places where your variables are not named correctly (e.g., p.area=list(formula=~sample.area)--there is no period (.) in samplearea in the column name so R will not read it correctly for use in MARK) and you also don't have Epsilon defined in all your models, which is necessary I think for this models structure.

I debugged your error down to the second model call (the one with Psi=Psi.time in it). All the others seem to work just fine on my laptop (see cleaned up version of your code below). I suspect the error you are seeing has something to do with the fact that estimates for some of the robust design occupancy models parameters are actually predicted based on the colonization/extinction estimates and an initial occupancy estimate for time 1, and you are trying to fit a model structure which predicts occupancy for each time period (which would mean you do not need both colonization and extinction estimates). But, I am not in my office right now so I cannot reference this models structure from the literature so I am kind of guessing, but I bet I am close.

Try running the below on your system (note I commented out the 2nd model call which tossed the error). This code works on mine. The epsilon parameter is not estimated worth a damn in any of the models (probably something to do with how I simulated the dataset, will add it to my to fix list) but the models run so that should help you out.

Bret

data(RDOccupancy)
run.RDExample=function()
{
Psi.area=list(formula=~samplearea)
Psi.cover=list(formula=~cover)
Psi.areabycover=list(formula=~samplearea*cover)
Psi.dot=list(formula=~1)
Psi.time=list(formula=~time)
p.dot=list(formula=~1)
p.occ=list(formula=~occ)
p.area=list(formula=~samplearea)
p.coverbyocc=list(formula=~occ*cover)
gam.area=list(formula=~samplearea)
epsilon.area=list(formula=~samplearea)
gam.dot=list(formula=~1)
epsilon.dot=list(formula=~1)
time_intervals=c(0,0,1,0,0,1,0,0)
RD_process=process.data(RDOccupancy, model="RDOccupEG", time.intervals=time_intervals)
RD_ddl=make.design.data(RD_process)
M1<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
#M2<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.time, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=FALSE)
M3<-mark(RD_process, RD_ddl,model.parameters=list(p=p.dot, Psi=Psi.area, Gamma=gam.area, Epsilon=epsilon.dot),invisible=TRUE)
M4<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.cover, Gamma=gam.area, Epsilon=epsilon.dot),invisible=TRUE)
M5<-mark(RD_process, RD_ddl, model.parameters=list(p=p.occ, Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
M6<-mark(RD_process, RD_ddl, model.parameters=list(p=p.occ, Psi=Psi.area, Gamma=gam.dot, Epsilon=epsilon.dot),invisible=TRUE)
return(collect.models())
}
robustexample<-run.RDExample()

[jlaake]

The sample.area instead of samplearea is a problem in the original example, but not relevant to the question. I don't know these models well but it appears that Psi cannot vary by time in the EG model. Thus, Brett's conclusion is correct that it is the Psi.time model that caused the problem. What Brett did is what you should do in debugging problems. The next step would be to look at RD_ddl$Psi and you would see that it has only one record. The error message you got

Error in `contrasts<-`(`*tmp*`, value = "contr.treatment") :
contrasts can be applied only to factors with 2 or more levels

was saying that it was trying to create contrasts for a factor level that only had one level. That occurred because RD_ddl$Psi has only one record and thus one value for time. Thus you cannot create a ~time model. It has been awhile since I put these models together but I think this is correct and there is not supposed to be time variation in Psi in this model.

--jeff

[Buco]

Hi,

First of all, I would like to say thank you, for the quick answers!
Now I understand what was the problem!

I have an other question:
How is it that I get another results in R, as in Mark?


Thanks!

[bacollier]

Buco,
RMark is a front end for MARK, not a stand alone program. If you are getting difference answers from RMark and MARK, then you have probably set something up differently between the 2 as MARK is doing all the analysis in the background in both cases.

Bret