www.phidot.org

Hi,

I noticed that when extracting results from a POPAN model (modelName$results$derived) that only the estimates for B_i were given and not N_i. I wrote a function to extract the estimated N_hat table from the markxxx.out file and thought I would pass it along in case others would find this useful.

Code: Select all

getNhat <- function(wd,outName, nocc) {

  prev.wd <- getwd()

  setwd(wd)

  mark.out <- readLines(outName)

  header.index <- grep("\f", mark.out)
  header.index <- c(header.index-1,header.index, header.index+1,header.index+2,header.index+3)

  mark.out <- mark.out[-header.index]

  Nhat.index <- grep("N-hat", mark.out) + 1
  Nhat.string.tab <- mark.out[(Nhat.index):(Nhat.index + nocc)]

  writeLines(Nhat.string.tab,"temp.txt")
  Nhat.tab <- read.table("temp.txt", header=TRUE)
  unlink("temp.txt")
  names(Nhat.tab) <- c("grp","occ","nhat","se","lcl", "ucl")
 
  setwd(prev.wd)

  return(Nhat.tab)

}


It's certainly not the most robust function, but it's done the job for so far.

Code: Select all

> getNhat(wd="C:\\Documents and Settings\\ .  . . . . . river\\",
+         outName="mark001.out",
+         nocc=mod.1$nocc)
   grp occ        nhat         se        lcl         ucl
1    1   1    2.120904   3.715364  -5.161211    9.403018
2    1   2  956.449100 303.977800 360.652600 1552.245600
3    1   3  725.686630 123.914280 482.814630  968.558620

. . . . .


I'm also not sure how similar or different this is comapred to what is done within RMark.

Thanks for your contribution. I was aware of that problem and have not worked out why I'm not retrieving those values from the mark binary file. That was the reason that I wrote popan.derived in the package. Note that there are many different values of N that can be estimated for this type of model.

--jeff

Ha, never mind then. It looks like popan.derived is the way to go ! I just didn't know about it.