-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path03_Dredge_hake_Larv.R
109 lines (75 loc) · 3.53 KB
/
03_Dredge_hake_Larv.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
rm(list=ls()) # Clear variables in workspace
library(MuMIn) # load package with dredge
##### prepare fish data ####
# main = getwd()
MainFile = "/Users/cdvestfals/Desktop/Hake_Analysis_for_GitHub"
DredgeFile = paste0(MainFile,"/03_Dredge")
dir.create(DredgeFile)
setwd(DredgeFile)
####################################################################
##### BRING IN HAKE DATA - FIRST-FEEDING LARVAL STAGE ####################
####################################################################
fish00 = data.frame(read.table(paste0(MainFile,'/00_Data/hake.data_2020_medians.csv'), header=TRUE, sep=","))
head(fish00)
max(fish00$year)
#years = 1980:max(fish00$year)
years = 1980:2010 # only have ROMS data 1980-2010
fish = fish00[fish00$year %in% years,]
roms = data.frame(read.csv(paste0(MainFile,'/01_ROMS_Prelim/Data_ROMS.for.hake.analysis.larv.csv'), header=TRUE))
roms = roms[roms$year %in% years,]
##########################################
#### RUN MODEL FITTING ###################
##########################################
# roms = the roms data after processing
#### Build equation ####
for(j in 2:ncol(roms)){# build equation
x1 = colnames(roms)[j]
if(j==ncol(roms)){x2 = x1}else{x2 = paste(x1,"+")}
if(j==2){ROMS = x2}else{ROMS = paste(ROMS,x2)}
}
ROMS
#### Set up Formula #### variables with AIC < 2 (AIC - AICq)
form_dredge = as.formula(paste('dev ~ I(CSTlarv.s^2) +', ROMS, sep=''))
#
form_dredge
saveRDS(form_dredge, "Formula_for_Dredge.hake_larv.rds")
#### Prep & combine data ####
# data_1 = # combine roms data and resids
data_1 = merge(fish,roms, 'year')
data_1 = na.omit(data_1)
write.table(data_1, 'Analyzed.data.hake_larv.csv',sep=',',col.names = TRUE, row.names = FALSE)
###############################################################################
##### run dredge twice ########################################################
##### once full, once minus the last five years ###############################
##### save output ###############################
###############################################################################
fit = lm(form_dredge, data=data_1, na.action = na.fail)
#### Run Dredge and save results ####
mtable = dredge(fit, rank = AICc, m.lim = c(NA,5),
subset= # block correlated terms from the model > 0.70
!(UWlarv.bn && UWlarv.bs) && # 0.90
!(UWlarv.bn && UWlarv.cn) && # 0.85
!(UWlarv.bs && UWlarv.cn) && # 0.79
!(UWlarv.bs && UWlarv.cs) && # 0.81
!(UWlarv.cn && UWlarv.cs) && # 0.77
!(DDlarv && SSH.ja.c) && # 0.82
!(DDlarv && SSH.ja.s) && # 0.80
!(DDlarv && EKE.ja.s) && # 0.71
!(SSH.ja.c && SSH.ja.s) && # 0.92
!(SSH.ja.c && UWlarv.bn) && # -0.79
!(SSH.ja.s && UWlarv.bn) && # -0.71
dc(CSTlarv.s, I(CSTlarv.s^2)),
extra = list(R2 = function(x)
summary(x)$r.squared[[1]],
F = function(y)
summary(y)$fstatistic[[1]]),
trace=2 )
mtable2 = subset(mtable, delta < 2)
mtable2
print(mtable2, abbrev.names = FALSE)
saveRDS(mtable, file='Table_model.fits.hake_larv.rds')
write.csv(mtable2,"R_Table_Delta2.hake_larv.csv")
mtable = readRDS('Table_model.fits.hake_larv.rds')
mtable2 = subset(mtable, delta < 2)
mtable2
setwd(MainFile)