length_based_index/model.R at master · DTUAqua/length_based_index · GitHub

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
## Default inputs
SPECIES  <- "Gadus morhua"
QUARTER  <- 4
KM       <- 50
MINSIZE  <- 4
MAXSIZE  <- 120
MINYEAR  <- 2000
MAXYEAR  <- 2015
BY       <- 2
DATFILE  <- "EBcod.RData"
OUTFILE  <- paste0("results", QUARTER, ".RData")

## For scripting
input <- parse(text=Sys.getenv("SCRIPT_INPUT"))
print(input)
eval(input)

## Load data
library(DATRAS)
d <- local({
    load(DATFILE)
    stopifnot(length(ls()) == 1)
    get(ls())
})
stopifnot( class(d) == "DATRASraw" )

## Make grid
library(gridConstruct)
grid <- gridConstruct(d,km=KM)
## plot(grid)
## map("worldHires",add=TRUE)

## Data subset
d <- addSpectrum(d,cm.breaks=seq(MINSIZE,MAXSIZE,by=BY))
d$haulid <- d$haul.id
d <- subset(d, Quarter == QUARTER, Gear != "GRT")
d <- subset(d, Year %in% MINYEAR:MAXYEAR )
d <- subset(d, 25<HaulDur & HaulDur<35 )
d <- as.data.frame(d)

library(mapdata)

## Set up time factor (careful with empty factor levels ! )
years <- as.numeric(levels(d$Year))
d$time <- factor(d$Year, levels=min(years):max(years))

## Set up spatial factor and grid
## grid <- gridConstruct(d,nearestObs=100)
d$position <- gridFactor(d,grid)
Q0 <- -attr(grid,"pattern")
diag(Q0) <- 0; diag(Q0) <- -rowSums(Q0)
I <- .symDiagonal(nrow(Q0))

## Set up design matrix
## TODO: Gear by sizeGroup
##A <- sparse.model.matrix( ~ sizeGroup:time + Gear - 1, data=d)
A0 <- sparse.model.matrix( ~ sizeGroup:time - 1, data=d)
A <- sparse.model.matrix( ~ Gear - 1, data=d)
A <- A[,-which.max(table(d$Gear)),drop=FALSE]

B <- cbind2(A,A0); B <- t(B)%*%B
if(min(eigen(B)$val)<1e-8)stop("Singular B")

data <- as.list(d)
data$A <- A
data$I <- I
data$Q0 <- Q0

data <- data[!sapply(data,is.character)]
data <- data[!sapply(data,is.logical)]

library(TMB)
compile("model.cpp")
dyn.load(dynlib("model"))
obj <- MakeADFun(
    data=data,
    parameters=list(
        logdelta= 0 ,
        logkappa= 0 ,
        tphi_time= 0 ,
        tphi_size = 0 ,
        logsigma= 0 ,
        beta= rep(0, ncol(A)) ,
        eta= array(0, c(nrow(Q0), nlevels(d$time), nlevels(d$sizeGroup) ) ),
        etanug= array(0, c(nlevels(d$sizeGroup), nlevels(d$haulid) ) ),
        etamean = array(0, c(nlevels(d$sizeGroup) , nlevels(d$time)) )
        ),
    DLL="model",
    random=c("eta","etanug","etamean","beta")
    )

print(obj$par)
runSymbolicAnalysis(obj)
system.time(obj$fn())

system.time(opt <- nlminb(obj$par, obj$fn, obj$gr))

system.time(sdr <- sdreport(obj))

system.time(sdr2 <- sdreport(obj,bias.correct=TRUE,hessian.fixed=solve(sdr$cov.fixed)))
mat <- summary(sdr2,"report")
rownames(mat) <- NULL
df <- as.data.frame(mat)
df$unbiased <- sdr2$unbiased$value
df <- cbind(df,
            expand.grid(sizeGroup=levels(d$sizeGroup),time=levels(d$time))
            )
xtabs(Estimate ~ sizeGroup + time,data=df)
x1 <- xtabs(N ~ sizeGroup + time,data=d) / xtabs( ~ sizeGroup + time,data=d)
x2 <- xtabs(Estimate ~ sizeGroup + time,data=df)
x3 <- xtabs(unbiased ~ sizeGroup + time,data=df)

save(sdr, df, file=OUTFILE)

## pdf("plotQ4.pdf")
## matplot(x1,type="l",main="Raw average")
## matplot(x2,type="l",main="Posterior mode")
## matplot(x3,type="l",main="Posterior mean")
## dev.off()