utilityFunctions.R 11.4 KB
Newer Older
1 2 3 4
#' @include mortalityTable.R
NULL


5
fitExtrapolationLaw = function(data, ages, data.ages = ages, Dx = NULL, Ex = NULL, qx = NULL, method = "LF2", law = "HP", fit = 75:99, extrapolate = 80:120, fadeIn = 80:90, fadeOut = NULL, verbose = FALSE) {
6 7 8 9 10 11 12 13 14 15 16
    # Add the extrapolate ages to the needed ages
    neededAges = union(ages, extrapolate)
    # constrain the fit and fade-in range to given ages
    fit = intersect(ages, fit)
    if (!is.null(fadeIn))
        fadeIn = intersect(ages, fadeIn)
    if (!is.null(fadeOut))
        fadeOut = intersect(ages, fadeOut)

    # Hohe Alter: Fitte Heligman-Pollard im Bereich 75-99
    fitLaw = MortalityLaw(
17
        x = data.ages, Dx = Dx, Ex = Ex, qx = qx,
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
        law = law, opt.method = method,
        fit.this.x = fit)
    # summary(fitAP.m.75.99)
    # plot(fitAP.m.75.99)
    qPredict = predict(fitLaw, extrapolate)

    weights = rep(0, length(neededAges))
    names(weights) = neededAges

    if (!is.null(fadeIn)) {
        weights[neededAges < min(fadeIn)] = 0
        fadeInLen = length(fadeIn);
        weights[match(fadeIn, neededAges)] = (0:(fadeInLen - 1)) / (fadeInLen - 1)
        weights[neededAges > max(fadeIn)] = 1
    } else if (!is.null(fadeOut)) {
        weights[neededAges < min(fadeOut)] = 1
        fadeOutLen = length(fadeOut);
        weights[match(fadeOut, neededAges)] = ((fadeOutLen - 1):0) / (fadeOutLen - 1)
        weights[neededAges > max(fadeOut)] = 0
    }

    probs = fillAges(qPredict, givenAges = extrapolate, neededAges = neededAges, fill = 0) * weights +
        fillAges(data, givenAges = ages, neededAges = neededAges, fill = 0) * (1 - weights)

    if (verbose) {
        list(probs = probs, law = fitLaw, weights = weights)
    } else {
        probs
    }
}




# Fit an exponential function exp(-A*(x-x0)) to the last value (f(100) and f'(100) need to coincide):
fitExpExtrapolation = function(data, idx, up = TRUE, verbose = FALSE) {
    # browser()
    # Anchor point of the extrapolation
    f0 = data[[idx]]
    if (up) {
        A = -(data[[idx]] - data[[idx - 1]]) / f0
    } else {
        A = -(data[[idx + 1]] - data[[idx]]) / f0
    }
    x0 = idx + (log(f0) / A)
    fun.extra = function(x) exp(-A*(x - x0))
    if (up) {
        newdata = c(data[1:idx], sapply((idx + 1):length(data), fun.extra))
    } else {
        newdata = c(sapply(1:(idx - 1), fun.extra), data[idx:length(data)])
    }
    if (verbose) {
        list(data = newdata, A = A, x0 = x0, idx = idx)
    } else {
        newdata
    }
}


#' @export
78
mT.setName = function(table, name) {
79 80 81 82 83
    if (is.array(table)) {
        return(array(
            lapply(table, mT.setName, name = name),
            dim = dim(table), dimnames = dimnames(table))
        )
84 85
    } else if (is.list(table)) {
        return(lapply(table, mT.setName, name = name))
86
    }
87
    if (!is(table, "mortalityTable"))
88
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
89 90 91 92 93 94 95 96 97

    table@name = name
    table
}



#' @export
mT.fillAges = function(table, neededAges, fill = 0) {
98 99 100 101 102
    if (is.array(table)) {
        return(array(
            lapply(table, mT.fillAges, neededAges = neededAges, fill = fill),
            dim = dim(table), dimnames = dimnames(table))
        )
103 104
    } else if (is.list(table)) {
        return(lapply(table, mT.fillAges, neededAges = neededAges, fill = fill))
105
    }
106
    if (!is(table, "mortalityTable"))
107
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127

    existingAges = ages(table)
    if (.hasSlot(table, "ages"))
        table@ages = neededAges
    if (.hasSlot(table, "deathProbs"))
        table@deathProbs = fillAges(table@deathProbs, givenAges = existingAges, neededAges = neededAges, fill = fill)
    if (.hasSlot(table, "exposures") && !is.null(table@exposures) && length(table@exposures) > 1)
        table@exposures = fillAges(table@exposures, givenAges = existingAges, neededAges = neededAges, fill = 0)
    if (.hasSlot(table, "trend") && !is.null(table@trend) && length(table@trend) > 1)
        table@trend = fillAges(table@trend, givenAges = existingAges, neededAges = neededAges, fill = 0)
    if (.hasSlot(table, "trend2") && !is.null(table@trend2) && length(table@trend2) > 1)
        table@trend2 = fillAges(table@trend2, givenAges = existingAges, neededAges = neededAges, fill = 0)
    if (.hasSlot(table, "loading") && !is.null(table@loading) && length(table@loading) > 1)
        table@loading = fillAges(table@loading, givenAges = existingAges, neededAges = neededAges, fill = 0)
    if (!is.null(table@data$deaths))
        table@data$deaths = fillAges(table@data$deaths, givenAges = existingAges, neededAges = neededAges, fill = 0)
    table
}

#' @export
128
mT.scaleProbs = function(table, factor = 1.0, name.postfix = "scaled", name = NULL) {
129 130 131 132 133
    if (is.array(table)) {
        return(array(
            lapply(table, mT.scaleProbs, factor = factor, name.postfix = name.postfix, name = name),
            dim = dim(table), dimnames = dimnames(table))
        )
134 135
    } else if (is.list(table)) {
        return(lapply(table, mT.scaleProbs, factor = factor, name.postfix = name.postfix, name = name))
136 137 138 139
    }
    if (!is(table, "mortalityTable"))
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")

140
    table@deathProbs = factor * table@deathProbs
141 142 143 144 145
    if (is.null(name)) {
        if (!is.null(name.postfix)) {
            name = paste(table@name, name.postfix)
        }
    }
146 147 148 149 150 151 152 153
    if (!is.null(name)) {
        table@name = name
    }
    table
}


#' @export
154
mT.setTrend = function(table, trend, trendages = NULL, baseYear = NULL, dampingFunction = identity) {
155 156 157 158 159
    if (is.array(table)) {
        return(array(
            lapply(table, mT.setTrend, trend = trend, trendages = trendages, baseYear = baseYear, dampingFunction = dampingFunction),
            dim = dim(table), dimnames = dimnames(table))
        )
160 161
    } else if (is.list(table)) {
        return(lapply(table, mT.setTrend, trend = trend, trendages = trendages, baseYear = baseYear, dampingFunction = dampingFunction))
162
    }
163
    if (!is(table, "mortalityTable"))
164
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
165 166 167

    t = mortalityTable.trendProjection(
        table,
168 169
        baseYear = if (is.null(baseYear)) table@baseYear else baseYear,
        trend = trend[match(table@ages, if (is.null(trendages)) ages(table) else trendages)],
170
        dampingFunction = dampingFunction
171 172 173 174 175 176 177 178 179 180 181
    )
    t
}
#' @describeIn mT.setTrend Add a trend to the mortality table (returns a mortalityTable.trendProjection obect)
#' @export
mT.addTrend = mT.setTrend



#' @export
mT.extrapolateTrendExp = function(table, idx, up = TRUE) {
182 183 184 185 186
    if (is.array(table)) {
        return(array(
            lapply(table, mT.extrapolateTrendExp, idx = idx, up = up),
            dim = dim(table), dimnames = dimnames(table))
        )
187 188
    } else if (is.list(table)) {
        return(lapply(table, mT.extrapolateTrendExp, idx = idx, up = up))
189
    }
190
    if (!is(table, "mortalityTable"))
191
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
192 193 194

    if (.hasSlot(table, "trend") && !is.null(table@trend))
        table@trend = fitExpExtrapolation(table@trend, idx = idx,up = up)
195
    if (.hasSlot(table, "trend2") && !is.null(table@trend2))
196 197 198 199 200 201
        table@trend2 = fitExpExtrapolation(table@trend2, idx = idx,up = up)
    table
}


#' @export
202
mT.translate = function(table, baseYear, name = NULL) {
203 204 205 206 207
    if (is.array(table)) {
        return(array(
            lapply(table, mT.translate, baseYear = baseYear, name = name),
            dim = dim(table), dimnames = dimnames(table))
        )
208 209
    } else if (is.list(table)) {
        return(lapply(table, mT.translate, baseYear = baseYear, name = name))
210
    }
211
    if (!is(table, "mortalityTable"))
212
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
213 214 215

    table@deathProbs = periodDeathProbabilities(table, Period = baseYear)
    table@baseYear = baseYear
216 217 218
    if (!is.null(name)) {
        table@name = name
    }
219 220 221 222 223 224
    table
}


#' @export
mT.extrapolateProbsExp = function(table, age, up = TRUE) {
225 226 227 228 229
    if (is.array(table)) {
        return(array(
            lapply(table, mT.extrapolateProbsExp, age = age, up = up),
            dim = dim(table), dimnames = dimnames(table))
        )
230 231
    } else if (is.list(table)) {
        return(lapply(table, mT.extrapolateProbsExp, age = age, up = up))
232
    }
233
    if (!is(table, "mortalityTable"))
234
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251

    if (.hasSlot(table, "deathProbs")) {
        idx = match(age, ages(table))
        fit = fitExpExtrapolation(table@deathProbs, idx = idx, up = up, verbose = TRUE)
        table@deathProbs = fit$data
        table@data$extrapolationData = c(
            table@data$extrapolationData,
            list(list(law = "Exp", idx = idx, up = up, fit = fit)))
    }
    table
}


#' @export
mT.fitExtrapolationLaw = function(table, method = "LF2", law = "HP",
                                  fit = 75:99, extrapolate = 80:120,
                                  fadeIn = 80:90, fadeOut = NULL) {
252 253 254
    if (!is(table, "mortalityTable"))
        stop("First argument must be a mortalityTable.")

255
    ages = ages(table)
256 257 258
    # if (!is.null(table@exposures) && !is.na(table@exposures)) {
        # Ex = table@exposures
        # qx = table@deathProbs
259 260 261 262 263
        # if (!is.null(table@data$deaths)) {
        #     Dx = table@data$deaths
        # } else {
        #     Dx = table@deathProbs * Ex
        # }
264 265
    # } else {
        # Ex = rep(1, length(ages))
266
        # Dx = table@deathProbs
267 268
        # qx = table@deathProbs
    # }
269 270 271
    table  = mT.fillAges(table, neededAges = union(ages, extrapolate), fill = 0)
    fitted = fitExtrapolationLaw(
        data = table@deathProbs, ages = ages(table),
272
        qx = table@deathProbs, data.ages = ages,
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
        method = method, law = law,
        fit = fit, extrapolate = extrapolate,
        fadeIn = fadeIn, fadeOut = fadeOut,
        verbose = TRUE
    )
    # Store all fit parameters in the data slot of the mortality table
    table@data$extrapolationData = c(
        table@data$extrapolationData,
        list(list(law = law, method = method, fit = fit,
                  extrapolate = extrapolate, fadeIn = fadeIn, fadeOut = fadeOut,
                  fit = fitted)))
    table@deathProbs = fitted$probs

    table
}

289 290
#' @export
mT.setDimInfo = function(table, ..., append = TRUE) {
291 292 293 294 295
    if (is.array(table)) {
        return(array(
            lapply(table, mT.setDimInfo, ..., append = append),
            dim = dim(table), dimnames = dimnames(table))
        )
296 297
    } else if (is.list(table)) {
        return(lapply(table, mT.setDimInfo, ..., append = append))
298
    }
299
    if (!is(table, "mortalityTable"))
300
        stop("First argument must be a mortalityTable or a list of mortalityTable objects.")
301 302 303 304 305 306 307 308 309

    if (append) {
        table@data$dim[names(list(...))] = list(...)
    } else {
        table@data$dim = list(...)
    }
    table
}

310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328

#' @export
pT.getSubTable = function(table, subtable = "qx") {
    if (is.array(table)) {
        return(array(
            lapply(table, pT.getSubTable, subtable = subtable),
            dim = dim(table), dimnames = dimnames(table))
        )
    } else if (is.list(table)) {
        return(lapply(table, pT.getSubTable, subtable = subtable))
    }
    if (!is(table, "pensionTable"))
        stop("First argument must be a pensionTable or a list of pensionTable objects.")

    if (.hasSlot(table, subtable))
        slot(table, subtable)
    else
        NULL
}
329