RMark Huggins with c=p+offset

[geof]

I have a two-occasion (ch=11, 10, or 01) dataset and I am fitting a Huggins model. Let us say that I have two covariates X and Y, where (at least) Y is continuous. Suppressing the link function here, the model I want to fit looks like this:

p ~ X
c ~ p + Y

In other words, I would like the recapture probability to equal the original capture probability, offset by some adjustment that depends on Y. Maybe denoted p(X)c(p,Y).

My naive attempt is like this:

Code: Select all

a=process.data(mydata,model="Huggins",groups=c("X"),allgroups=T)
b=make.design.data(a)
pform=list(formula=~X,share=F) 
cform=list(formula=~X+Y)
mm=make.mark.model(data=a,ddl=b,parameters=list(p=pform,c=cform))
myfit=run.mark.model(mm,realvcv=T)

However, this appears to allow c to depend on X in a different manner (i.e., different parameters) than how p depends on X. (Right? The output gives different coefficients for X for p vs c). Does anyone know how to force the coefficient(s) for X to be equal in the p and c models while allowing for Y to participate only in the model for c?

Thanks.

Geof

[jlaake]

When c and p are related you use the p formula and share=TRUE. What share does is to merge the design data for the two parameters and adds a field that is 0 for the dominant field (p in this case) and 1 for the subordinate field (c in this case). The name of the 0/1 field is the name of the subordinate parameter. To get what you want you use ~X+c:Y as shown below. That uses X for p and c and only Y for c. --jeff

Code: Select all

a=process.data(mydata,model="Huggins",groups=c("X"),allgroups=T)
b=make.design.data(a)
pform=list(formula=~X+c:Y,share=TRUE) 
mm=make.mark.model(data=a,ddl=b,parameters=list(p=pform))
myfit=run.mark.model(mm,realvcv=T)


[geof]

Jeff,

Thanks, I think that's perfect. I never would have figured out the "c:Y" notation. With apologies for switching to some "real" output,

Code: Select all

Parameter                    Beta       
 -------------------------  -------------
    1:p:(Intercept)         0.7610368   
    2:p:viscatFA            -0.4328788   
    3:p:viscatGO            -0.1468221   
    4:p:c:distdivsd1        1.0573326   
    5:p:c:distdivsd2        0.4254220   
    6:p:c:distdivsd3        -0.2944771   
    7:p:c:distdivsd4        -0.6596013   

let me check my understanding. For viscat=FA and distdivs=d2, the capture probability is invlogit(0.76-0.432) and the recapture probability is invlogit(0.76-0.432+0.425). Is that right?

[jlaake]

I never would have figured out the "c:Y" notation.

It is quite useful for restricting interactions by setting some design data variable x=0/1 and then using x:y to restrict the effect of y where x=1.

let me check my understanding. For viscat=FA and distdivs=d2, the capture probability is invlogit(0.76-0.432) and the recapture probability is invlogit(0.76-0.432+0.425). Is that right?

Yes that is correct.

--jeff

[geof]

It is quite useful for restricting interactions by setting some design data variable x=0/1 and then using x:y to restrict the effect of y where x=1.

Of course. It is the presence of the "c" in "c:Y" that surprised me here, not the general notation for interactions in R.

Anyway, thanks for your help; maybe I'll see you in Panama.

Geof

[jlaake]

Ah! Didn't see it was you. Yes, my parameter sharing approach could use some more documentation. In fact I need a lot of work on the documentation. Even with as much as I've written I need more and it needs to be updated. Won't be in Panama. --jeff

[Snook]

Hi Jeff,

I am going to tag along on this thread if that is ok. I am using a Barker model in a simulation right now. Just setting things up in a simple fashion for now (everything is time and group independent with no covariates).
I've been trying to figure out how to set F=F'.

Following the above example I did:
a <- process.data(BARKER.ch,model="Barker") #Processes input data into Barker form
b <- make.design.data(BARKER.processed)

F.shared<-list(formula=~1,share=TRUE)

mark(a,b,model="Barker",model.parameters=list(F=F.shared,FPrime=F.shared))

What am I doing wrong here that doesn't end up fixing F=F'?


-Andrew

[jlaake]

You aren't doing anything wrong. Only certain parameters have been specified to be shared and these have not been setup to be shared. You can make it work by modifying the values in parameters.txt. You'll find that file under RMark subdirectory where your R library is located. Typically library under your version of R under Program Files.

Replace

Barker F 0 -1 1 Triang all logit ~1
Barker FPrime 0 -1 0 Triang all logit ~1

with

Barker F 0 -1 1 Triang all logit ~1 FALSE Fprime
Barker FPrime 0 -1 0 Triang all logit ~1 F

The first parameter is dominant(primary). When you are sharing, you only give the formula for the dominant parameter and specify share=TRUE - the default is FALSE as specified in parameters.txt.

F.shared<-list(formula=~1,share=TRUE)
mark(a,b,model="Barker",model.parameters=list(F=F.shared)

Let me know if this works ok for you. I'll change the default value in the next release.

--jeff

[jlaake]

It appears that the tabs have been stripped in what I sent. Send me your email address to jeff.laake@noaa.gov and I'll send you a new file.

--jeff

[timkdavies]

Hi there,

I am using the Huggins model (Full het) for analysing small mammal data. When I am setting the formulas to be run in the model, p.dot.shared=list(formula=~1,share=TRUE) is working fine. But when I try to make p and c different but constant I keep on getting outputs with quite ridiculous estimates, SEs and Ci. I presume my code must be incorrect, I have tried several different ways listed below:

Code: Select all

p.dot=list(formula=~1,share=FALSE)
c.dot=list(formula=~1)

and

Code: Select all

p.dot=list(formula=~1)
c.dot=list(formula=~1)

P and c are different but I get estimates of 0.267960E+010, SE- 0.3102223E+011, CI lower - 34895013, 0.2057E+012 (for a group where only 1 individual was caught), for all other groups I got similar values too.

When I run -p.dot.shared=list(formula=~1,share=TRUE), I get much more realistic values - -such as N-hat - 2.24, SE - 0.39, CI lower - 2.008113 and upper - 4.5266676

The models with the crazily high estimates have the lowest AIC value as well. What am I doing wrong?