Models dropped from model averaging

[Snail_recapture]

Hi,

Newbie here again... apologies

So, I've managed to import my data (with a group and an individual covariate for each record) and run some models. However, I get a LOT of models getting dropped from the model averaging 'because one or more beta variances are not positive', and i'm not entirely sure what the best course of action to take is. I've tried adjust=FALSE, but the output still says the model is dropped... anyway, here is my code from the start so you can see where i'm at. Treatment denotes whether animals were transplanted 'within' or 'between' habitats, and 'PC1' is a shape component. I don't expect PC1 to affect recapture rates, but I expect it to affect survival in the 'between' group. Also, in the tbh.results I removed the deviance column as it messed up the formatting... sorry about that.

Code: Select all

>tbh<-convert.inp("C:/Users/...~ALL_PC1.inp",covariates=c("treatment","PC1"))

>tbh$treatment=factor(tbh$treatment,labels=c("within","between"))

> tbh[1:5,]
                       ch freq treatment     PC1
1 00000000100001000100000   -1    within -0.0269
2 00000000010110001000000   -1    within -0.0383
3 00000000010100001000000   -1    within -0.0052
4 00000000010000000000000    1    within -0.0628
5 00000000010100010000000   -1    within -0.0078

>tbh.model<-mark(tbh)

>tbh.process=process.data(tbh,model="CJS",time.intervals=c(1,1,1,1,1,9,1,1,1,1,9,1,2,1,11,1,3,3,1,27,16,8),groups="treatment")

>tbh.ddl<-make.design.data(tbh.process)

>tbh.model.set<-function(){
  Phi.dot=list(formula=~1)
  Phi.tment=list(formula=~treatment)
  Phi.time=list(formula=~time)
  Phi.time.tment=list(formula=~time*treatment)
  Phi.PC1=list(formula=~PC1)
  Phi.tment.PC1=list(formula=~treatment+PC1)
  p.dot=list(formula=~1)
  p.time=list(formula=~time)
  p.tment=list(formula=~treatment)
  p.time.tment=list(formula=~time*treatment)
  cml=create.model.list("CJS")
  results=mark.wrapper(cml,data=tbh.process,ddl=tbh.ddl,adjust=FALSE)
  return(results)
}

>tbh.results=tbh.model.set()
> tbh.results
                                        model npar     AICc   DeltaAICc       weight
15 Phi(~time * treatment)p(~time * treatment)   44 4768.444   0.0000000 4.799891e-01
11             Phi(~time)p(~time * treatment)   40 4769.194   0.7497838 3.299270e-01
19        Phi(~treatment)p(~time * treatment)   36 4772.084   3.6404737 7.775217e-02
3                 Phi(~1)p(~time * treatment)   35 4772.729   4.2850543 5.633050e-02
23  Phi(~treatment + PC1)p(~time * treatment)   37 4773.656   5.2117553 3.544158e-02
7               Phi(~PC1)p(~time * treatment)   36 4774.745   6.3008737 2.055957e-02
14             Phi(~time * treatment)p(~time)   29 4798.572  30.1278310 1.377387e-07
10                         Phi(~time)p(~time)   23 4810.775  42.3314100 3.083773e-10
18                    Phi(~treatment)p(~time)   19 4811.418  42.9740517 2.236323e-10
22              Phi(~treatment + PC1)p(~time)   20 4812.997  44.5530443 1.015455e-10
2                             Phi(~1)p(~time)   18 4814.885  46.4414397 3.950034e-11
6                           Phi(~PC1)p(~time)   19 4816.862  48.4175517 1.470597e-11
16        Phi(~time * treatment)p(~treatment)   11 5411.708 643.2643514 0.000000e+00
12                    Phi(~time)p(~treatment)    6 5418.927 650.4832649 0.000000e+00
13                Phi(~time * treatment)p(~1)    9 5427.084 658.6401697 0.000000e+00
9                             Phi(~time)p(~1)    4 5430.073 661.6290246 0.000000e+00
20               Phi(~treatment)p(~treatment)    4 5477.326 708.8824246 0.000000e+00
4                        Phi(~1)p(~treatment)    3 5477.346 708.9025037 0.000000e+00
24         Phi(~treatment + PC1)p(~treatment)    5 5478.570 710.1264751 0.000000e+00
8                      Phi(~PC1)p(~treatment)    4 5479.066 710.6215246 0.000000e+00
17                       Phi(~treatment)p(~1)    3 5484.472 716.0277037 0.000000e+00
21                 Phi(~treatment + PC1)p(~1)    4 5485.522 717.0783246 0.000000e+00
1                                Phi(~1)p(~1)    2 5492.590 724.1456027 0.000000e+00
5                              Phi(~PC1)p(~1)    3 5494.497 726.0532037 0.000000e+00

> tbh.model.avg.p=model.average(tbh.results,"p")

Model  2 dropped from the model averaging because one or more beta variances are not positive

Model  3 dropped from the model averaging because one or more beta variances are not positive

Model  6 dropped from the model averaging because one or more beta variances are not positive

Model  10 dropped from the model averaging because one or more beta variances are not positive

Model  11 dropped from the model averaging because one or more beta variances are not positive

Model  14 dropped from the model averaging because one or more beta variances are not positive

Model  15 dropped from the model averaging because one or more beta variances are not positive

Model  18 dropped from the model averaging because one or more beta variances are not positive

Model  19 dropped from the model averaging because one or more beta variances are not positive

Model  23 dropped from the model averaging because one or more beta variances are not positive
Warning messages:
1: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

2: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

3: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

4: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

> tbh.model.avg.Phi=model.average(tbh.results,"Phi")

Model  9 dropped from the model averaging because one or more beta variances are not positive

Model  10 dropped from the model averaging because one or more beta variances are not positive

Model  11 dropped from the model averaging because one or more beta variances are not positive

Model  12 dropped from the model averaging because one or more beta variances are not positive

Model  13 dropped from the model averaging because one or more beta variances are not positive

Model  14 dropped from the model averaging because one or more beta variances are not positive

Model  15 dropped from the model averaging because one or more beta variances are not positive

Model  16 dropped from the model averaging because one or more beta variances are not positive
Warning messages:
1: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

2: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

3: In get.real(model, parameter, data = Covdata, se = TRUE, vcv = vcv) :
Improper V-C matrix for beta estimates. Some variances non-positive.

4: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

5: In get.real(model, parameter, design = model$design.matrix, se = TRUE,  :
 
Improper V-C matrix for beta estimates. Some variances non-positive.

6: In get.real(model, parameter, data = Covdata, se = TRUE, vcv = vcv) :
Improper V-C matrix for beta estimates. Some variances non-positive.


Any help would be grand, I'm pretty new to all this and unfortunately have nobody here in my lab who has used this before! Thanks in advance

[jlaake]

Sophie-

No reason to be embarrassed. It will help if you read the documentation for any function that is giving you problems but in this case if you just type ?model.average you would then need to know to look under model.average.marklist to get help. One of the issues with S3 classes in R. In the help for model.average.marklist you would see an argument drop which you can set to FALSE to include all models.

But that isn't the entire solution and I think you have been misled by my recent posting about adjust. In fact your example is exactly why I chose adjust=TRUE as the default. Notice that when you used model.average for Phi that it is the models with time*treatment for Phi that are being reported as a problem. Undoubtedly you have lots of parameters at boundaries probably phi=1. Assuming my math is correct at 6:40 am the correct parameter count for your "top" model is actually 86 and not 44 which would push it out of contention. If you set adjust=TRUE I believe RMark will give 88 which is high by 2 but much closer. You can use adjust.parameter.count to set them to the correct value. Alternatively you can use the sin link in this case by using Phi.time.tment=list(formula=~-1+time:treatment,link="sin") and likewise for p. In fact you can use the sin link for all of the models you have specified as long as you specify it so it does not have an intercept.
Another thing you might try is to change the time units. Presumably many are at 1 because it is a short time interval. If the time interval is days, you can divide by say 7 to make it a weekly rate or some other time unit. This will shift the parameters away from 1 and may help. Not sure if it is appropriate for your data but you may also want to consider using time+treatment if the temporal pattern might be the same for the treatments. For that model time+treatment you cannot use the sin link.

regards --jeff

[Snail_recapture]

Apologies for the delayed reply.... thanks again Jeff, that's really helpful. I've been doing some reading to try and get my head round things, so hopefully i'm doing things right now.