time-varying individual covariate

[abreton]

From function convert.inp,

Code: Select all

LYC=convert.inp("C:/Users/abreton.CNRDOM/Documents/CSU_Project/MARK_Analyses/LYC_2004_2008/EH_File_LYC_2004_10_Bins.inp",
 group.df=data.frame(age.yr=c("SA04","A04","SA05","A05","SA06","A06","SA07","A07","SA08","A08","SA09","A09","SA10","A10"),
 aclass=rep(c("SA","A"),7),
 syear=c(rep("2004",2),rep("2005",2),rep("2006",2),rep("2007",2),rep("2008",2),rep("2009",2),rep("2010",2))),
 covariates=c("L1","L2","L3","L4","L5","L6","L7","L8","L9","L10","L11","may1Bin","RMiBin"),use.comments=FALSE)

Covariates L1-L11 are fish Lengths on occasions 1-11. Elsewhere, I used a function to predict lengths on ALL passes from available lengths (usually just one). Here's a sample showing the first 5 fish in the dataset,

Code: Select all

> LYC.huggins$data[1:5,]
              ch freq age.yr aclass syear  L1  L2  L3  L4  L5 L6 L7 L8 L9 L10 L11 may1Bin RMiBin group
1:6  00001......    1   SA04     SA  2004 177 177 177 177 182  0  0  0  0   0   0     170    115     8
1:11 00001......    1   SA04     SA  2004 144 144 144 144 150  0  0  0  0   0   0     140    115     8
1:12 00001......    1   SA04     SA  2004 181 181 181 181 186  0  0  0  0   0   0     180    112     8
1:17 00001......    1   SA04     SA  2004 162 162 162 162 168  0  0  0  0   0   0     160    113     8
1:20 00011......    1   SA04     SA  2004 170 170 170 170 170  0  0  0  0   0   0     170    115     8


Note: there were only 5 occasions in this year (group), other years in this dataset had 11 occasions.

In the design matrix in MARK I'm able to specify L1-L11 as time-varying individual covariates. This is accomplished by populating a single beta column with the covariate labels. For example (intercept+Length model),

Code: Select all

B1   B2
1   L1
1   L2
1   L3
1   L4
1   L5
1   L6
1   L7
1   L8
1   L9
1   L10
1   L11
1   L1
1   L2


I'm fitting data to Huggins without mixtures. Thus, L1-L11 will be repeated for groups associated with p. And L2-L11 will be repeated for the same groups associated with c.

In the language of RMark, how do I incorporate L1-L11 as a (single column vector) time-varying individual covariate in the design data? I've searched the forum, read the FM appendix (not enough apparently), and I'm still stumped. Note page C.78 in the Mark book, "Now let’s consider time-varying individual covariates..." I'm still stumped. Thanks for any assistance, and the previously mentioned resources which allowed me to solve every other challenge I encountered. - andre

[jlaake]

Andre-

As long as those times 1-11 match the times in the design data for that field, then you simply use ~L in the formula as shown in Appendix C example (~td) with td1981 to td1986. Your page numbers don't agree with mine but on mine it shows on C-76 how it creates the DM just like you would in MARK. RMark recognizes that there is no covariate named L (or td in the example) so it looks for covariates of the form Lt where t is the time values for that parameter. If it finds a matching covariate for each time then it will construct a column with the appropriate values (L1, L2...L11). If there are any missing covariates it will complain to you. You can also have a time*L interaction as well or any other interaction with the indivdual covariate.

--jeff

[abreton]

I see the text you're referring to (C.76-66), thank you. That's straightforward, I should have focused here rather than text starting on page C.78 (time-varying ind cov).

Your suggestion brought me much closer to the finish line (!), but I'm still doing something wrong.

Code: Select all

# convert
LYC=convert.inp("C:/Users/abreton.CNRDOM/Documents/CSU_Project/MARK_Analyses/LYC_2004_2008/EH_File_LYC_2004_10_Bins.inp",
 group.df=data.frame(age.yr=c("SA04","A04","SA05","A05","SA06","A06","SA07","A07","SA08","A08","SA09","A09","SA10","A10"),
 aclass=rep(c("SA","A"),7),
 syear=c(rep("2004",2),rep("2005",2),rep("2006",2),rep("2007",2),rep("2008",2),rep("2009",2),rep("2010",2))),
 covariates=c("L1","L2","L3","L4","L5","L6","L7","L8","L9","L10","L11","may1Bin","RMiBin"),use.comments=FALSE)

# process
LYC.huggins=process.data(LYC,model="Huggins",begin.time=1,mixtures=1,
 time.intervals=c(1,1,1,1,1,1,1,1,1,1), groups=("age.yr"))

# make design data
LYC.ddl=make.design.data(LYC.huggins)

# specify parameters to fix
p.fix.indices=c(83:88,6:11,93:99,16:22,107:110,30:33,119:121,42:44,129:132,52:55,154,77)
c.fix.indices=c(229:234,159:164,238:244,168:174,251:254,181:184,262:264,192:194,271:274,201:204,275:277,205:207,285:286,
 294,215:216,224)
p.fix.count=rep(0,length(p.fix.indices))
c.fix.count=rep(0,length(c.fix.indices))

# specify pim
p.L.fix.shared=list(formula=~L,share=TRUE,fixed=list(index=c(p.fix.indices,c.fix.indices),value=c(p.fix.count,c.fix.count)))

# fit length model
LYC.huggins.M.L=mark(LYC.huggins,LYC.ddl, model="Length",
 model.parameters=list(p=p.L.fix.shared),output=FALSE, simplify=FALSE)

When I run this I get an error,

Error in extract.mark.output(out, model, adjust, realvcv) :
MARK did not run properly. If error message was not shown, re-run MARK with invisible=FALSE
********Following model failed to run : p(~L)c() **********

andre

[abreton]

When I incorporate just (e.g.) the individual covariate L1 (length on the first occasion) the model runs in RMark and the results match MARK (betas to three decimal places),

Code: Select all

# specify pim
p.L1.shared=list(formula=~L1,share=TRUE)

# fit length model
LYC.huggins.M.L1=mark(LYC.huggins,LYC.ddl, model="Length",
 model.parameters=list(p=p.L1.shared),output=FALSE, simplify=FALSE)

Here are the first 14 rows from the design matrix for p and c,

Code: Select all

> LYC.huggins.M.L1$design.matrix[1:14,]
                p:(Intercept) p:L1
p gA04 aA04 t1              1   L1
p gA04 aA04 t2              1   L1
p gA04 aA04 t3              1   L1
p gA04 aA04 t4              1   L1
p gA04 aA04 t5              1   L1
p gA04 aA04 t6              1   L1
p gA04 aA04 t7              1   L1
p gA04 aA04 t8              1   L1
p gA04 aA04 t9              1   L1
p gA04 aA04 t10             1   L1
p gA04 aA04 t11             1   L1
p gA05 aA05 t1              1   L1
p gA05 aA05 t2              1   L1
p gA05 aA05 t3              1   L1
> LYC.huggins.M.L1$design.matrix[155:168,]
                p:(Intercept) p:L1
c gA04 aA04 t2              1   L1
c gA04 aA04 t3              1   L1
c gA04 aA04 t4              1   L1
c gA04 aA04 t5              1   L1
c gA04 aA04 t6              1   L1
c gA04 aA04 t7              1   L1
c gA04 aA04 t8              1   L1
c gA04 aA04 t9              1   L1
c gA04 aA04 t10             1   L1
c gA04 aA04 t11             1   L1
c gA05 aA05 t2              1   L1
c gA05 aA05 t3              1   L1
c gA05 aA05 t4              1   L1
c gA05 aA05 t5              1   L1

When I replace ~L1 in the pim specification with ~L,

Code: Select all

# specify pim
p.L.shared=list(formula=~L,share=TRUE)

# fit length model
LYC.huggins.M.L=mark(LYC.huggins,LYC.ddl, model="Length",
 model.parameters=list(p=p.L.shared),output=FALSE, simplify=FALSE)

I get the error message that I reported previously. The p structure in the design matrix from the ~L model matches the ~L1 model with the exception that L1 L1 L1(etc) are replaced by L1 L2 L3 (etc). However, the c structure does not, instead L1 L1 L1 (etc) is replaced by LNA LNA LNA (etc). I WONDER IF THIS WHY THE MODEL FAILS TO RUN IN MARK? Here are the first 14 rows for c from the ~L model,

Code: Select all

> LYC.huggins.M.L$design.matrix[155:168,]
                p:(Intercept) p:L
c gA04 aA04 t2              1 LNA
c gA04 aA04 t3              1 LNA
c gA04 aA04 t4              1 LNA
c gA04 aA04 t5              1 LNA
c gA04 aA04 t6              1 LNA
c gA04 aA04 t7              1 LNA
c gA04 aA04 t8              1 LNA
c gA04 aA04 t9              1 LNA
c gA04 aA04 t10             1 LNA
c gA04 aA04 t11             1 LNA
c gA05 aA05 t2              1 LNA
c gA05 aA05 t3              1 LNA
c gA05 aA05 t4              1 LNA
c gA05 aA05 t5              1 LNA

Note this ~L model is the only model I've been unable to run. All the others I've tried work perfectly, no doubt I'm doing something wrong. Thanks for any assistance.

(The option to edit previous posts, add code, etc. are very very helpful. Thanks Evan for the gentle reminder about inserting code.)

*** Update: Jeff resolved this problem on 10-9-11. ***