extracting p-value from occu output

[Shannon]

I'll preface this by saying I'm primarily a SAS person (over R), but
I've been able to do this with "standard" analyses (e.g., lm or glm).

Basically, I want to pull a single number from the output, in this
case the p-value (yes, I know, "significance testing" is passe, but
some people still want it and at this point it's mostly curiosity...
and that nagging sense of "this just HAS to be possible...").

In an "ordinary" analysis, I can get things like that with the [$]
operator. Not so in S4 classes. I have figured out that [@] works
like [$] in S4 classes, but only sort of.

For the sake of illustration:
####################################################################
> fm <- occu(~HabType ~HabType, data=umf)
> summary(fm)

Call:
occu(formula = ~HabType ~ HabType, data = umf)

Occupancy (logit-scale):
Estimate SE z P(>|z|)
(Intercept) 0.245 0.248 0.988 0.323
HabType -0.489 0.351 -1.392 0.164

Detection (logit-scale):
Estimate SE z P(>|z|)
(Intercept) 2.30 0.334 6.90 5.28e-12
HabType 0.49 0.570 0.86 3.90e-01

AIC: 294.4622
Number of sites: 132
optim convergence code: 0
optim iterations: 18
Bootstrap iterations: 0
####################################################################

I want to extract the p-value for "HabType" in occupancy (0.164) and
store it as a variable.

I have tried using str(fm) and through that managed to get down to:
####################################################################
fmest <- fm@estimates
fmestst <- fmest@estimates$state
> fmestst@estimates
(Intercept) HabType
0.2453471 -0.4885544
####################################################################
fm@estimates@estimates$state@estimates did not work.


But I don't even see a "slot" for p-values
####################################################################
> str(fmestst)
Formal class 'unmarkedEstimate' [package "unmarked"] with 7 slots
..@ name : chr "Occupancy"
..@ short.name : chr "psi"
..@ estimates : Named num [1:2] 0.245 -0.489
.. ..- attr(*, "names")= chr [1:2] "(Intercept)" "HabType"
..@ covMat : num [1:2, 1:2] 0.0616 -0.0616 -0.0616 0.1232
..@ covMatBS : NULL
..@ invlink : chr "logistic"
..@ invlinkGrad: chr "logistic.grad"
####################################################################

Wondering if this was intentional (i.e., "we don't want you using p-
values") or if it's more about my ignorance of S4 classes.

Thanks for your time!
Shannon

[Andy Royle]

 hi Shannon, if you look at str(summary(fm)) you will see a component for p-values -- see this:
regards
andy




> str(summary(fm))

Call:
occu(formula = ~obsvar1 ~ 1, data = pferUMF)

Occupancy (logit-scale):
 Estimate   SE     z P(>|z|)

     8.16 21.7 0.377   0.707

Detection (logit-scale):
            Estimate    SE       z  P(>|z|)

(Intercept)   -1.922 0.164 -11.745 7.51e-32
obsvar1       -0.113 0.169  -0.671 5.02e-01

AIC: 262.8795
Number of sites: 130
optim convergence code: 0
optim iterations: 86
Bootstrap iterations: 0

List of 2
 $ state:'data.frame':  1 obs. of  4 variables:
  ..$ Estimate: num 8.16
  ..$ SE      : num 21.7
  ..$ z       : num 0.377
  ..$ P(>|z|) : num 0.707
 $ det  :'data.frame':  2 obs. of  4 variables:
  ..$ Estimate: num [1:2] -1.922 -0.113
  ..$ SE      : num [1:2] 0.164 0.169
  ..$ z       : num [1:2] -11.745 -0.671
  ..$ P(>|z|) : num [1:2] 7.51e-32 5.02e-01
>

From:	Shannon <shannon...@gmail.com>
To:	unmarked <unma...@googlegroups.com>
Date:	02/14/2012 12:04 PM
Subject:	[unmarked] extracting p-value from occu output
Sent by:	unma...@googlegroups.com

---

[Pablo García Díaz]

Dear users,

I just want to share a concern about the algoritm pcountOpen: is it any improvement in the speed of the models planned? Running pcountOpen with covariates is somewhat exasperating. This afternoon-evening I was trying to run the following model:

columnas<-c("2006", "2007", "2008", "2009", "2010", "2011")

filas<-c("Camaces", "Saldeana", "Uces", "Huebra", "Monleras", "Saucelle", "Aldeadavila", "Agueda", "Trabanca", "Villagonzalo", 

"Riolobos", "Salamanca","FA", "Juzbado", "Pelayos", "Santa", "Cristo", "Castraz", "CR", "Béjar", "Francia")

y<-matrix(c(1,1,1,1,1,1,

2,2,2,1,2,2,

1,1,1,1,1,1,

2,2,2,1,1,2,

1,1,0,1,0,1,

1,2,0,1,2,1,

2,2,2,1,2,2,

2,2,2,2,2,2,

2,2,1,2,1,1,

2,1,0,2,1,1,

1,1,0,1,2,2,

2,2,1,2,2,2,

4,4,3,3,3,3,

3,3,2,3,2,2,

1,0,1,0,0,1,

1,1,1,2,1,1,

1,0,1,0,1,2,

2,2,2,2,2,2,

2,2,2,3,3,3,

1,1,2,2,2,1,

2,2,2,2,2,2), nrow=21,ncol=6,dimnames=list(filas, columnas), byrow=TRUE)

library(unmarked)

prey<-data.frame(prey=c(7,24,19,36,7,5,10,23,4,12,6,12,27,16,10,5,20,31,11,4,7))

umf<- unmarkedFramePCO(y=y, siteCovs=prey,numPrimary=6)

fit2<-pcountOpen(~1, ~1, ~prey, ~1, umf, mixture ="P", K=150)

summary(fit2)

I was doing the work in and independent powerful computer, only doing such work. After 4 hours I was required to quite R and turn off the computer. The model was not computed yet. 

Of course, this is just a comment.

Pablo

[Richard Chandler]

Hi Pablo,

Thanks for your comment. This is what the help file says:

Warning

This function can be extremely slow, especially if there are covariates of gamma or omega. Consider testing the timing on a small subset of the data, perhaps with se=FALSE. Finding the lowest value of K that does not affect estimates will also help with speed.


So I lowered K to 30 and it only took 54s to fit the model. Then I did a quick check, and it seems as though K=30 is sufficient for these data.

system.time({
fit2<-pcountOpen(~1, ~1, ~prey, ~1, umf, mixture ="P", K=30,
                 control=list(trace=TRUE, REPORT=1))
})
(fit2check <- update(fit2, K=40, se=FALSE))


Richard



_____________________________________
Richard Chandler, post-doc
USGS Patuxent Wildlife Research Center
301-497-5696

From:	Pablo García Díaz <her...@hotmail.com>
To:	<unma...@googlegroups.com>
Date:	02/14/2012 01:41 PM
Subject:	[unmarked] comment on pcountOpen
Sent by:	unma...@googlegroups.com

---

[michael...@mail.huji.ac.il]

Dear users,

I stumbled upon this post while trying to extract p-values from a model fitted using function distsamp. Unlike with occu, this seems to be impossible. For example:

data(linetran)

ltUMF <- with(linetran, {

   unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4),

   siteCovs = data.frame(Length, area, habitat),

   dist.breaks = c(0, 5, 10, 15, 20),

   tlength = linetran$Length * 1000, survey = "line", unitsIn = "m")

   })

ltUMF

summary(ltUMF)

hist(ltUMF)

(fm1 <- distsamp(~ 1 ~ 1, ltUMF))

Now, it seems as though summary returns NULL when operating on models fitted with distsamp. Executing x=summary(fm1) just print the summary on screen, while x is NULL.

Is there any other way to extract the p-values from a distsamp model?

Thanks,

Michael

[Richard Chandler]

Hi Michael,

This is pretty clunky, but it should work:

summary(fm1@estimates)$state

--
You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Richard Chandler

Assistant Professor

Warnell School of Forestry and Natural Resources

University of Georgia

706-542-5815

Lab website

[michael...@mail.huji.ac.il]

Dear Richard,

This works, thank you very much!

Michael