Skip to content

Commit 92fa801

Browse files
committed
fixing bug in mixed and overall simplifying censoring of distances
1 parent 25c334d commit 92fa801

File tree

1 file changed

+39
-47
lines changed

1 file changed

+39
-47
lines changed

R/parameterize.R

Lines changed: 39 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -71,12 +71,11 @@ si_lnorm2 <- function(ttree, params, cutoff = NULL) {
7171
#' @export
7272
dist_gamma_mixed <- function(ttree, params, cutoff = NULL) {
7373
if(is.null(cutoff)) {
74-
ttree[, dist_prob := fcase(owned & dist_diff > 100,
75-
dgamma(dist_diff, shape= params$DK_shape, scale = params$DK_scale),
76-
!owned & dist_diff > 100,
77-
dgamma(dist_diff, shape= params$DK2_shape, scale = params$DK2_scale),
78-
dist_diff <= 100,
79-
pgamma(100, shape = params$DK_shape, scale = params$DK_scale)/100)]
74+
ttree[, dist_diff_c := ifelse(dist_diff < 100, 100, dist_diff)]
75+
ttree[, dist_prob := fifelse(owned,
76+
dgamma(dist_diff_c, shape= params$DK_shape, scale = params$DK_scale),
77+
dgamma(dist_diff_c, shape= params$DK2_shape, scale = params$DK2_scale))]
78+
ttree[, dist_diff_c := NULL]
8079
} else {
8180
# return the cutoff value given a prob (either length 1 or length of the ttree)
8281
ifelse(ttree$owned, qgamma(cutoff, shape = params$DK_shape, scale = params$DK_scale),
@@ -91,15 +90,14 @@ dist_gamma_mixed <- function(ttree, params, cutoff = NULL) {
9190
#' @export
9291
dist_weibull_mixed <- function(ttree, params, cutoff = NULL) {
9392
if(is.null(cutoff)) {
94-
ttree[, dist_prob := fcase(owned & dist_diff > 100,
95-
dweibull(dist_diff, shape= params$DK_shape_weibull,
96-
scale = params$DK_scale_weibull),
97-
!owned & dist_diff > 100,
98-
dweibull(dist_diff, shape= params$DK2_shape_weibull,
99-
scale = params$DK2_scale_weibull),
100-
dist_diff <= 100,
101-
pweibull(100, shape = params$DK_shape_weibull,
102-
scale = params$DK_scale_weibull)/100)]
93+
94+
ttree[, dist_diff_c := fifelse(dist_diff < 100, 100, dist_diff)]
95+
ttree[, dist_prob := fifelse(owned,
96+
dweibull(dist_diff_c, shape= params$DK_shape_weibull,
97+
scale = params$DK_scale_weibull),
98+
dweibull(dist_diff_c, shape= params$DK2_shape_weibull,
99+
scale = params$DK2_scale_weibull))]
100+
ttree[, dist_diff_c := NULL]
103101
} else {
104102
# return the cutoff value given a prob (either length 1 or length of the ttree)
105103
ifelse(ttree$owned, qweibull(cutoff, shape = params$DK_shape_weibull,
@@ -115,12 +113,12 @@ dist_weibull_mixed <- function(ttree, params, cutoff = NULL) {
115113
#' @export
116114
dist_lnorm_mixed <- function(ttree, params, cutoff = NULL) {
117115
if(is.null(cutoff)) {
118-
ttree[, dist_prob := fcase(owned & dist_diff > 100,
119-
dlnorm(dist_diff, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog),
120-
!owned & dist_diff > 100,
121-
dlnorm(dist_diff, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog),
122-
dist_diff <= 100,
123-
plnorm(100, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog)/100)]
116+
117+
ttree[, dist_diff_c := fifelse(dist_diff < 100, 100, dist_diff)] # censor
118+
ttree[, dist_prob := fifelse(owned,
119+
dlnorm(dist_diff_c, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog),
120+
dlnorm(dist_diff_c, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog))]
121+
ttree[, dist_diff_c := NULL] # ditching censored distance difference
124122
} else {
125123
# return the cutoff value given a prob (either length 1 or length of the ttree)
126124
ifelse(ttree$owned, qlnorm(cutoff, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog),
@@ -135,10 +133,9 @@ dist_lnorm_mixed <- function(ttree, params, cutoff = NULL) {
135133
#' @export
136134
dist_gamma1<- function(ttree, params, cutoff = NULL) {
137135
if(is.null(cutoff)) {
138-
ttree[, dist_prob := fcase(dist_diff > 100,
139-
dgamma(dist_diff, shape= params$DK_shape, scale = params$DK_scale),
140-
dist_diff <= 100,
141-
pgamma(100, shape = params$DK_shape, scale = params$DK_scale)/100)]
136+
ttree[, dist_prob := fifelse(dist_diff > 100,
137+
dgamma(dist_diff, shape= params$DK_shape, scale = params$DK_scale),
138+
dgamma(100, shape = params$DK_shape, scale = params$DK_scale))]
142139
} else {
143140
# return the cutoff value given a prob (either length 1 or length of the ttree)
144141
qgamma(cutoff, shape = params$DK_shape, scale = params$DK_scale)
@@ -152,12 +149,11 @@ dist_gamma1<- function(ttree, params, cutoff = NULL) {
152149
#' @export
153150
dist_weibull1<- function(ttree, params, cutoff = NULL) {
154151
if(is.null(cutoff)) {
155-
ttree[, dist_prob := fcase(dist_diff > 100,
156-
dweibull(dist_diff, shape= params$DK_shape_weibull,
152+
ttree[, dist_prob := fifelse(dist_diff > 100,
153+
dweibull(dist_diff, shape= params$DK_shape_weibull,
157154
scale = params$DK_scale_weibull),
158-
dist_diff <= 100,
159-
pweibull(100, shape = params$DK_shape_weibull,
160-
scale = params$DK_scale_weibull)/100)]
155+
dweibull(100, shape = params$DK_shape_weibull,
156+
scale = params$DK_scale_weibull))]
161157
} else {
162158
# return the cutoff value given a prob (either length 1 or length of the ttree)
163159
qweibull(cutoff, shape = params$DK_shape_weibull, scale = params$DK_scale_weibull)
@@ -171,10 +167,9 @@ dist_weibull1<- function(ttree, params, cutoff = NULL) {
171167
#' @export
172168
dist_lnorm1 <- function(ttree, params, cutoff = NULL) {
173169
if(is.null(cutoff)) {
174-
ttree[, dist_prob := fcase(dist_diff > 100,
175-
dlnorm(dist_diff, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog),
176-
dist_diff <= 100,
177-
plnorm(100, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog)/100)]
170+
ttree[, dist_prob := fifelse(dist_diff > 100,
171+
dlnorm(dist_diff, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog),
172+
dlnorm(100, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog))]
178173
} else {
179174
# return the cutoff value given a prob (either length 1 or length of the ttree)
180175
qlnorm(cutoff, meanlog = params$DK_meanlog, sdlog = params$DK_sdlog)
@@ -188,10 +183,9 @@ dist_lnorm1 <- function(ttree, params, cutoff = NULL) {
188183
#' @export
189184
dist_gamma2 <- function(ttree, params, cutoff = NULL) {
190185
if(is.null(cutoff)) {
191-
ttree[, dist_prob := fcase(dist_diff > 100,
192-
dgamma(dist_diff, shape= params$DK2_shape, scale = params$DK2_scale),
193-
dist_diff <= 100,
194-
pgamma(100, shape = params$DK2_shape, scale = params$DK2_scale)/100)]
186+
ttree[, dist_prob := fifelse(dist_diff > 100,
187+
dgamma(dist_diff, shape= params$DK2_shape, scale = params$DK2_scale),
188+
dgamma(100, shape = params$DK2_shape, scale = params$DK2_scale))]
195189
} else {
196190
# return the cutoff value given a prob (either length 1 or length of the ttree)
197191
qgamma(cutoff, shape = params$DK2_shape, scale = params$DK2_scale)
@@ -205,12 +199,11 @@ dist_gamma2 <- function(ttree, params, cutoff = NULL) {
205199
#' @export
206200
dist_weibull2 <- function(ttree, params, cutoff = NULL) {
207201
if(is.null(cutoff)) {
208-
ttree[, dist_prob := fcase(dist_diff > 100,
209-
dweibull(dist_diff, shape= params$DK2_shape_weibull,
202+
ttree[, dist_prob := fifelse(dist_diff > 100,
203+
dweibull(dist_diff, shape= params$DK2_shape_weibull,
210204
scale = params$DK2_scale_weibull),
211-
dist_diff <= 100,
212-
pweibull(100, shape = params$DK2_shape_weibull,
213-
scale = params$DK2_scale_weibull)/100)]
205+
dweibull(100, shape = params$DK2_shape_weibull,
206+
scale = params$DK2_scale_weibull))]
214207
} else {
215208
# return the cutoff value given a prob (either length 1 or length of the ttree)
216209
qweibull(cutoff, shape = params$DK2_shape_weibull, scale = params$DK2_scale_weibull)
@@ -224,10 +217,9 @@ dist_weibull2 <- function(ttree, params, cutoff = NULL) {
224217
#' @export
225218
dist_lnorm2 <- function(ttree, params, cutoff = NULL) {
226219
if(is.null(cutoff)) {
227-
ttree[, dist_prob := fcase(dist_diff > 100,
228-
dlnorm(dist_diff, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog),
229-
dist_diff <= 100,
230-
plnorm(100, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog)/100)]
220+
ttree[, dist_prob := fifelse(dist_diff > 100,
221+
dlnorm(dist_diff, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog),
222+
dlnorm(100, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog))]
231223
} else {
232224
# return the cutoff value given a prob (either length 1 or length of the ttree)
233225
qlnorm(cutoff, meanlog = params$DK2_meanlog, sdlog = params$DK2_sdlog)

0 commit comments

Comments
 (0)