how to modavgpred with obsCov

278 views
Skip to first unread message

Christina

unread,
Feb 23, 2021, 6:30:13 PM2/23/21
to unmarked

I am trying to model average the predictions for Model 3. I keep getting the error below and I don’t know how to fix it. coy_cov is an obsCov with 1 ,0 or NA for each survey. Am I missing something obvious? Thank you

 

M<-42

J<-246

T <- 3

 

season<- matrix(c('1','2','3'), nrow(raccoon), T, byrow=TRUE)

coy_cov<-coyoteMS.data[,2:247]

coy_cov[is.na(coy_cov) != is.na(coy_cov)]<-NA

hcov<-humans[,2:247]

hcov[is.na(hcov) != is.na(hcov)]<-NA

raccoonMS.umf <- unmarkedMultFrame(y=raccoon[,2:247],

                                   siteCovs = NULL,

                                   obsCovs= list(coy=coy_cov,hcov=hcov),

                                   yearlySiteCovs=list(season=season),

                                   numPrimary = T)

 

Cand.models<-list()

Cand.models[[1]]<-colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~season-1, raccoonMS.umf, se=TRUE) 

Cand.models[[2]]<-colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~coy, raccoonMS.umf, se=TRUE) 

Cand.models[[3]]<-colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~season-1+ coy, raccoonMS.umf, se=TRUE) 

Cand.models[[4]]<-colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~hcov, raccoonMS.umf, se=TRUE) 

 

Modnames <- paste("mod", 1:length(Cand.models), sep = "")

 

new.dat<-data.frame(season=factor(season,levels=unique(season)),coy=coy_cov)

 

modavgPred(cand.set = Cand.models, modnames = Modnames,newdata='new.dat',parm.type='detect')

Error in if (ncolumns == 1) newdata$blank.fake.column.NAs <- NA :

  argument is of length zero 

Marc J. Mazerolle

unread,
Feb 23, 2021, 7:27:25 PM2/23/21
to unma...@googlegroups.com
Hi Christina,

You should not put the name of your new.dat object between quotes. You
just need to supply the object. The following should work:

modavgPred(cand.set = Cand.models, modnames = Modnames,
newdata = new.dat, parm.type = 'detect')

Best,

Marc
--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

Veuillez noter que je suis en télétravail. Le meilleur moyen de me rejoindre est par courriel.

-------- Message initial --------
De: Christina <christina.a...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: [unmarked] how to modavgpred with obsCov
Date: Tue, 23 Feb 2021 15:30:13 -0800

[Externe UL*]

Christina

unread,
Mar 6, 2021, 5:59:12 PM3/6/21
to unmarked
Thanks for your response. However, when I do that I get a different error shown below that I also have not been able to figure out

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels 

Christina

Marc J. Mazerolle

unread,
Mar 8, 2021, 7:45:55 AM3/8/21
to unma...@googlegroups.com
Christina,

You should provide a minimum reproducible example so we can try to
debug the problem. At least, what does str( ) give you for coy_cov,
season, and new.dat? The variables specified in new.dat must be of the
same type as was used in your analysis. If you simply want a prediction
for each year for a given value of coy_cov, you need to supply a single
value of coy_cov (i.e., new.dat should only have 4 rows). If you want
predictions across a range of coy_cov values for each year, you should
create a vector of coy_cov values of interest to you and replicate this
vector for each year (e.g., expand.grid( ) ).

Best,

Marc
--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

Veuillez noter que je suis en télétravail. Le meilleur moyen de me rejoindre est par courriel.

-------- Message initial --------
De: Christina <christina.a...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: Re: [unmarked] how to modavgpred with obsCov
Date: Sat, 06 Mar 2021 14:59:12 -0800

[Externe UL*]

Christina

unread,
Mar 9, 2021, 10:55:44 AM3/9/21
to unmarked
Marc, 

I am looking for the latter. I have also already tried for the first season (which includes 82 surveys)

coy1<-coy_cov[,1:82] 
coy1_DF<-as.data.frame(coy1) 
new.dat<-data.frame(season=season[1],coy='coy1_DF') 
modavgPred(cand.set = Cand.models, modnames = Modnames,newdata=new.dat,parm.type='detect') 

Which gives the error 
Error in `[[<-.data.frame`(`*tmp*`, nm, value = c(1L, 1L, 1L, 1L, 1L,  : 
  replacement has 126 rows, data has 42 

I guess I am really just stuck on how to reformat this properly. Even when I tried this 

coy1<-cbind(season[1],coy1) 

It did not work. 

I can see that there is a mismatch in rows but I cannot figure out where or how. My guess would be that even though I am separating out one season it is still trying to run with all 3 seasons. 

str() for when I tried all seasons gave the following: 

str(coy_cov)
tibble [42 x 246] (S3: tbl_df/tbl/data.frame) 

str(season)
 chr [1:42, 1:3] 

str(new.dat)
'data.frame': 42 obs. of  4 variables:
 $ season.V1: Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
 $ season.V2: Factor w/ 1 level "2": 1 1 1 1 1 1 1 1 1 1 ...
 $ season.V3: Factor w/ 1 level "3": 1 1 1 1 1 1 1 1 1 1 ...
 $ coy      : chr  "coy_cov" "coy_cov" "coy_cov" "coy_cov" ...


for one season new.dat it gave 

 str(new.dat)
'data.frame': 42 obs. of  2 variables:
 $ V1 : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
 $ coy: chr  "coy1_DF" "coy1_DF" "coy1_DF" "coy1_DF" ...  

Happy to send data outside of the group if that helps. 

Thank you, 
Christina 




Marc J. Mazerolle

unread,
Mar 9, 2021, 1:02:53 PM3/9/21
to unma...@googlegroups.com
Christina,

Send me offlist your data and models saved in an R object file, with
save( ).

Marc
--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

Veuillez noter que je suis en télétravail. Le meilleur moyen de me rejoindre est par courriel.

-------- Message initial --------
De: Christina <christina.a...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: Re: [unmarked] how to modavgpred with obsCov
Date: Tue, 09 Mar 2021 07:55:44 -0800

[Externe UL*]

Marc J. Mazerolle

unread,
Mar 12, 2021, 8:47:07 AM3/12/21
to unma...@googlegroups.com
For those interested, the issue with Christina's error was related to
the mis-specification of the data frame and exclusion of a detection
variable that appeared in some models.

The following resolved the issue:

##produce all combinations of season and coyote detection
> new.dat <- expand.grid(season = factor(c("1", "2", "3")),
coy = c(0, 1),
hcov = 1)
> new.dat
season coy hcov
1 1 0 1
2 2 0 1
3 3 0 1
4 1 1 1
5 2 1 1
6 3 1 1

> modavgPred(cand.set = Cand.models, modnames = Modnames,
newdata = new.dat, parm.type = 'detect')

Model-averaged predictions on the response scale
based on entire model set and 95% confidence interval:

mod.avg.pred uncond.se lower.CL upper.CL
1 0.295 0.010 0.276 0.314
2 0.289 0.012 0.265 0.313
3 0.302 0.013 0.277 0.329
4 0.222 0.034 0.163 0.293
5 0.217 0.034 0.158 0.289
6 0.228 0.035 0.168 0.301


Marc
--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

Veuillez noter que je suis en télétravail. Le meilleur moyen de me rejoindre est par courriel.

-------- Message initial --------
De: Christina <christina.a...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: Re: [unmarked] how to modavgpred with obsCov
Date: Tue, 09 Mar 2021 07:55:44 -0800

[Externe UL*]

Tanvi Gurjar

unread,
Mar 13, 2024, 12:27:31 PM3/13/24
to unmarked
Hi Marc,

I am trying to get averaged prediction with AICcmodavg basis c-hat values between 1.9 and 2.1 for my top models for Lissotriton vulgaris occupancy. When I try to use "modavgPred", I am get the same error as Christina. The dataframe I have specified for 'modavgPred' has the same covariates as the dataframe used to run my candidate models. I would be very grateful for any input on this.
Please find the code below:

#specifying a dataframe, same as the original
newdata <- unmarkedFrameOccu(y=BM, siteCovs = siteCovs, obsCovs = obsCovs)
#Creating a list containing the candidate models
Cand.list <- list(BM.m26, BM.m27, BM.m30)

#input-ing this into the "modavgPred" function:
modavgPred(cand.set = Cand.list, modnames = NULL, newdata=newdata, parm.type = 'psi', c.hat = 2)

#Alternatively, I tried giving specifying the range for covariates used in the respective models:
modavgCustom( cand.set = Cand.list, modnames = NULL, newdata=newdata,  (years=seq(-1.04,3.72, length=20),size = seq(-0.55,2.84, length=20),water.depth=seq(-1.89,1.22, length=20), FISH = factor("A", levels=c("P","A")), parm.type = "psi", c.hat = 2, type="response"), parm.type = 'psi', c.hat=2)

The error message: 
Error in if (ncolumns == 1) newdata$blank.fake.column.NAs <- NA : argument is of length zero

Here is the original data frame:
BM3 <- unmarkedFrameOccu (y = BM, siteCovs = siteCovs, obsCovs = obsCovs)

And the individual models:
BM.m26 <- occu(~1 ~years+size+FISH, data=BM3) 
BM.m27 <- occu(~1 ~years+water.depth+FISH, data=BM3) 
BM.m30 <- occu(~1 ~years+FISH, data=BM3) 

Here is the "newdata" dataframe:
> str(newdata) Formal class 'unmarkedFrameOccu' [package "unmarked"] with 5 slots ..@ y : num [1:88, 1:5] 0 0 0 0 0 0 0 0 0 0 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:5] "o1" "o2" "o3" "o4" ... ..@ obsCovs :'data.frame': 440 obs. of 3 variables: .. ..$ MONTH : Factor w/ 5 levels "EaAu","EAp","EJ",..: 2 4 3 5 1 2 4 3 5 1 ... .. ..$ ACCESS: Factor w/ 3 levels "EASY","HARD",..: 1 1 1 1 1 1 1 1 1 1 ... .. ..$ HYDROP: Factor w/ 2 levels "DRY","WET": 2 2 2 2 2 2 2 2 2 2 ... ..@ siteCovs: tibble [88 × 21] (S3: tbl_df/tbl/data.frame) .. ..$ MANAGEMENT : num [1:88] 2 2 3 2 2 3 2 2 2 3 ... .. .. ..- attr(*, "levels")= chr [1:3] "1" "2" "3" .. ..$ AGE : Factor w/ 2 levels "Alt","Jung": 2 2 2 2 2 1 1 1 1 2 ... .. ..$ SUBSTRATE : Factor w/ 4 levels "Gra","Gyp","Lim",..: 4 4 4 4 4 4 4 4 4 4 ... .. ..$ FISH : Factor w/ 2 levels "A","P": 1 1 2 2 2 2 2 2 2 2 ... .. ..$ RACOON : Factor w/ 2 levels "A","P": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ LATITUDE : num [1:88] 50.9 50.9 50.9 50.9 50.9 ... .. ..$ LONGITUDE : num [1:88] 6.77 6.77 6.77 6.78 6.78 ... .. ..$ size : num [1:88, 1] -0.551 -0.551 -0.551 -0.551 -0.551 ... .. .. ..- attr(*, "scaled:center")= num 219823 .. .. ..- attr(*, "scaled:scale")= num 398721 .. ..$ years : num [1:88, 1] -0.572 -0.954 -0.858 -0.189 -0.189 ... .. .. ..- attr(*, "scaled:center")= num 11 .. .. ..- attr(*, "scaled:scale")= num 10.5 .. ..$ open.water : num [1:88, 1] -2.83 0.705 -1.527 -1.155 -1.155 ... .. .. ..- attr(*, "scaled:center")= num 81.1 .. .. ..- attr(*, "scaled:scale")= num 26.9 .. ..$ woods : num [1:88, 1] -0.745 -0.593 -0.782 -0.745 -0.593 ... .. .. ..- attr(*, "scaled:center")= num 20.6 .. .. ..- attr(*, "scaled:scale")= num 26.4 .. ..$ emerse.veg : num [1:88, 1] 2.996 -0.702 1.634 1.244 1.244 ... .. .. ..- attr(*, "scaled:center")= num 18 .. .. ..- attr(*, "scaled:scale")= num 25.7 .. ..$ reed : num [1:88, 1] -0.72 -0.72 2.41 1.79 2.25 ... .. .. ..- attr(*, "scaled:center")= num 23 .. .. ..- attr(*, "scaled:scale")= num 31.9 .. ..$ rushes : num [1:88, 1] 3.199 -0.858 -0.858 -1.072 -0.858 ... .. .. ..- attr(*, "scaled:center")= num 25.1 .. .. ..- attr(*, "scaled:scale")= num 23.4 .. ..$ avg.veg.ht : num [1:88, 1] -0.5254 -0.7554 0.5097 -0.2953 -0.0653 ... .. .. ..- attr(*, "scaled:center")= num 106 .. .. ..- attr(*, "scaled:scale")= num 87 .. ..$ total.area : num [1:88, 1] 0.167 0.167 0.167 0.167 0.167 ... .. .. ..- attr(*, "scaled:center")= num 150 .. .. ..- attr(*, "scaled:scale")= num 85.7 .. ..$ ac.area : num [1:88, 1] -0.0895 -0.0895 -0.0895 -0.0895 -0.0895 ... .. .. ..- attr(*, "scaled:center")= num 42.4 .. .. ..- attr(*, "scaled:scale")= num 16.9 .. ..$ inac.area : num [1:88, 1] 0.0895 0.0895 0.0895 0.0895 0.0895 ... .. .. ..- attr(*, "scaled:center")= num 57.6 .. .. ..- attr(*, "scaled:scale")= num 16.9 .. ..$ long_temp : num [1:88, 1] 1.09 1.09 1.09 1.09 1.09 ... .. .. ..- attr(*, "scaled:center")= num 10.1 .. .. ..- attr(*, "scaled:scale")= num 0.597 .. ..$ long_prec : num [1:88, 1] 0.485 0.485 0.485 0.485 0.485 ... .. .. ..- attr(*, "scaled:center")= num 721 .. .. ..- attr(*, "scaled:scale")= num 94.8 .. ..$ water.depth: num [1:88, 1] 0.485 0.485 0.485 0.485 0.485 ... .. .. ..- attr(*, "scaled:center")= num 721 .. .. ..- attr(*, "scaled:scale")= num 94.8 ..@ mapInfo : NULL ..@ obsToY : num [1:5, 1:5] 1 0 0 0 0 0 1 0 0 0 ...

And, here is the original dataframe: > str(BM3) Formal class 'unmarkedFrameOccu' [package "unmarked"] with 5 slots ..@ y : num [1:88, 1:5] 0 0 0 0 0 0 0 0 0 0 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:5] "o1" "o2" "o3" "o4" ... ..@ obsCovs :'data.frame': 440 obs. of 3 variables: .. ..$ MONTH : Factor w/ 5 levels "EaAu","EAp","EJ",..: 2 4 3 5 1 2 4 3 5 1 ... .. ..$ ACCESS: Factor w/ 3 levels "EASY","HARD",..: 1 1 1 1 1 1 1 1 1 1 ... .. ..$ HYDROP: Factor w/ 2 levels "DRY","WET": 2 2 2 2 2 2 2 2 2 2 ... ..@ siteCovs: tibble [88 × 21] (S3: tbl_df/tbl/data.frame) .. ..$ MANAGEMENT : num [1:88] 2 2 3 2 2 3 2 2 2 3 ... .. .. ..- attr(*, "levels")= chr [1:3] "1" "2" "3" .. ..$ AGE : Factor w/ 2 levels "Alt","Jung": 2 2 2 2 2 1 1 1 1 2 ... .. ..$ SUBSTRATE : Factor w/ 4 levels "Gra","Gyp","Lim",..: 4 4 4 4 4 4 4 4 4 4 ... .. ..$ FISH : Factor w/ 2 levels "A","P": 1 1 2 2 2 2 2 2 2 2 ... .. ..$ RACOON : Factor w/ 2 levels "A","P": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ LATITUDE : num [1:88] 50.9 50.9 50.9 50.9 50.9 ... .. ..$ LONGITUDE : num [1:88] 6.77 6.77 6.77 6.78 6.78 ... .. ..$ size : num [1:88, 1] -0.551 -0.551 -0.551 -0.551 -0.551 ... .. .. ..- attr(*, "scaled:center")= num 219823 .. .. ..- attr(*, "scaled:scale")= num 398721 .. ..$ years : num [1:88, 1] -0.572 -0.954 -0.858 -0.189 -0.189 ... .. .. ..- attr(*, "scaled:center")= num 11 .. .. ..- attr(*, "scaled:scale")= num 10.5 .. ..$ open.water : num [1:88, 1] -2.83 0.705 -1.527 -1.155 -1.155 ... .. .. ..- attr(*, "scaled:center")= num 81.1 .. .. ..- attr(*, "scaled:scale")= num 26.9 .. ..$ woods : num [1:88, 1] -0.745 -0.593 -0.782 -0.745 -0.593 ... .. .. ..- attr(*, "scaled:center")= num 20.6 .. .. ..- attr(*, "scaled:scale")= num 26.4 .. ..$ emerse.veg : num [1:88, 1] 2.996 -0.702 1.634 1.244 1.244 ... .. .. ..- attr(*, "scaled:center")= num 18 .. .. ..- attr(*, "scaled:scale")= num 25.7 .. ..$ reed : num [1:88, 1] -0.72 -0.72 2.41 1.79 2.25 ... .. .. ..- attr(*, "scaled:center")= num 23 .. .. ..- attr(*, "scaled:scale")= num 31.9 .. ..$ rushes : num [1:88, 1] 3.199 -0.858 -0.858 -1.072 -0.858 ... .. .. ..- attr(*, "scaled:center")= num 25.1 .. .. ..- attr(*, "scaled:scale")= num 23.4 .. ..$ avg.veg.ht : num [1:88, 1] -0.5254 -0.7554 0.5097 -0.2953 -0.0653 ... .. .. ..- attr(*, "scaled:center")= num 106 .. .. ..- attr(*, "scaled:scale")= num 87 .. ..$ total.area : num [1:88, 1] 0.167 0.167 0.167 0.167 0.167 ... .. .. ..- attr(*, "scaled:center")= num 150 .. .. ..- attr(*, "scaled:scale")= num 85.7 .. ..$ ac.area : num [1:88, 1] -0.0895 -0.0895 -0.0895 -0.0895 -0.0895 ... .. .. ..- attr(*, "scaled:center")= num 42.4 .. .. ..- attr(*, "scaled:scale")= num 16.9 .. ..$ inac.area : num [1:88, 1] 0.0895 0.0895 0.0895 0.0895 0.0895 ... .. .. ..- attr(*, "scaled:center")= num 57.6 .. .. ..- attr(*, "scaled:scale")= num 16.9 .. ..$ long_temp : num [1:88, 1] 1.09 1.09 1.09 1.09 1.09 ... .. .. ..- attr(*, "scaled:center")= num 10.1 .. .. ..- attr(*, "scaled:scale")= num 0.597 .. ..$ long_prec : num [1:88, 1] 0.485 0.485 0.485 0.485 0.485 ... .. .. ..- attr(*, "scaled:center")= num 721 .. .. ..- attr(*, "scaled:scale")= num 94.8 .. ..$ water.depth: num [1:88, 1] 0.485 0.485 0.485 0.485 0.485 ... .. .. ..- attr(*, "scaled:center")= num 721 .. .. ..- attr(*, "scaled:scale")= num 94.8 ..@ mapInfo : NULL ..@ obsToY : num [1:5, 1:5] 1 0 0 0 0 0 1 0 0 0 ...

Please let me know if you need any other data as well.

Thank you very much and best wishes!
Tanvi

Marc J. Mazerolle

unread,
Mar 13, 2024, 1:17:55 PM3/13/24
to unma...@googlegroups.com
Hi Tanvi,

The newdata object in modavgPred( ) must be a data.frame where columns
correspond to the different explanatory variables of occupancy
appearing in the different models. In the code you sent, you copied the
unmarkedFrameOccu( ) object. Thus, you should provide a data frame with
the different columns coding for the explanatory variables on occupancy
with the same name and type as they appear in the unmarkedFrameOccu( ).
Another piece of advice is to run modavgPred( ) separately for each
variable of interest. For example, if you want to show the variation in
occupancy vs water depth, it is easier to create a data frame where
water depth varies across values of interest and hold all the other
explanatory variables constant. You can then get predictions and plot
them easily. Repeat the same steps for another variable as needed.

Best,

Marc

--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et
aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

-------- Message initial --------
De: Tanvi Gurjar <tannu4...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: Re: [unmarked] how to modavgpred with obsCov
Date: Wed, 13 Mar 2024 09:27:31 -0700
> str(newdata)Formal class 'unmarkedFrameOccu' [package "unmarked"]
And, here is the original dataframe:> str(BM3)Formal class

Tanvi Gurjar

unread,
Mar 13, 2024, 3:53:17 PM3/13/24
to unmarked
Dear Marc, 

Thank you! It worked!
I have another associated question:
I have a total of 5 sampling occasions and 88 waterbodies for which I am modeling amphibian occupancy in mining quarries, and selected the following three models for L. vulgaris  (detection constant). 
1. Waterbody Age + water column depth + Fish pr/ab .............. c-hat = 2.05   P = 0.06; AIC = 213.09  
2. Waterbody Age + waterbody size + Fish pr/ab ...................... c-hat = 2.14  P = 0.04;  AIC = 212.12  
3. Waterbody Age + Fish pr/ab .................................................... c-hat = 1.99  P = 0.04;   AIC = 211.28

These did not look good to me therefore I simulated a small data set and checked the MacKenzie-Bailey GOF statistic for the 3 models with variation in the experimental design:
1. 100 sites, 8 sampling occassions (J)
2. 100 sites, 5 sampling occassions (J)
3. 88 sites, 8 sampling occassions (J)
Models fit on the simulated data returned c-hat values just under 1 (between 0.89 and 0.9), and p-values between 0.8-0.9.
In addition to the simulation, I also checked AIC and QAIC values and occupancy predictions by accounting for c-hat. 

My question is:
1.) If the simulated data return decent c-hat and p values for a specific set of covariates, is it best to stop there, 
2.) Or, is it better practive to select a model based on the original data by accounting for overdispersion (model selection based on QAIC for example)?

Sorry for the length of this query. 

I'd be grateful for any feedback on this forum :)!

Many thanks and best wishes,
Tanvi

Marc J. Mazerolle

unread,
Mar 14, 2024, 8:39:00 AM3/14/24
to unma...@googlegroups.com
Hi Tanvi,

If you simulate data, then it is normal not to have overdispersion when
you assess the fit of the simulated data, because you are drawing
observations from the correct distribution. Simulation is informative
if you want to compare the fit of the model on the observed data to
what you obtain when the model truly fits the data (simulated data). In
your case, you should be using the estimate of c-hat you obtained from
the goodness of fit test to adjust your inferences.

Marc
--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et
aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

-------- Message initial --------
De: Tanvi Gurjar <tannu4...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: Re: [unmarked] how to modavgpred with obsCov
Date: Wed, 13 Mar 2024 12:53:17 -0700
> 6.77 6.77 6.77 6.78 6.78 ... .. ..$ size : num [1:88, 1] -0.551 -
> 0.551
> 0.0895

Tanvi Gurjar

unread,
Mar 14, 2024, 9:28:34 AM3/14/24
to unmarked
Dear Marc, 

Thanks alot for your speedy response!
As a final question, when I try to use aictab to get QAIC values, I always get AIC values only. I don't know what I am doing wrong. The output is the same for both of these.  
Tab <- aictab(cand.set = (Cand.list), modnames = NULL, sort = FALSE) #<-- for AIC values
Model selection based on AIC: K AIC Delta_AIC AICWt Cum.Wt LL Mod1 6 409.02 0.00 0.98 0.98 -198.51 Mod2 5 417.22 8.20 0.02 1.00 -203.61 Mod3 2 426.77 17.75 0.00 1.00 -211.38

TabC <- aictab(cand.set = (Cand.list), modnames = NULL, c_hat=2.97, sort = TRUE) #for QAIC values
Model selection based on AIC:
K AIC Delta_AIC AICWt Cum.Wt LL Mod1 6 409.02 0.00 0.98 0.98 -198.51 Mod2 5 417.22 8.20 0.02 1.00 -203.61 Mod3 2 426.77 17.75 0.00 1.00 -211.38

Cand.list <- list(TM.m28,TM.m33,TM.m0)
Many thanks for your time!
Tanvi

Marc J. Mazerolle

unread,
Mar 14, 2024, 10:02:04 AM3/14/24
to unma...@googlegroups.com
Hi Tanvi,

Be careful to use the correct arguments. See ?aictab and other
functions. To adjust inferences for overdispersion, the correct
argument is c.hat = yourValue.

Marc
--
____________________________________
Marc J. Mazerolle
Professeur agrégé et directeur du bac. en environnements naturels et
aménagés
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 407120
Email: marc.ma...@sbf.ulaval.ca

-------- Message initial --------
De: Tanvi Gurjar <tannu4...@gmail.com>
Répondre à: unma...@googlegroups.com
À: unmarked <unma...@googlegroups.com>
Objet: Re: [unmarked] how to modavgpred with obsCov
Date: Thu, 14 Mar 2024 06:28:34 -0700
Reply all
Reply to author
Forward
0 new messages