-
Notifications
You must be signed in to change notification settings - Fork 1
/
fill_funcs.R
189 lines (175 loc) · 6.62 KB
/
fill_funcs.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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
#' Combine prior and data for transition matrix
#'
#' \code{fill_transition} returns the expected value of the transition
#' matrix combining observed transitions for one time step and a prior
#'
#' @param TF A list of two matrices, T and F, as ouput by \code{\link[popbio]{projection.matrix}}.
#' @param N A vector of observed transitions.
#' @param P A matrix of the priors for each column. Defaults to uniform.
#' @param priorweight total weight for each column of prior as a percentage of sample size or 1 if negative
#' @param returnType A character vector describing the desired return value. Defaults to "T" the transition matrix.
#'
#' @return The return value depends on parameter returnType.
#' \itemize{
#' \item A - the summed matrix "filled in" using a dirichlet prior
#' \item T - just the filled in transition matrix
#' \item TN - the augmented matrix of fates -- use in calculating the CI or for simulation
#' }
#'
#' @export
#'
fill_transitions <- function(TF, N, P = NULL, priorweight = -1, returnType = "T") {
check_TF(TF)
Tmat <- TF$T
Fmat <- TF$F
order <- dim(Tmat)[1]
if (missing(P)) {
# fill in with a uniform prior <- <- <-
P <- matrix(1 / (order + 1), nrow = order + 1, ncol = order)
} else {
if (ncol(P) != order | nrow(P) != (order + 1)) {
stop("Bad dimensions on P")
}
}
Tfilled <- matrix(NA, nrow = order, ncol = order)
TN <- matrix(NA, nrow = order + 1, ncol = order)
for (i in 1:order) {
observed <- Tmat[, i] * N[i]
if (priorweight > 0 & N[i] > 0) {
P[, i] <- P[, i] * priorweight * N[i]
}
allfates <- c(observed, N[i] - sum(observed)) + P[, i]
# missing <- allfates == 0
# allfates[missing] <- 1
# allfates[!missing] <- allfates[!missing] + sum(missing)
Tfilled[, i] <- allfates[1:order] / sum(allfates)
TN[, i] <- allfates
}
if (returnType == "A") {
return(Tfilled + Fmat)
} else if (returnType == "T") {
return(Tfilled)
} else if (returnType == "TN") {
return(TN)
} else {
stop("Bad returntype in fill_transitions()")
}
}
#' Combine prior and data for fertility matrix
#'
#' \code{fill_fertility} returns the expected value of the fertility
#' matrix combining observed recruits for one time step and a Gamma prior for each column.
#'
#' Assumes that only one stage reproduces ... needs generalizing.
#'
#' @param TF A list of two matrices, T and F, as ouput by \code{\link[popbio]{projection.matrix}}.
#' @param N A vector of observed stage distribution.
#' @param alpha A matrix of the prior parameter for each stage. Impossible stage combinations marked with NA_real_.
#' @param beta A matrix of the prior parameter for each stage. Impossible stage combinations marked with NA_real_.
#' @param priorweight total weight for each column of prior as a percentage of sample size or 1 if negative
#' @param returnType A character vector describing the desired return value. Defaults to "F" the fertility matrix
#'
#' @return The return value depends on parameter returnType.
#' \itemize{
#' \item A - the summed matrix "filled in" using a Gamma prior
#' \item F - just the filled in fertility matrix
#' \item ab - the posterior parameters alpha and beta as a list.
#' }
#'
#' @export
#'
fill_fertility <- function(TF, N, alpha = 0.00001, beta = 0.00001, priorweight = -1, returnType = "F") {
check_TF(TF)
Tmat <- TF$T
Fmat <- TF$F
order <- dim(Tmat)[1]
if (length(N) != order | sum(is.na(N)) > 0) {
stop("N isn't the correct length or has missing values.")
}
if (is.null(dim(alpha)) & length(alpha) != 1 | is.null(dim(beta)) & length(beta) != 1) {
stop("alpha or beta is not a matrix or a single value.")
}
if (!(is.numeric(alpha) & is.numeric(beta))) {
stop("alpha or beta must be numeric matrices or single values.")
}
if ((length(alpha) != order^2 | length(beta) != order^2)) {
warning("length(alpha | beta) != order^2: only using first value of alpha and beta")
alpha <- matrix(rep(alpha[1], order^2), nrow = order, ncol = order)
beta <- matrix(rep(beta[1], order^2), nrow = order, ncol = order)
}
Ffilled <- matrix(NA, nrow = order, ncol = order)
babies_next_year <- sweep(Fmat, 2, N, FUN = "*")
# test reproducing stages with beta, because of divide by zero issues
reproducing_stages <- apply(beta, 2, function(x) sum(!is.na(x))) > 0
# matrix multiplication doesn't preserve the column structure
if ((all(N[reproducing_stages] > 0) | sum(beta, na.rm = TRUE) > 0)) {
if (priorweight > 0) {
alpha_post <- sweep(alpha, 2, N, FUN = "*") * priorweight + babies_next_year
beta_post <- sweep(beta, 2, N, FUN = "*") * priorweight + N
} else {
alpha_post <- alpha + babies_next_year
beta_post <- sweep(beta, 2, N, FUN = "+")
}
Ffilled <- alpha_post / beta_post
Ffilled[is.na(Ffilled)] <- 0
} else {
stop("No reproducing stages in N, and no positive values in beta.")
}
if (returnType == "A") {
return(Tmat + Ffilled)
} else if (returnType == "F") {
return(Ffilled)
} else if (returnType == "ab") {
return(list(alpha = alpha_post, beta = beta_post))
} else {
stop("Bad returntype in fill_fertility()")
}
}
#' Helper functions for generating different priors
#'
#' Extract the number of individuals in each stage from a dataframe
#' of transitions.
#'
#' @param transitions a dataframe of observations of individuals in different stages
#' @param stage the name of the variable with the stage information
#' @param sort a vector of stage names in the desired order. Default is the order of levels in stage.
#'
#' @return a vector of the counts of observations in each level of stage.
#' @export
#'
get_state_vector <- function(transitions, stage = NULL,
sort = NULL) {
if (missing(stage)) {
stage <- "stage"
}
nl <- as.list(1:ncol(transitions))
names(nl) <- names(transitions)
stage <- eval(substitute(stage), nl, parent.frame())
if (is.null(transitions[, stage])) {
stop("No stage column matching ", stage)
}
if (missing(sort)) {
sort <- levels(transitions[[stage]])
}
tf <- table(transitions[[stage]])[sort]
return(as.vector(tf))
}
check_TF <- function(TF) {
if(typeof(TF) != "list") {
stop("TF must be a list of 2 matrices with equal dimensions")
}
if(length(TF) != 2) {
stop("TF must be a list of 2 matrices with equal dimensions")
}
if(is.null(TF$T) || is.null(TF$T)) {
stop("the matrices in TF must be named 'T' and 'F'")
}
Tmat <- TF$T
Fmat <- TF$F
if(!identical(dim(Tmat), dim(Fmat))) {
stop("The transition matrix's dimensions don't match the fertility matrix's dimensions")
}
if(dim(Tmat)[1] != dim(Tmat)[2]) {
stop("T and F must be square matrices")
}
}