groups/covariates and model specification

[geof]

Hi everybody.

I'm a novice user of RMark/MARK and have a two questions.

The first is about groups vs. individual covariates. I have
individual records of a simple independent observer survey. The
"catch" histories are either 10, 01, or 11. A continuous covariate is
"estdist". There is also a categorical covariate "simvis" related to
visibility. It is a factor with levels "fine" and "yuck". The
dataframe is "poddat".

Below are three ways I ran the analysis (Huggins model). They are:

1: groups=simvis and simvis not included in formula
2: groups=simvis and simvis included in formula
3: no groups, simvis converted to 0/1 and included in formula

I had expected all three to give the same answer, with #2 being
redundant. Instead #1 gives a different answer. Upon inspection, it
is obvious that simvis is not being fit in this model (e.g., no
parameter estimate and smaller Npar). However, I thought that the
whole point of MARK book section 11.7 was that you could do it either
way. Version #3 should clearly (?) do what I want, so I'm inferring
that #2 is correct. But I don't understand what groups are doing if
#1 isn't matching.

I've included the code and results at the end of this message. I'm
showing the step-by-step functions (rather than the wrapper).


The second question pertains to my overall model. My goal is to
estimate detection probabilities.

The observers are adjacent (but independent), they operate
simultaneously, and sighting has no effect on the animals. (And I
know whether they independently reported a sighting of the same
animal...that's a recap) Thus, I think that (i) I want no effect for
time, (ii) the covariates affect p, and (iii) p=c, requiring
"share=T". Is that the sensible model?

So, to sum up...(a) what is "groups" doing/not doing, and (b) am I
specifying the right model?

Thanks for your help.

Geof

--------------
Version 1: groups=simvis and simvis not in formula

a=process.data(poddat,model="Huggins",groups="simvis")
b=make.design.data(a)
pform=list(formula=~estdist,share=T)
cc=make.mark.model(data=a,ddl=b,parameters=list(p=pform))
d=run.mark.model(cc)
summary.mark(d)

Name : p(~estdist)c()
Npar : 2
-2lnL: 4440.479
AICc : 4444.481
Beta
estimate se lcl ucl
p:(Intercept) 0.4174550000 9.761000e-02 0.2261394000 0.6087707000
p:estdist -0.0002962318 4.473163e-05 -0.0003839058 -0.0002085578
Real Parameter p
Group:simvisfine 0.453644 0.453644
Group:simvisyuck 0.453644 0.453644

---------------
Version 2: groups=simvis and simvis in formula

a=process.data(poddat,model="Huggins",groups="simvis")
b=make.design.data(a)
pform=list(formula=~simvis+estdist,share=T)
cc=make.mark.model(data=a,ddl=b,parameters=list(p=pform))
d=run.mark.model(cc)
summary.mark(d)

Output summary for Huggins model
Name : p(~simvis + estdist)c()
Npar : 3
-2lnL: 4440.012
AICc : 4446.018
Beta
estimate se lcl ucl
p:(Intercept) 0.4582622000 0.1143664000 0.2341040000 0.6824204000
p:simvisyuck -0.0673546000 0.0985758000 -0.2605632000 0.1258540000
p:estdist -0.0002977809 0.0000447766 -0.0003855431 -0.0002100188
Real Parameter p
1 2
Group:simvisfine 0.4629912 0.4629912
Group:simvisyuck 0.4462928 0.4462928

----------------
Version 3: no groups, simvis converted to 0/1 and in formula

poddat$yuckvis=ifelse(poddat$simvis=="yuck",1,0)

a=process.data(poddat,model="Huggins")
b=make.design.data(a)
pform=list(formula=~yuckvis+estdist,share=T)
cc=make.mark.model(data=a,ddl=b,parameters=list(p=pform))
d=run.mark.model(cc)
summary.mark(d)

Output summary for Huggins model
Name : p(~yuckvis + estdist)c()
Npar : 3
-2lnL: 4440.012
AICc : 4446.018
Beta
estimate se lcl ucl
p:(Intercept) 0.4582622000 0.1143664000 0.2341040000 0.6824204000
p:yuckvis -0.0673546000 0.0985758000 -0.2605632000 0.1258540000
p:estdist -0.0002977809 0.0000447766 -0.0003855431 -0.0002100188
Real Parameter p
1 2
0.4536137 0.4536137
Real Parameter c
2
0.4536137

[jlaake]

Specifying groups simply creates the structure that allows you to fit models with those covariates. It is like specification of groups in MARK to create different PIMS for each group. Then you can fit models with or without the group covariate. Typically, you run process.data and make.design.data once and run different models based from the processed data list and design data. So you could have run the first 2 models without repeating the calls to process.data and make.design data.

a=process.data(poddat,model="Huggins",groups="simvis")
b=make.design.data(a)
pform=list(formula=~estdist,share=T)
cc=make.mark.model(data=a,ddl=b,parameters=list(p=pform))
d=run.mark.model(cc)
summary.mark(d)
pform=list(formula=~simvis+estdist,share=T)
cc=make.mark.model(data=a,ddl=b,parameters=list(p=pform))
d=run.mark.model(cc)
summary.mark(d)

Or more simply because mark will call process.data and make.design.data you could have done

pform=list(formula=~estdist,share=T)
c=mark(poddat,model="Huggins",groups="simvis,model.parameters=list(p=pform))
summary(c)
pform=list(formula=~simvis+estdist,share=T)
cc=mark(poddat,model="Huggins",groups="simvis,model.parameters=list(p=pform))
summary(c)

The reason for having the separate steps of process.data and make.design.data is for efficiency when the data set is large and also to enable addition of other design data which cannot be handled with a single call to mark.

With regard to your other questions, time is equivalent to observer in your example and that should be considered in your models as not all observers are the same in their abilities. With regard to p=c, it sounds like you understand how to do that. Apparently MARK will let you fit c for the second occasion but all 3 parameters should not be estimable because you only have 3 bits of data (01,10,11) and you condition on the total number of captures so only 2 are can vary. That is why with 2 occasions you have to assume independence which means that your observers act independently and there is no remaining heterogeneity. If this some type of visual survey, then I suggest that you look at the literature on mark-recapture distance sampling.

--jeff

[geof]

Thanks Jeff.

#1. Got it. Had I reverted to mark() I might not have confused myself!

#2. Yes, I know that time would be equivalent to "observer". Actually "observers" are observation stations, with different observers rotating in/out of each station. I've got covariates for observer name and #observers @ station. I think I want to assume that the only reasons the stations differ are attributable to covariates such as those (and others). So I think I'm happy with what I've got, but yes, I can also devote another parameter if I want and I'll take another look at mrds.

Thanks for the clear and *swift* response!

Geof