Commit 6fd5fd4c authored by Renaud Lancelot's avatar Renaud Lancelot 🌍
Browse files

A few changes in function trend(): corrected a mistake in the computation of...

A few changes in function trend(): corrected a mistake in the computation of daily mortality peak, and changes the definition of a few indicators
parent f2955779
......@@ -327,7 +327,7 @@ trend <- function(x){
## 75%, 90% after the peak dr (vector), dr0 (max) daily death rate
Y <- lo <- hi <- r19 <- a <- b <- d2peak <- ldpeak <- NA
a19 <- b19 <- maxgr <- dr0 <- dr <- NA
maxgr <- dr0 <- dr <- NA
T95 <- T90 <- T80 <- T70 <- T60 <- T50 <- NA
T40 <- T30 <- T20 <- T10 <- NA
......@@ -343,6 +343,8 @@ trend <- function(x){
pred <- predict(fm, se = T)
Y <- exp(pred$fit)
## these are counts we want rates
dmr <- 100 * Y / x$pop
## 95% CI
z <- qnorm(p = .975, lower.tail = T)
......@@ -361,16 +363,14 @@ trend <- function(x){
## time from the first reported death
day <- as.numeric(x$date - min(x$date))
## average daily mortality growth rate
admgr <- mean(dY[day < dlock])
l0 <- length(admgr)
## daily mortality growth rate for the early stage of the epidemic
gr11 <- (dY[day == dlock + 11])
l0 <- length(gr11)
## relative change in mortality growth rate, reference value = admgr
rgr11 <- 100 * dY[day = dlock + 11] / admgr
l11 <- length(rgr11)
rgr19 <- 100 * dY[day = dlock + 19] / admgr
## relative change in mortality growth rate, reference value = gr11
rgr19 <- 100 * dY[day = dlock + 19] / gr11
l19 <- length(rgr19)
rgr46 <- 100 * dY[day = dlock + 46] / admgr
rgr46 <- 100 * dY[day = dlock + 46] / gr11
l46 <- length(rgr46)
## negative growth rate
......@@ -392,11 +392,12 @@ trend <- function(x){
l2 <- length(ldpeak)
## max fitted death rate
dr0 <- Y[day == d2peak]
dr0 <- dmr[day == d2peak]
l3 <- length(dr0)
## decay
Y2 <- Y[day >= d2peak]
Ymax <- Y[day == d2peak]
l4 <- length(Y2)
d2 <- day[day >= d2peak] - d2peak
l5 <- length(d2)
......@@ -404,11 +405,13 @@ trend <- function(x){
## collects the results
fittedTrend <- data.frame(
name = x$name,
date = x$date,
day = day,
date1st = rep(min(x$date), nrow(x)),
obsdr = 100 * tapply(x$deaths,
x$day,
sum) / x$pop,
obsdr = 100 * tapply(x$deaths,
x$day,
sum) / x$pop,
pop = unique(x$pop),
lockdown = rep(xlock, nrow(x)),
d2peak = rep(if(l1 > 0) d2peak else NA, nrow(x)),
ldpeak = rep(if(l2 > 0) ldpeak else NA, nrow(x)),
......@@ -421,16 +424,16 @@ trend <- function(x){
## mortality decay
size <- unique(x$pop)
cond10 <- Y2 <= .90 * dr0
cond20 <- Y2 <= .80 * dr0
cond30 <- Y2 <= .70 * dr0
cond40 <- Y2 <= .60 * dr0
cond50 <- Y2 <= .50 * dr0
cond60 <- Y2 <= .40 * dr0
cond70 <- Y2 <= .30 * dr0
cond80 <- Y2 <= .20 * dr0
cond90 <- Y2 <= .10 * dr0
cond95 <- Y2 <= .05 * dr0
cond10 <- Y2 <= .90 * Ymax
cond20 <- Y2 <= .80 * Ymax
cond30 <- Y2 <= .70 * Ymax
cond40 <- Y2 <= .60 * Ymax
cond50 <- Y2 <= .50 * Ymax
cond60 <- Y2 <= .40 * Ymax
cond70 <- Y2 <= .30 * Ymax
cond80 <- Y2 <= .20 * Ymax
cond90 <- Y2 <= .10 * Ymax
cond95 <- Y2 <= .05 * Ymax
fittedParam <- data.frame(
name = unique(x$name),
......@@ -438,8 +441,7 @@ trend <- function(x){
d2peak = if(l1 > 0) d2peak else NA,
ldpeak = if(l2 > 0) ldpeak else NA,
maxdr = if(l3 > 0) dr0 else NA,
admgr = if(l0 > 0) admgr else NA,
rgr11 = if(l11 > 0) rgr11 else NA,
gr11 = if(l0 > 0) gr11 else NA,
rgr19 = if(l19 > 0) rgr19 else NA,
rgr46 = if(l46 > 0) rgr46 else NA,
T10 = if(any(cond10)) min(d2[cond10]) else NA,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment