Title: MSH surveys over posterior; Phylo lm

Name: Tammy L. Elliott

Date: August 29, 2017

R version 3.1

Losses related to threat categories

Data

#Show binary data.drame
head(conservation.losses)

##        Losses Threatened
## AbiBal      0          0
## AcaRho      0          0
## AcePen      0          0
## AceRub      0          0
## AceSac      0          0
## AceSpi      0          0

#Phylogenetic regression with phylolm 
species.losses.phylolm<-phylolm(conservation.losses$Losses~conservation.losses$Threatened, data=conservation.losses, phy=MSH.trees[[1]], model=("BM"), boot=100)

# iterate phylogenetic regression over 100 trees
species.comm.phylolm.100<-NULL
for(a in 1:length(MSH.trees)){
    species.comm.phylolm.100[[a]]<-
phylolm(conservation.losses$Losses~conservation.losses$Threatened, data=conservation.losses, phy=MSH.trees[[a]], model=("BM"), boot=100)
    }

# Extract p.value coefficients
species.comm.p.value<-NULL
for(b in 1:length(species.comm.phylolm.100)){
    species.comm.p.value[[b]]<-
        coef(summary(species.comm.phylolm.100[[b]]))[2,6]
}

#get mean
(species.comm.p.mean<-mean(species.comm.p.value))

## [1] 0.2544424

#Get CI's
species.comm.p.value.error <- qt(0.975, df=length(species.comm.p.value))*sd(species.comm.p.value)/sqrt(length(species.comm.p.value))
(species.comm.p.value.left<-species.comm.p.mean-species.comm.p.value.error)

## [1] 0.1888866

(species.comm.p.value.right<-species.comm.p.mean+species.comm.p.value.error)

## [1] 0.3199982

#Show min and max p-values; these range from 0 to 0.97 which seems very odd.
min(species.comm.p.value)

## [1] 0

max(species.comm.p.value)

## [1] 0.9761349

Phylolm with Native/non-native as predictor

native.nonnative.comm<-Quebec_comm_split$Introduced+Quebec_comm_split$Excluded+Quebec_comm_split$Ephemeral

#Create binary dataframe for losses and introductions
native.comm<-as.data.frame(cbind(losses.binary[,1], native.nonnative.comm))
colnames(native.comm)<-c("Gains", "Introductions")

#Show binary data.frame
head(native.comm)

##        Gains Introductions
## AbiBal     1             0
## AcaRho     1             0
## AcePen     1             0
## AceRub     1             0
## AceSac     1             0
## AceSpi     1             0

# iterate phylogenetic regression over 100 trees
native.comm.phylolm.100<-NULL
for(c in 1:length(MSH.trees)){
    native.comm.phylolm.100[[c]]<-
phylolm(native.comm$Gains~native.comm$Introductions, data=native.comm, phy=MSH.trees[[c]], model=("BM"), boot=100)
    }

# Extract p.value coefficients
native.comm.p.value<-NULL
for(d in 1:length(native.comm.phylolm.100)){
    native.comm.p.value[[d]]<-
        coef(summary(native.comm.phylolm.100[[d]]))[2,6]
}

native.comm.p.value

##   [1] 3.718691e-02 4.756015e-01 2.351690e-02 1.087240e-01 1.713091e-01
##   [6] 2.112648e-01 1.012663e-01 4.001877e-02 6.691027e-02 7.128125e-01
##  [11] 1.779987e-01 1.758079e-01 3.463491e-01 3.507254e-01 3.663887e-02
##  [16] 1.616716e-01 5.699342e-01 2.264445e-01 3.126525e-01 1.855065e-01
##  [21] 1.011161e-01 5.172563e-02 1.615243e-01 1.889579e-01 7.408735e-04
##  [26] 3.447815e-01 2.679576e-01 5.754870e-02 1.933266e-01 2.935309e-01
##  [31] 4.430349e-01 7.442255e-01 4.770920e-01 8.763057e-01 4.994008e-02
##  [36] 1.357261e-02 2.697909e-01 8.523317e-02 2.605850e-01 4.891217e-01
##  [41] 1.414688e-01 9.548334e-01 5.796724e-01 9.851396e-02 1.037120e-02
##  [46] 8.319001e-01 4.592048e-01 5.439470e-01 2.900548e-01 2.110543e-01
##  [51] 3.630504e-01 3.681948e-01 6.500416e-02 3.854283e-01 1.747735e-01
##  [56] 2.449159e-01 4.811574e-01 1.219782e-01 1.460818e-01 1.881263e-01
##  [61] 1.191935e-01 5.605826e-02 1.545792e-01 1.549567e-01 1.231021e-01
##  [66] 3.037826e-01 1.618706e-01 1.028833e-01 4.381171e-02 5.774541e-01
##  [71] 5.008421e-01 7.313069e-01 9.536656e-06 4.423367e-01 4.442908e-01
##  [76] 2.836564e-01 2.054338e-01 1.231735e-01 2.956902e-01 4.226726e-01
##  [81] 3.408382e-02 6.592621e-02 3.652300e-01 5.159239e-01 4.173356e-01
##  [86] 2.657095e-08 3.368214e-01 6.394464e-01 9.874880e-02 8.121002e-02
##  [91] 3.319448e-01 1.665097e-01 2.026667e-01 1.202082e-01 1.151768e-01
##  [96] 7.546019e-01 9.489337e-03 4.232076e-02 2.771078e-01 1.812639e-01

#get mean
(native.comm.p.mean<-mean(native.comm.p.value))

## [1] 0.265253

#Get CI's
native.comm.p.value.error <- qt(0.975, df=length(native.comm.p.value))*sd(native.comm.p.value)/sqrt(length(native.comm.p.value))
(native.comm.p.value.left<-native.comm.p.mean-native.comm.p.value.error)

## [1] 0.2219172

(native.comm.p.value.right<-native.comm.p.mean+native.comm.p.value.error)

## [1] 0.3085889