Occupancy in RMAR: how to extract specific results?

[Triciaserow]

Hi,

I'm currently working on occupancy models using RMARK. I have a huge number (about 500) of input files for occupancy modeling. For example:

>OV_C91_1=convert.inp("H:\\My Documents\\Analysis\\inp\\OV_C91_1.txt")
>OV_C91_1.p.dot=mark(OV_C91_1, model="Occupancy")
Output summary for Occupancy model
Name : p(~1)Psi(~1)

Npar : 2
-2lnL: 59.24371
AICc : 65.64371

Beta
estimate se lcl ucl
p:(Intercept) -1.597419 0.4406167 -2.461028 -0.7338101
Psi:(Intercept) 1.059176 1.1936992 -1.280475 3.3988264


Real Parameter p
1 2 3 4 5 6 7 8 9 10
0.1683427 0.1683427 0.1683427 0.1683427 0.1683427 0.1683427 0.1683427 0.1683427 0.1683427 0.1683427

Real Parameter Psi
1
0.742533


The question is, how do I select the outputs that I want (p, psi, se of p, se of psi, AICc...) from all my models and summarize them into a table?

The column of the table will be the name of the parameters, and the rows are all the models (OV_C91_1, OV_C91_2, OV_C91_3....................).

Thank you.


Tricia
[/code]

[jlaake]

It depends on what you want betas or reals? If it is betas use

summary(model)$beta

if reals

summary(model)$real$Psi or $p

You can use rbind,cbind as you need to combine the results.

--jeff

[Triciaserow]

It depends on what you want betas or reals? If it is betas use

summary(model)$beta

if reals

summary(model)$real$Psi or $p

You can use rbind,cbind as you need to combine the results.

--jeff

Thank you for your reply. I will try those commands.

--Tricia

[Triciaserow]

Hi,

Another question is, how can I pull out the standard errors of psi and p? From the summary(model) in RMark, it only gave me the se of Beta. I have to type the name of the model, which gave me a new window of notepad and that's where I find the se of p and psi (from the real function of parameters).

For example:

>summary(A.p.dot)

Beta
estimate se lcl ucl
p:(Intercept) -1.350850 0.3323397 -2.002236 -0.6994641
Psi:(Intercept) 3.552349 5.5165369 -7.260063 14.3647620


>A.p.dot
LOGIT Link Function Parameters of { p(~1)Psi(~1) }
95% Confidence Interval
Parameter Beta Standard Error Lower Upper
------------------------- -------------- -------------- -------------- --------------
1:p:(Intercept) -1.3508499 0.3323397 -2.0022358 -0.6994641
2:Psi:(Intercept) 3.5523492 5.5165369 -7.2600634 14.364762


Real Function Parameters of { p(~1)Psi(~1) }
95% Confidence Interval
Parameter Estimate Standard Error Lower Upper
------------------------- -------------- -------------- -------------- --------------
1:p g1 a0 t1 0.2057315 0.0543063 0.1189684 0.3319311
2:Psi g1 a0 t1 0.9721411 0.1494031 0.7025695E-003 0.9999994


I need the real se rather than se of betas. I have thousands of data sets to run and if RMark can pull that value out, it will save me a lot of time.

Thank you.

Tricia

[jlaake]

You need to spend some time working through the Appendix and the workshop material and learning some more R. But in the meantime,here is an example that uses the weta example for occupancy model. Only the code below the #### row is relevant to your question. Everything above that is from the weta example.

data(weta)
# Create function to fit the 18 models in the book
fit.weta.models=function()
{
# use make.time.factor to create time-varying dummy variables Obs1 and Obs2
# observer 3 is used as the intercept
weta=make.time.factor(weta,"Obs",1:5,intercept=3)
# Process data and use Browse covariate to group sites; it could have also
# been used an individual covariate because it is a 0/1 variable.
weta.process=process.data(weta,model="Occupancy",groups="Browse")
weta.ddl=make.design.data(weta.process)
# time factor variable copied to Day to match names used in book
weta.ddl$p$Day=weta.ddl$p$time
# Define p models
p.dot=list(formula=~1)
p.day=list(formula=~Day)
p.obs=list(formula=~Obs1+Obs2)
p.browse=list(formula=~Browse)
p.day.obs=list(formula=~Day+Obs1+Obs2)
p.day.browse=list(formula=~Day+Browse)
p.obs.browse=list(formula=~Obs1+Obs2+Browse)
p.day.obs.browse=list(formula=~Day+Obs1+Obs2+Browse)
# Define Psi models
Psi.dot=list(formula=~1)
Psi.browse=list(formula=~Browse)
# Create model list
cml=create.model.list("Occupancy")
# Run and return marklist of models
return(mark.wrapper(cml,data=weta.process,ddl=weta.ddl))
}
weta.models=fit.weta.models()
##################################################################################
# create dataframe of real estimates for all weta models
estimates=sapply(weta.models[1:nrow(weta.models$model.table)],function(x)
c(summary(x,se=T)$real$p$estimate,summary(x,se=T)$real$Psi$estimate))
estimates=t(estimates)
colnames(estimates)=c(row.names(summary(weta.models[[1]],se=T)$real$p),row.names(summary(weta.models[[1]],se=T)$real$Psi))
# create dataframe of std error of real estimates for all weta models
se.estimates=sapply(weta.models[1:nrow(weta.models$model.table)],function(x)
c(summary(x,se=T)$real$p$se,summary(x,se=T)$real$Psi$se))
se.estimates=t(se.estimates)
colnames(se.estimates)=c(row.names(summary(weta.models[[1]],se=T)$real$p),row.names(summary(weta.models[[1]],se=T)$real$Psi))