Title: MSH surveys over posterior; Phylo lm - phylogenetic binomial logistic regression

Name: Tammy L. Elliott

Date: August 29, 2017

R version 3.1

#see model description here:http://127.0.0.1:11389/library/phylolm/html/phyloglm.html
#note fitting of a logistic regression because our response is binary
#Iterate over 100 trees
model.gains<-NULL
for(a in 1:length(model.trees)){
    model.gains[[a]]<-
        phyloglm(gains~natives, data=addition.data, phy=model.trees[[a]], method=("logistic_MPLE"), boot=0)
}

summary(model.gains[[1]])

## 
## Call:
## phyloglm(formula = gains ~ natives, data = addition.data, phy = model.trees[[a]], 
##     method = ("logistic_MPLE"), boot = 0)
##        AIC     logLik Pen.logLik 
##      776.0     -385.0     -381.1 
## 
## Method: logistic_MPLE
## Mean tip height: 413.3026
## Parameter estimate(s):
## alpha: 0.1184042 
## 
## Coefficients:
##             Estimate   StdErr z.value p.value    
## (Intercept)  0.28975  0.16456  1.7608 0.07828 .  
## natives     -1.73879  0.20243 -8.5897 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Note: Wald-type p-values for coefficients, conditional on alpha=0.1184042

Phylogenetic binomial regression results for gains due to natives

(model.gains.natives.estimate.mean<-mean(model.gains.natives.estimate))

## [1] -1.788931

#Get CI's
model.gains.natives.estimate.error <- qnorm(0.975)*sd(model.gains.natives.estimate)/sqrt(length(model.gains.natives.estimate))
(model.gains.natives.estimate.left<-model.gains.natives.estimate.mean-model.gains.natives.estimate.error)

## [1] -1.808908

(model.gains.natives.estimate.right<-model.gains.natives.estimate.mean+model.gains.natives.estimate.error)

## [1] -1.768955

#natives standard error
(model.gains.natives.stderr.mean<-mean(model.gains.natives.stderr))

## [1] 0.2036202

#Get CI's
model.gains.natives.stderr.error <- qnorm(0.975)*sd(model.gains.natives.stderr)/sqrt(length(model.gains.natives.stderr))
(model.gains.natives.stderr.left<-model.gains.natives.stderr.mean-model.gains.natives.stderr.error)

## [1] 0.2031398

(model.gains.natives.stderr.right<-model.gains.natives.stderr.mean+model.gains.natives.stderr.error)

## [1] 0.2041006

#natives z values
(model.gains.natives.z.mean<-mean(model.gains.natives.z))

## [1] -8.781413

#Get CI's
model.gains.natives.z.error <- qnorm(0.975)*sd(model.gains.natives.z)/sqrt(length(model.gains.natives.z))
(model.gains.natives.z.left<-model.gains.natives.z.mean-model.gains.natives.z.error)

## [1] -8.858446

(model.gains.natives.z.right<-model.gains.natives.z.mean+model.gains.natives.z.error)

## [1] -8.704381

#get mean of p.value for natives
(model.gains.p.mean<-mean(model.gains.p.value))

## [1] 1.444937e-17

#Get CI's
model.gains.p.value.error <- qnorm(0.975)*sd(model.gains.p.value)/sqrt(length(model.gains.p.value))
(model.gains.p.value.left<-model.gains.p.mean-model.gains.p.value.error)

## [1] 7.625819e-18

(model.gains.p.value.right<-model.gains.p.mean+model.gains.p.value.error)

## [1] 2.127293e-17

Phylogenetic binomial regression results for losses with special conservation staus

#note fitting of a logistic regression because our response is binary
#Iterate over 100 trees
model.losses<-NULL
for(f in 1:length(losses.trees)){
    model.losses[[f]]<-
        phyloglm(losses~conservation, data=not.found.data, phy=losses.trees[[f]],method=("logistic_MPLE"), boot=0)
    }

summary(model.losses[[1]])

## 
## Call:
## phyloglm(formula = losses ~ conservation, data = not.found.data, 
##     phy = losses.trees[[f]], method = ("logistic_MPLE"), boot = 0)
##        AIC     logLik Pen.logLik 
##      471.4     -232.7     -229.9 
## 
## Method: logistic_MPLE
## Mean tip height: 413.3026
## Parameter estimate(s):
## alpha: 0.05590981 
## 
## Coefficients:
##              Estimate   StdErr  z.value   p.value    
## (Intercept)  -2.26480  0.16808 -13.4747 < 2.2e-16 ***
## conservation  1.75761  0.38673   4.5448 5.499e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Note: Wald-type p-values for coefficients, conditional on alpha=0.05590981

#conservation estimate
(model.losses.conservation.estimate.mean<-mean(model.losses.conservation.estimate))

## [1] 1.845282

#Get CI's
model.losses.conservation.estimate.error <- qnorm(0.975)*sd(model.losses.conservation.estimate)/sqrt(length(model.losses.conservation.estimate))
(model.losses.conservation.estimate.left<-model.losses.conservation.estimate.mean-model.losses.conservation.estimate.error)

## [1] 1.818594

(model.losses.conservation.estimate.right<-model.losses.conservation.estimate.mean+model.losses.conservation.estimate.error)

## [1] 1.87197

#conservation standard error
(model.losses.conservation.stderr.mean<-mean(model.losses.conservation.stderr))

## [1] 0.3855088

#Get CI's
model.losses.conservation.stderr.error <- qnorm(0.975)*sd(model.losses.conservation.stderr)/sqrt(length(model.losses.conservation.stderr))
(model.losses.conservation.stderr.left<-model.losses.conservation.stderr.mean-model.losses.conservation.stderr.error)

## [1] 0.3850706

(model.losses.conservation.stderr.right<-model.losses.conservation.stderr.mean+model.losses.conservation.stderr.error)

## [1] 0.385947

#conservation z values
(model.losses.conservation.z.mean<-mean(model.losses.conservation.z))

## [1] 4.787925

#Get CI's
model.losses.conservation.z.error <- qnorm(0.975)*sd(model.losses.conservation.z)/sqrt(length(model.losses.conservation.z))
(model.losses.conservation.z.left<-model.losses.conservation.z.mean-model.losses.conservation.z.error)

## [1] 4.715818

(model.losses.conservation.z.right<-model.losses.conservation.z.mean+model.losses.conservation.z.error)

## [1] 4.860031

#get mean of p.value of conservation
(model.losses.p.mean<-mean(model.losses.p.value))

## [1] 8.181746e-06

#Get CI's
model.losses.p.value.error <- qnorm(0.975)*sd(model.losses.p.value)/sqrt(length(model.losses.p.value))
(model.losses.p.value.left<-model.losses.p.mean-model.losses.p.value.error)

## [1] 3.780588e-06

(model.losses.p.value.right<-model.losses.p.mean+model.losses.p.value.error)

## [1] 1.258291e-05