PAMPwrProject.Rmd

---
title: "Power Simulation Notes"
author: "Eiren Jacobson"
date: "April 18, 2017"
output: pdf_document
---

header-includes:
   - \usepackage{caption}
output: 
  pdf_document:
    latex_engine: xelatex

---
\begin{center}
\Large{Eiren Jacobson} 
\end{center}

\begin{center}
\large{`r format(Sys.time(), '%d %B %Y')`}
\end{center}

```{r ChunkOpts, echo=FALSE}
knitr::opts_chunk$set(echo=TRUE, results="hide", message=FALSE, warning=FALSE)
```

# Variability in the 3-year passive acoustic dataset

* The mean PPS values measured across all stations and years is shown in the following table

```{r Setup, results="hide", echo=FALSE}
 
 require(data.table)
 require(ggplot2)
 require(reshape2)
 require(nlme)
 require(plyr)
 
 load("./Data/CpodPPS.RData") # PPS.Data created by cpodParse
 load("./Data/CXY.RData") # CXY created by spDensity
 load("./Data/predGrid.RData") #gridOut created by spDensity

```


```{r, results="markup"}

# Need to change names to match 
names(CXY)[4]<- "DENSITY"
names(CXY)[1] <- "MOORING"

# Get annual mean PPS values for each station
total <-  PPS.Data[,.(Total = sum(PPS)), 
                           by = .(YEAR, MOORING)]

 vartab <-  total[,.(CV = sd(Total)/mean(Total)), 
                            by = .(MOORING)]
 m <- merge(vartab, CXY[,c(1,4)], by="MOORING")
 plot(m$DENSITY, m$CV)


 
 mod <- lm(CV ~ DENSITY, data=m)
# Attach the acoustic data to the aerial survey density estimates
total <- merge(CXY[,c(1,4)], total, by="MOORING")

all.pps <- merge(CXY[,c(1,4)], PPS.Data, by="MOORING")
 
# Create a spatial model of PPS as a function of year and underlying mean density
D.model <- glm(log(Total) ~ factor(YEAR, ordered=TRUE) + log(DENSITY), 
                data = total)


Dm <- glm(log(Total) ~ 0 + factor(YEAR, ordered=TRUE) + log(DENSITY), 
          family=quasi(link="identity", variance="mu"), data=total)

summary(Dm)
plot(Dm)

# glm(Total ~ factor(YEAR) + DENSITY, family=poisson(link="log"), data=total) AIC = 465658

#

# save(D.model, file="./Data/Dmodel.RData")



# Use model to predict expected PPS values as a function of density over the three years

pred.grid <- expand.grid("DENSITY"=seq(0.001, 3, by=0.001), "YEAR"=as.factor(2013:2015))

p.vals <- predict(Dm, pred.grid, se.fit=TRUE)

pred.grid <- cbind.data.frame(pred.grid, PredPPS=p.vals)

ggplot()+
  geom_line(data=pred.grid, 
            aes(x=DENSITY, y=exp(PredPPS.fit), color=YEAR))+
  geom_point(data=total, 
             aes(x=DENSITY, y=Total, color=as.factor(YEAR)))+
  ylab("Total (91 day) PPS")+
  theme_bw()



```

```{r}




```

```{r}
source("./Scripts/simSpatial")
n.moorings <- c(10, 25, 50, 75, 100)
p.change <- seq(-0.5, 0.5, by = 0.05)
nrep <- 1000
tab.res <- data.frame()

for (no in n.moorings){
  
  for (pc in p.change){
    
      r.est <- vector()
  
      for (i in 1:nrep){
    
        nd <- as.matrix(simSpatial(no, 10, pc, b = 0))

        nd <- as.data.frame(nd, names=c("MOORING", "YEAR", "DENSITY", "PPS"))
    
        nm <- lme(fixed = log(PPS) ~ YEAR, 
             random = ~ 1 | MOORING, data = nd)
        
        ye <- ifelse(summary(nm)$tTable[10]<0.05 & sign(summary(nm)$tTable[2])==sign(pc), 1, 0)
    
        r.est <- c(r.est, ye)
    
      } # end i in n.rep
      
      tab.res <- rbind.data.frame(tab.res, 
                                  cbind.data.frame("No.M" = no, "P.Change" = pc, 
                                                   "P.Dect" = sum(r.est)/length(r.est)))
      
  } # end pc in p.change
    
    
  } # end nm in no. moorings




cbPalette <- c("#E69F00", "#56B4E9", "#009E73","#0072B2",  "#CC79A7")
ggsave("./Figures/ChangeVPower.pdf",
ggplot()+
  geom_line(data=data.frame("P.Change"=c(-0.5, 0.5), "P.Dect"=c(0.8, 0.8)), 
            aes(x=P.Change, y=P.Dect), linetype=1, size=1)+
  geom_line(data=tab.res, aes(x=P.Change, y=P.Dect, 
                              group=as.factor(No.M), color=as.factor(No.M)), size=1)+
  scale_color_manual(name="No. Moorings", values=cbPalette)+
  
 # geom_line(data=data.frame("P.Change"=c(-0.5, -0.5), "P.Dect"=c(0.14, 0.33)), 
 #          aes(x=P.Change, y=P.Dect) ,size =3, color = "red")+
  xlab("Percent Change in Population")+
  ylab("Power")+
  scale_x_continuous(labels=scales::percent)+
  theme_bw()+
  theme(legend.key=element_blank()),
height=4, width=5, units="in")

table <- dcast(tab.res, No.M ~ P.Change)

print(table)


```


```{r}
require(reshape2)
load("./Data/simSpatialb0Resultsx10000.RData")
table <- dcast(tab.res, No.M ~ P.Change, value.var = "P.Dect")


```