www.phidot.org

I have a very basic Rmark code question on how to Fix parameters.
I have 'singular' issues for some parameters because of small populations.
Thus, I want to fix some parameters (survival in this example) to the s. value.

This code, which seems like it matches examples, adds another parameter rather than fixing the one I want. I see in the design matrix that I am not quite doing what I want.

Any help is much appreciated!


Code:

Code: Select all

timint = c(0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,0)

P.PCdotshared = list(formula=~1, share=T) #parameter constant but shared
GPGDP.fixed = list(formula=~1, fixed=0)  #fix parameter to 0

s.time.fixed=list(formula=~time,fixed=list(index=c(6),value=c(0.4)))


tst1=mark(data = rob_m14col, model = "RDHuggins", time.intervals=time.intervals,
 model.parameters=list(S=s.time.fixed,GammaDoublePrime=GPGDPfixed, GammaPrime=GPGDPfixed,p=PCdotshared), filename=paste('tst','3',sep=''))


Code: Select all

From out file:
  INPUT ---    fixed=2;
  INPUT ---        parm(11)=0.4 ;
  INPUT ---        parm(12)=0 ;

  INPUT ---    group=1 S    rows=9 cols=9 Triang ;
  INPUT ---        1 2 3 4 5 11 7 8 9 ;
  INPUT ---        2 3 4 5 6 7 8 9 ;
  INPUT ---        3 4 5 6 7 8 9 ;
  INPUT ---        4 5 6 7 8 9 ;
  INPUT ---        5 6 7 8 9 ;
  INPUT ---        6 7 8 9 ;
  INPUT ---        7 8 9 ;
  INPUT ---        8 9 ;
  INPUT ---        9 ;


And results from out file:

Code: Select all

         LOGIT Link Function Parameters of { S(~time)Gamma''(~1)Gamma'(~1)p(~1)c() }
                                                              95% Confidence Interval
 Parameter                    Beta         Standard Error      Lower           Upper
 -------------------------  --------------  --------------  --------------  --------------
    1:S:(Intercept)         0.0594393       0.5717224       -1.0611366      1.1800151     
    2:S:time2               0.0683505       0.7594010       -1.4200754      1.5567765     
    3:S:time3               -0.3472055      0.6981709       -1.7156205      1.0212096     
    4:S:time4               -1.2493524      0.8087020       -2.8344083      0.3357036     
    5:S:time5               -2.1196401      1.2068735       -4.4851123      0.2458320     
    6:S:time6               -10.916672      1986.0899       -3903.6530      3881.8197     
    7:S:time7               -34.566659      0.3890998E-007  -34.566659      -34.566659   
    8:S:time8               3.0709713       5.4708243       -7.6518445      13.793787     
    9:S:time9               -1.6043899      1.2605026       -4.0749749      0.8661952     
   10:p:(Intercept)         0.2142441       0.0970186       0.0240876       0.4044006     

  Program  MARK  - Survival Rate Estimation with Capture-Recapture Data
   gfortran(Win64) Vers. 8.2 Sep 2017   23-Feb-2019 23:32:42    Page  006
 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 


          Real Function Parameters of { S(~time)Gamma''(~1)Gamma'(~1)p(~1)c() }
                                                               95% Confidence Interval
  Parameter                  Estimate       Standard Error      Lower           Upper
 --------------------------  --------------  --------------  --------------  --------------
     1:S g1 c1 a0 t1         0.5148554       0.1428044       0.2570923       0.7649505                           
     2:S g1 c1 a1 t2         0.5319040       0.1123405       0.3193935       0.7334387                           
     3:S g1 c1 a2 t3         0.4285508       0.0982207       0.2546680       0.6220683                           
     4:S g1 c1 a3 t4         0.2332745       0.1022764       0.0902403       0.4827263                           
     5:S g1 c1 a4 t5         0.1130257       0.1065646       0.0156172       0.5058092                           
     6:S g1 c2 a4 t6         0.1926439E-004  0.0382601       0.1899443E-308  1.0000000                           
     7:S g1 c1 a6 t7         0.1032060E-014  0.5900518E-015  -.1244415E-015  0.2188562E-014                     
     8:S g1 c1 a7 t8         0.9581299       0.2183639       0.5320681E-003  0.9999990                           
     9:S g1 c1 a8 t9         0.1758168       0.1627350       0.0230647       0.6584100                           
    10:p g1 s1 t1            0.5533571       0.0239784       0.5060216       0.5997445                           
    11:S g1 c1 a5 t6         0.4000000       0.0000000       0.4000000       0.4000000       Fixed               
    12:Gamma'' g1 c1 a0 t1   0.0000000       0.0000000       0.0000000       0.0000000       Fixed               


If you are trying to fix survival for that time you will need to add more indices because you have only specified the index for the first cohort. See PIM structure. You are specifying the indices based on the all different PIMs and not the simplified PIMs.

Thanks much for your help on both my questions!

I spent some time trying to figure this out, and got two approaches working (I think), and I get the same results from both-which is a good sign(!).

I am looking for syntax/approach that is easy to put into a function so I can batch run various models. I can explain more background if needed but thought I'd try to keep this simple.

Any suggestions on how to simplify this, or is this the best approach? I am leaning toward the index approach because it follows how I have the rest of the function set up, and seems like I could pass parameters more easily (though I haven't got to that step yet).

Here is an example of a basic model I am running. In this example, I keep everything constant, except allowing S to vary with time.

Code: Select all

fil = convert.inp('C:/data/inpfiles/rob/site1_tst.inp')
time.intervals = c(0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,0) 

PCdotshared = list(formula=~1, share=T) #Mo, p=c, parameter constant but shared
GPGDPfixed = list(formula=~1, fixed=0)  #no temp emigration, fix parameter to 0
Stime= list(formula=~time)  #vary by time

m.stime.pcMo.gnonetst2=mark(data = fil, model = "RDHuggins", time.intervals=time.intervals,  model.parameters=list(S=s.time.fixed,GammaDoublePrime=GPGDPfixed,GammaPrime=GPGDPfixed,p=PCdotshared), filename=paste(filnam,'_stime.pcMo.gnonetst2',sep=''))


Now, I want to fix S for session 6 (to the S.dot value).
This first try uses the indexes.

Code: Select all

#first get design info to know the indices
fil.proc=process.data(fil,model="RDHuggins", time.intervals=time.intervals)
tst.dll=make.design.data(fil.proc)

#create variable
s.time.fixed=list(formula=~time,fixed=list(index=tst.dll$S[tst.dll$S$time==6,c('par.index')],value=c(rep(.4,length(tst.dll$S[tst.dll$S$time==6,c('par.index')])))))

m.stime.pcMo.gnonetst2=mark(data = fil, model = "RDHuggins", time.intervals=time.intervals,
model.parameters=list(S=s.time.fixed,GammaDoublePrime=GPGDPfixed,GammaPrime=GPGDPfixed,p=PCdotshared), filename=paste(filnam,'_stime.pcMo.gnonetst2',sep=''))


This second try modifies the design matrix

Code: Select all

tst2.ddl=tst.dll
tst2.ddl$S$fix=NA
tst2.ddl$S$fix[ddl$S$time==6]=.4
m.stime.pcMo.gnonetst5=mark(data = fil.proc, tst2.ddl, model = "RDHuggins", time.intervals=time.intervals, model.parameters=list(S=Stime, GammaDoublePrime=GPGDPfixed,GammaPrime=GPGDPfixed,p=PCdotshared), filename=paste(filnam,'_stime.pcMo.gnonetst5',sep=''))


And the PIM:

Code: Select all

  INPUT ---    group=1 S    rows=9 cols=9 Triang ;
  INPUT ---        1 2 3 4 5 10 6 7 8 ;
  INPUT ---        2 3 4 5 10 6 7 8 ;
  INPUT ---        3 4 5 10 6 7 8 ;
  INPUT ---        4 5 10 6 7 8 ;
  INPUT ---        5 10 6 7 8 ;
  INPUT ---        10 6 7 8 ;
  INPUT ---        6 7 8 ;
  INPUT ---        7 8 ;
  INPUT ---        8 ;


And the Estimates - time 6 was set to 0.4. (In my example,the problem got kicked to the next
time period. I have very few individuals during this 'era' of sampling. Once I get the method
down, I can play around with the model to try to solve this. All suggestions on general strategy are
welcome, though my specific question at the moment has to do with coding in RMark):

Code: Select all

         Real Function Parameters of { S(~time)Gamma''(~1)Gamma'(~1)p(~1)c() }
                                                               95% Confidence Interval
  Parameter                  Estimate       Standard Error      Lower           Upper
 --------------------------  --------------  --------------  --------------  --------------
     1:S g1 c1 a0 t1         0.5148332       0.1428026       0.2570786       0.7649315                           
     2:S g1 c1 a1 t2         0.5319099       0.1123405       0.3193982       0.7334436                           
     3:S g1 c1 a2 t3         0.4285512       0.0982204       0.2546689       0.6220679                           
     4:S g1 c1 a3 t4         0.2332757       0.1022768       0.0902409       0.4827280                           
     5:S g1 c1 a4 t5         0.1124529       0.1060254       0.0155498       0.5040460                           
     6:S g1 c1 a6 t7         0.1101896E-005  0.0025814       0.1086435E-309  1.0000000                           
     7:S g1 c1 a7 t8         0.9581136       0.2183748       0.5336837E-003  0.9999990                           
     8:S g1 c1 a8 t9         0.1758185       0.1627367       0.0230648       0.6584140                           
     9:p g1 s1 t1            0.5533868       0.0239753       0.5060572       0.5997679                           
    10:S g1 c1 a5 t6         0.4000000       0.0000000       0.4000000       0.4000000       Fixed               
    11:Gamma'' g1 c1 a0 t1   0.0000000       0.0000000       0.0000000       0.0000000       Fixed               


Any suggestions on how to simplify this?

Thanks very much!

Both the index approach and adding fix to the design data. The latter is best when you want the parameters to be fixed for all models and the former for just specific models. For your case, the index approach is the better one. However, I'm questioning your logic of fixing it to some value like 0.4. I understand you have sparse data but what about binning time so the sparse year is held constant with a consecutive year(s). Or use an environmental covariate. Anything seems better than plugging in some number you think is reasonable.

With regard to using a function, I'm not sure what you want to do. There are lots of examples in RMark help that show how to use create.model.list and mark.wrapper to run a series of models and return a marklist which has each model and a model selection table. Also unless you have a real good reason to name individual files, that is not something I would choose to do. The output file name is stored in the mark object.

Thank you very much for your prompt reply!

Regarding what value to fix the parameters at, I would welcome any advice. These data have years when no animals were found or only 1 or 2, and so I am getting lots of singular values and some non-convergence (Va09 error). Rather than just tossing out these models, one recommendation was to try fixing these years at the value for the constant model. Thus, 0.4 is S for the S.constant model. And where p is singular, I would fix that at the p.constant value. I suppose I could use the value from the adjacent year. An environment covariate would be best and for some situations, I may be able to do this, but not all. I tried using initial as explained in the RMark manual, but couldn't get it to work, so hard-coded it. (I am still a relative RMark newbie and R is not my most proficient programming language, but I am getting there.)

Re the coding: A while ago, I built functions based off examples I pulled from RMark help etc. that run a series of models, average them, and format output. Not necessarily the most elegant but it works. I pass parameters like datasets or other information. I use this 'stock ' code for various projects, tweaking as needed. I haven't quite got there but I can see possibly wanting to pass which indices to fix or other information. To date my needs have been pretty simple, but now I need to do more complicated modeling. I hadn't known about create.model.wrapper until this go round of analysis. With a limited amount of time, I couldn't get it to work, and since I already had code written and with a lot of other stuff to learn, I took the path of least resistance. I will pursue it eventually, it looks very useful.

Thanks again! This forum is really helpful!