rm(list=ls())


#install unmarked if not installed
#install.packages("unmarked")
#load library
library(unmarked)
## Loading required package: reshape
## Loading required package: lattice
## Loading required package: parallel
## Loading required package: Rcpp
bfly<-read.csv("../Butterfly.csv")
bfly<-subset(bfly,select=c(Field,Border,Transect,Treatment,Date,Visit,CLLESSU))
library(reshape)
bfly$Transect=as.factor(bfly$Transect)

junk<-melt(bfly,id.var=c("Transect","Visit"),measure.var="CLLESSU")

j2=cast(junk,Transect ~ Visit)

covs<-subset(bfly,select=c(Transect,Treatment))
covs<-as.data.frame(with(covs,table(Transect,Treatment)))
covs<-subset(covs,Freq>0)

site.cov<-merge(j2,covs,by="Transect")[,8]
site.cov<-data.frame(treatment=site.cov)


y<-j2[,2:7]

y<-cbind(j2$"1",j2$"2",j2$"3",j2$"4",j2$"5",j2$"6")

S <- nrow(y) # number of sites
J <- 2 # number of secondary sampling occasions
T <- 3 # number of primary period
#dummy variables to createseasonal effects (need 2 for 3 occasions)
season<-as.factor(rep(c(1,2,3),S))

yearly.cov<-data.frame(season)
#additive time effects within season. to model effects change between seasons create an interaction term
secondary<-as.factor(rep(c(1,2,1,2,1,2),S))

second.cov<-data.frame(secondary)
butterfly<-unmarkedMultFrame(y=y,siteCovs=site.cov, numPrimary=T,yearlySiteCovs=yearly.cov,obsCovs=second.cov)
fm1<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, 
    pformula = ~1, data=butterfly) 
#back transformations
backTransform(fm1,'det')
## Backtransformed linear combination(s) of Detection estimate(s)
## 
##  Estimate    SE LinComb (Intercept)
##     0.497 0.033 -0.0116           1
## 
## Transformation: logistic
backTransform(fm1,"psi")
## Backtransformed linear combination(s) of Initial estimate(s)
## 
##  Estimate     SE LinComb (Intercept)
##     0.771 0.0607    1.22           1
## 
## Transformation: logistic
backTransform(fm1,"col")
## Backtransformed linear combination(s) of Colonization estimate(s)
## 
##  Estimate   SE LinComb (Intercept)
##     0.437 0.11  -0.254           1
## 
## Transformation: logistic
backTransform(fm1,"ext")
## Backtransformed linear combination(s) of Extinction estimate(s)
## 
##  Estimate     SE LinComb (Intercept)
##    0.0909 0.0529    -2.3           1
## 
## Transformation: logistic
fm2<-colext(psiformula = ~treatment, gammaformula = ~1, epsilonformula = ~1, 
    pformula = ~1, data=butterfly) 
fm3<-colext(psiformula= ~1, gammaformula = ~treatment, epsilonformula = ~1, 
    pformula = ~1, data=butterfly) 

fm4<-colext(psiformula= ~1, gammaformula = ~1, epsilonformula = ~treatment, 
    pformula = ~1, data=butterfly) 

fm1
## 
## Call:
## colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, 
##     pformula = ~1, data = butterfly)
## 
## Initial:
##  Estimate    SE    z  P(>|z|)
##      1.22 0.344 3.54 0.000406
## 
## Colonization:
##  Estimate    SE      z P(>|z|)
##    -0.254 0.447 -0.568    0.57
## 
## Extinction:
##  Estimate   SE    z P(>|z|)
##      -2.3 0.64 -3.6 0.00032
## 
## Detection:
##  Estimate    SE       z P(>|z|)
##   -0.0116 0.132 -0.0879    0.93
## 
## AIC: 1024.421
fm2
## 
## Call:
## colext(psiformula = ~treatment, gammaformula = ~1, epsilonformula = ~1, 
##     pformula = ~1, data = butterfly)
## 
## Initial:
##                       Estimate     SE     z P(>|z|)
## (Intercept)              0.243  0.619 0.392   0.695
## treatmentburncontrol     0.581  0.828 0.702   0.482
## treatmentdisk            7.775 18.476 0.421   0.674
## treatmentdiskcontrol     1.263  0.856 1.476   0.140
## treatmentfieldcontrol    0.298  0.817 0.365   0.715
## 
## Colonization:
##  Estimate    SE      z P(>|z|)
##   -0.0662 0.403 -0.164   0.869
## 
## Extinction:
##  Estimate   SE     z  P(>|z|)
##     -2.07 0.57 -3.63 0.000279
## 
## Detection:
##  Estimate   SE     z P(>|z|)
##    0.0528 0.13 0.406   0.685
## 
## AIC: 1022.907
fm3
## 
## Call:
## colext(psiformula = ~1, gammaformula = ~treatment, epsilonformula = ~1, 
##     pformula = ~1, data = butterfly)
## 
## Initial:
##  Estimate   SE    z  P(>|z|)
##      1.19 0.33 3.59 0.000329
## 
## Colonization:
##                       Estimate   SE      z P(>|z|)
## (Intercept)              -2.94 3.48 -0.844   0.399
## treatmentburncontrol      3.07 3.53  0.868   0.385
## treatmentdisk             2.95 3.94  0.749   0.454
## treatmentdiskcontrol      2.04 3.48  0.585   0.558
## treatmentfieldcontrol     3.86 3.62  1.068   0.285
## 
## Extinction:
##  Estimate    SE     z  P(>|z|)
##     -2.42 0.687 -3.52 0.000437
## 
## Detection:
##  Estimate    SE      z P(>|z|)
##    -0.022 0.131 -0.168   0.867
## 
## AIC: 1025.79
fm4
## 
## Call:
## colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~treatment, 
##     pformula = ~1, data = butterfly)
## 
## Initial:
##  Estimate    SE    z  P(>|z|)
##      1.19 0.341 3.49 0.000487
## 
## Colonization:
##  Estimate    SE      z P(>|z|)
##     -0.18 0.432 -0.417   0.677
## 
## Extinction:
##                       Estimate    SE      z P(>|z|)
## (Intercept)             -1.112 0.958 -1.161   0.246
## treatmentburncontrol    -1.803 1.717 -1.050   0.294
## treatmentdisk           -0.589 1.280 -0.460   0.645
## treatmentdiskcontrol    -2.822 3.511 -0.804   0.422
## treatmentfieldcontrol   -0.539 1.176 -0.459   0.646
## 
## Detection:
##  Estimate    SE      z P(>|z|)
##   0.00403 0.132 0.0304   0.976
## 
## AIC: 1029.132