-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathAMRO_InteractionModels_Combined.R
127 lines (112 loc) · 4.14 KB
/
AMRO_InteractionModels_Combined.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
####################################
#
# Nest success code
# Created by: Max Henschell
# Modified: 2 June 2015
#
####################################
#srm(list=ls()) #clean out the system
####################################
#
# Install required packages
#
####################################
packages.list<-c('MuMIn','BMA','xtable','RCurl', 'reshape', 'ggplot2', 'corrplot')
required.packages <- packages.list
new.packages <- required.packages[!(required.packages %in% installed.packages()[,'Package'])]
if(length(new.packages)) install.packages(new.packages, repos = 'http://cran.rstudio.com/')
# packages.list2<-c('coefplot2')
# required.packages2 <- packages.list2
# new.packages2 <- required.packages2[!(required.packages2 %in% installed.packages()[,'Package'])]
# if(length(new.packages2)) install.packages(new.packages2, ,
# repos='http://www.math.mcmaster.ca/bolker/R',
# type='source')
# packages.list<-c('BMA','xtable','RCurl', 'bbmle', 'reshape', 'ggplot2', 'BEST', 'arm', 'MCMCpack', 'popbio', 'extrafont', 'coefplot2')
lapply(packages.list, require, character.only=T)
#
# library(BMA)
# library(xtable)
# Logistical Exposure Link Function
# See Shaffer, T. 2004. A unifying approach to analyzing nest success.
# Auk 121(2): 526-540.
# library(MASS)
# logexp <- function(exposure = 1)
# {
# linkfun <- function(mu) qlogis(mu^(1/exposure))
# linkinv <- function(eta) plogis(eta)^exposure
# mu.eta <- function(eta) exposure * plogis(eta)^(exposure-1) *
# .Call(stats:::C_logit_mu_eta, eta, PACKAGE = 'stats')
# valideta <- function(eta) TRUE
# link <- paste('logexp(', deparse(substitute(exposure)), ')',
# sep='')
# structure(list(linkfun = linkfun, linkinv = linkinv,
# mu.eta = mu.eta, valideta = valideta,
# name = link),
# class = 'link-glm')
# }
#
# panel.cor <- function(x, y, digits = 2, cex.cor, ...)
# {
# usr <- par("usr"); on.exit(par(usr))
# par(usr = c(0, 1, 0, 1))
# # correlation coefficient
# r <- cor(x, y)
# txt <- format(c(r, 0.123456789), digits = digits)[1]
# txt <- paste("r= ", txt, sep = "")
# text(0.5, 0.6, txt)
#
# # p-value calculation
# p <- cor.test(x, y)$p.value
# txt2 <- format(c(p, 0.123456789), digits = digits)[1]
# txt2 <- paste("p= ", txt2, sep = "")
# if(p<0.01) txt2 <- paste("p= ", "<0.01", sep = "")
# text(0.5, 0.4, txt2)
# }
####################################
#
# Source logistic exposure and correlation plot code from github
#
####################################
source_https <- function(url, ...) {
# load package
require(RCurl)
# parse and evaluate each .R script
sapply(c(url, ...), function(u) {
eval(parse(text = getURL(u, followlocation = TRUE, cainfo = system.file('CurlSSL', 'cacert.pem', package = 'RCurl'))), envir = .GlobalEnv)
})
}
source_https('https://raw.github.com/maxhenschell/R-functions/master/logexp.R',
'https://raw.github.com/maxhenschell/R-functions/master/CorrPlots.R')
####################################
#
# AMRO
#
####################################
AMRO <- read.csv('DSR-AMRO.csv', header = T, na.strings = 'NA')
AMRO <- AMRO[AMRO$EXP != 0,]#Remove exp = 0 (first visit or visits after termination)
AMRO <- AMRO[complete.cases(AMRO),]
AMRO <- AMRO[AMRO$Study %in% c("MAH","MJP"),]
AMRO$Study <- factor(AMRO$Study)
# AMRO$Year <- factor(AMRO$Year)
AMRO$GPH[AMRO$GPH == 0.00] <- (min(AMRO$GPH[AMRO$GPH > 0]))/2
AMRO$OrdDate <- AMRO$OrdDate-min(AMRO$OrdDate)+1
AMRO$OrdDate2 <- AMRO$OrdDate^2
AMRO$Year <- as.factor(AMRO$Year)
AMRO1 <- AMRO
AMRO.nl <- AMRO1[!duplicated(AMRO[,1]),]
#Center and scale
for (i in 11:23){
AMRO[,i] <- (AMRO[,i] - mean(AMRO[,i]))/sd(AMRO[,i])
}
####################################
#
# Bayesian Model Averaging
# Develop a full model, then run code
#
####################################
AMRO.BMA <- bic.glm(FATE ~
CoreArea*Study + Edge60m*Study +
WBH*GPH*Study + TrailDist*GPH*Study, glm.family = binomial(logexp(exposure = AMRO$EXP)), data = AMRO)
imageplot.bma(AMRO.BMA)
summary(AMRO.BMA)
BMA.df <- data.frame(summary(AMRO.BMA))