Impact Calculation Tool for R: Difference between revisions
No edit summary |
(lifetable function with data.frames works!) |
||
Line 1: | Line 1: | ||
{{method|moderator=Virpi Kollanus}} | |||
=Life table calculations= | =Life table calculations= | ||
==R calculations== | |||
{{defend|# |Lifetable function works with data.frames!|--[[User:Jouni|Jouni]] 23:38, 9 July 2013 (EEST)}} | |||
{{attack|# |But it does not work with ovariables for an unknown colnames reason.|--[[User:Jouni|Jouni]] 23:38, 9 July 2013 (EEST)}} | |||
<rcode> | |||
library(OpasnetUtils) | |||
library(reshape2) | |||
# lifetable is a life table function. | |||
lifetable <- function( | |||
pop = data.frame(Birth = 2000, Result = 1000), # population by birth year. | |||
mort = data.frame(Age = 0:100, Result = 0.02), # mortality by age. | |||
followup = 1:100, # follow-up time starting from year 1. | |||
keep = (0:20) * 5, # Which years to save? | |||
starttime = 2000 # What year is the starting time or reference? | |||
) { | |||
poptemp <- pop | |||
colnames(mort)[colnames(mort) == "Result"] <- "Mort" | |||
# mort <- fillna(orbind(mort, data.frame(Age = -100:-1, Result = 0))) # There could | |||
# be an automatic protection for unborn people but that would mess with indices. | |||
for(i in followup) { | |||
poptemp$Age <- starttime - poptemp$Birth + i | |||
poptemp <- merge(poptemp, mort, all.x = TRUE) | |||
poptemp$Mort[is.na(poptemp$Mort)] <- 1 | |||
if(i == followup[1]) {out <- data.frame("V1" = poptemp$Result)} | |||
poptemp$Result <- poptemp$Result * (1 - poptemp$Mort) | |||
poptemp <- poptemp[colnames(poptemp) != "Mort"] | |||
if(i %in% keep) { | |||
out[[match(i, keep)]] <- poptemp$Result | |||
} | |||
} | |||
indices <- colnames(pop)[colnames(pop) != "Result"] | |||
out <- cbind(poptemp[indices], out) | |||
out <- melt(out, id.vars = indices, variable.name = "Time") | |||
out$Time <- starttime + keep[as.numeric(substr(out$Time, 2, 100))] | |||
return(out) | |||
} | |||
setGeneric("lifetable") | |||
setMethod( | |||
f = "lifetable", | |||
signature = signature(pop = "ovariable", mort = "ovariable"), | |||
definition = function( | |||
pop, | |||
mort, | |||
followup = 1:100, | |||
keep = (0:20) * 5, | |||
starttime = 2000 | |||
) { | |||
v = FALSE | |||
if (ncol(pop@output) == 0) pop <- EvalOutput(pop, verbose = v) | |||
rescolpop <- paste(pop@name, "Result", sep = "") | |||
colnames(pop@output)[colnames(pop@output) == rescolpop] <- "Result" | |||
# Here we should remove all non-index columns. But so far the user has to do it by hand. | |||
if (ncol(mort@output) == 0) mort <- EvalOutput(mort, verbose = v) | |||
rescolmort <- paste(mort@name, "Result", sep = "") | |||
colnames(pop)[colnames(pop) == rescolmort] <- "Result" | |||
callGeneric( | |||
pop = pop@output, | |||
mort = mort@output, | |||
followup = followup, | |||
keep = keep, | |||
starttime = starttime | |||
) | |||
} | |||
) | |||
pop <- opbase.data("Op_en5994", subset = "population") | |||
oprint(head(pop)) | |||
mort <- opbase.data("Op_en5994", subset = "mortality") | |||
oprint(head(mort)) | |||
lifetable(mort = data.frame(Age = 0:100, Result = 0.05)) | |||
pop <- EvalOutput(Ovariable(name = "pop", ddata = "Op_en5994.population")) | |||
pop@output$Age <- as.integer(levels(pop@output$Age)[pop@output$Age]) | |||
mort <- EvalOutput(Ovariable(name = "mort", ddata = "Op_en5994.mortality")) | |||
mort@output$Age <- as.integer(levels(mort@output$Age)[mort@output$Age]) | |||
mort <- mort / pop | |||
mort@output <- mort@output[c("Age", "Result")] | |||
oprint(mort) | |||
pop@output$Birth <- 2010 - pop@output$Age | |||
pop@output <- pop@output[c("Birth", "popResult")] | |||
oprint(pop) | |||
out <- lifetable(pop = pop, mort = mort, starttime = 2012, followup = 1:10) | |||
oprint(out) | |||
</rcode> | |||
==Input data== | ==Input data== |
Revision as of 20:38, 9 July 2013
[show] |
---|
Life table calculations
R calculations
←--#: . Lifetable function works with data.frames! --Jouni 23:38, 9 July 2013 (EEST) (type: truth; paradigms: science: defence)
⇤--#: . But it does not work with ovariables for an unknown colnames reason. --Jouni 23:38, 9 July 2013 (EEST) (type: truth; paradigms: science: attack)
Input data
Population structure in the beginning of the assessment follow-up period (pop_data)
Obs | Age | Result | Kuvaus |
---|---|---|---|
1 | 0 | 56683 | |
2 | 1 | 56683 | |
3 | 2 | 56683 | |
4 | 3 | 5668 | |
5 | 4 | 56683 | |
6 | 5 | 60615 | |
7 | 6 | 60615 | |
8 | 7 | 60615 | |
9 | 8 | 60615 | |
10 | 9 | 60615 | |
11 | 10 | 66167 | |
12 | 11 | 66167 | |
13 | 12 | 66167 | |
14 | 13 | 6616 | |
15 | 14 | 66167 | |
16 | 15 | 63786 | |
17 | 16 | 63786 | |
18 | 17 | 63786 | |
19 | 18 | 63786 | |
20 | 19 | 63786 | |
21 | 20 | 66423 | |
22 | 21 | 66423 | |
23 | 22 | 66423 | |
24 | 23 | 66424 | |
25 | 24 | 66423 | |
26 | 25 | 65882 | |
27 | 26 | 65882 | |
28 | 27 | 65882 | |
29 | 28 | 65882 | |
30 | 29 | 65882 | |
31 | 30 | 61495 | |
32 | 31 | 61495 | |
33 | 32 | 61495 | |
34 | 33 | 61495 | |
35 | 34 | 61495 | |
36 | 35 | 72474 | |
37 | 36 | 72474 | |
38 | 37 | 72474 | |
39 | 38 | 72474 | |
40 | 39 | 72474 | |
41 | 40 | 75917 | |
42 | 41 | 75917 | |
43 | 42 | 75917 | |
44 | 43 | 75917 | |
45 | 44 | 75917 | |
46 | 45 | 76977 | |
47 | 46 | 76977 | |
48 | 47 | 76977 | |
49 | 48 | 76977 | |
50 | 49 | 76977 | |
51 | 50 | 80206 | |
52 | 51 | 80206 | |
53 | 52 | 80206 | |
54 | 53 | 80206 | |
55 | 54 | 80206 | |
56 | 55 | 80291 | |
57 | 56 | 80291 | |
58 | 57 | 80291 | |
59 | 58 | 80291 | |
60 | 59 | 80291 | |
61 | 60 | 54300 | |
62 | 61 | 54300 | |
63 | 62 | 54300 | |
64 | 63 | 54300 | |
65 | 64 | 54300 | |
66 | 65 | 48077 | |
67 | 66 | 48077 | |
68 | 67 | 48077 | |
69 | 68 | 48077 | |
70 | 69 | 48077 | |
71 | 70 | 41475 | |
72 | 71 | 41475 | |
73 | 72 | 41475 | |
74 | 73 | 41475 | |
75 | 74 | 41475 | |
76 | 75 | 34987 | |
77 | 76 | 34987 | |
78 | 77 | 34987 | |
79 | 78 | 34987 | |
80 | 79 | 34987 | |
81 | 80 | 23300 | |
82 | 81 | 23300 | |
83 | 82 | 23300 | |
84 | 83 | 23300 | |
85 | 84 | 23300 | |
86 | 85 | 11292 | |
87 | 86 | 11292 | |
88 | 87 | 11292 | |
89 | 88 | 11292 | |
90 | 89 | 11292 | |
91 | 90 | 4394 | |
92 | 91 | 4394 | |
93 | 92 | 4394 | |
94 | 93 | 4394 | |
95 | 94 | 4394 | |
96 | 95 | 886 | |
97 | 96 | 886 | |
98 | 97 | 886 | |
99 | 98 | 886 | |
100 | 99 | 886 |
Annual birth rate (birth_rate)
Obs | Follow-up period | Result | Kuvaus |
---|---|---|---|
1 | 2010 | 57000 | |
2 | 2011 | 57000 | |
3 | 2012 | 57000 | |
4 | 2013 | 57000 | |
5 | 2014 | 57000 | |
6 | 2015 | 57000 | |
7 | 2016 | 57000 | |
8 | 2017 | 57000 | |
9 | 2018 | 57000 | |
10 | 2019 | 57000 | |
11 | 2020 | 57000 | |
12 | 2021 | 57000 | |
13 | 2022 | 57000 | |
14 | 2023 | 57000 | |
15 | 2024 | 57000 | |
16 | 2025 | 57000 | |
17 | 2026 | 57000 | |
18 | 2027 | 57000 | |
19 | 2028 | 57000 | |
20 | 2029 | 57000 |
Annual mortality rate (mort_rate)
Obs | Age | Result | Kuvaus |
---|---|---|---|
1 | 0 | 49.8 | |
2 | 1 | 49.8 | |
3 | 2 | 49.8 | |
4 | 3 | 49.8 | |
5 | 4 | 49.8 | |
6 | 5 | 10.2 | |
7 | 6 | 10.2 | |
8 | 7 | 10.2 | |
9 | 8 | 10.2 | |
10 | 9 | 10.2 | |
11 | 10 | 10.6 | |
12 | 11 | 10.6 | |
13 | 12 | 10.6 | |
14 | 13 | 10.6 | |
15 | 14 | 10.6 | |
16 | 15 | 30.6 | |
17 | 16 | 30.6 | |
18 | 17 | 30.6 | |
19 | 18 | 30.6 | |
20 | 19 | 30.6 | |
21 | 20 | 50 | |
22 | 21 | 50 | |
23 | 22 | 50 | |
24 | 23 | 50 | |
25 | 24 | 50 | |
26 | 25 | 47 | |
27 | 26 | 47 | |
28 | 27 | 47 | |
29 | 28 | 47 | |
30 | 29 | 47 | |
31 | 30 | 48.6 | |
32 | 31 | 48.6 | |
33 | 32 | 48.6 | |
34 | 33 | 48.6 | |
35 | 34 | 48.6 | |
36 | 35 | 99.4 | |
37 | 36 | 99.4 | |
38 | 37 | 99.4 | |
39 | 38 | 99.4 | |
40 | 39 | 99.4 | |
41 | 40 | 156.6 | |
42 | 41 | 156.6 | |
43 | 42 | 156.6 | |
44 | 43 | 156.6 | |
45 | 44 | 156.6 | |
46 | 45 | 247.6 | |
47 | 46 | 247.6 | |
48 | 47 | 247.6 | |
49 | 48 | 247.6 | |
50 | 49 | 247.6 | |
51 | 50 | 395.8 | |
52 | 51 | 395.8 | |
53 | 52 | 395.8 | |
54 | 53 | 395.8 | |
55 | 54 | 395.8 | |
56 | 55 | 555 | |
57 | 56 | 555 | |
58 | 57 | 555 | |
59 | 58 | 555 | |
60 | 59 | 555 | |
61 | 60 | 534 | |
62 | 61 | 534 | |
63 | 62 | 534 | |
64 | 63 | 534 | |
65 | 64 | 534 | |
66 | 65 | 702.4 | |
67 | 66 | 702.4 | |
68 | 67 | 702.4 | |
69 | 68 | 702.4 | |
70 | 69 | 702.4 | |
71 | 70 | 975 | |
72 | 71 | 975 | |
73 | 72 | 975 | |
74 | 73 | 975 | |
75 | 74 | 975 | |
76 | 75 | 1372.4 | |
77 | 76 | 1372.4 | |
78 | 77 | 1372.4 | |
79 | 78 | 1372.4 | |
80 | 79 | 1372.4 | |
81 | 80 | 1619.6 | |
82 | 81 | 1619.6 | |
83 | 82 | 1619.6 | |
84 | 83 | 1619.6 | |
85 | 84 | 1619.6 | |
86 | 85 | 1399.8 | |
87 | 86 | 1399.8 | |
88 | 87 | 1399.8 | |
89 | 88 | 1399.8 | |
90 | 89 | 1399.8 | |
91 | 90 | 934 | |
92 | 91 | 934 | |
93 | 92 | 934 | |
94 | 93 | 934 | |
95 | 94 | 934 | |
96 | 95 | 313 | |
97 | 96 | 313 | |
98 | 97 | 313 | |
99 | 98 | 313 | |
100 | 99 | 313 |
Mortality risk (mort_risk)
mort_rate / pop_data
Start year (start-year)
2010
Follow-up time in years (followup_time)
20
Analytica codes
Variables, which need to be translated into ovariables:
Population in time, child (pop_in_time_child)
var k: Birth_rate[Fu_year=Year_lt];
k:= if k = null then 0 else k;
var a:= if @Year_lt = 1 then Pop_data else (if @Age=1 then k else 0);
a:= a[Age=age_child];
var j:= Mort_risk[Age=age_child];
j:=j[Fu_period=Period_lt];
j:= Si_pi(j, 5, Period_lt, Year_lt, Year_help)*5;
j:= if j = null then j[Period_lt=max(Fu_period)] else j;
j:= if j < 0 then 0 else j;
j:= if j > 1 then 1 else j;
j:= 1-j;
var x:= 1;
while x<= min([size(age_child),size(Year_lt)]) do (
var b:= a*j; b:= b[@age_child=@age_child-1, @Year_lt=@Year_lt-1];
a:= if b=null then a else b;
x:= x+1);
sum(if Year_lt = period_vs_year then a else 0,Year_lt)
Population in time, beginning of time step (pop_in_time_beg)
var a:= sum(if floor(Age/5)+1 = @Age_cat then Pop_data else 0 , Age);
a:= if @Age_cat=1 then sum(Pop_in_time_child, Age_child) else (if @period_lt = 1 then a else 0);
var j:= Mort_risk;
j:=j[Fu_period=Period_lt];
j:= if j = null then j[Period_lt=max(Fu_period)] else j;
j:= if j < 0 then 0 else j;
j:= if j > 1 then 1 else j;
j:= 1-j;
j:= sum(if floor(Age/5)+1 = @Age_cat then j else 0 , Age)/5;
var m:=j[@Age_cat=@Age_cat+1];
m:= if m=null then 0 else m;
var n:=((j^5)+(j^4*m)+(j^3*m^2)+(j^2*m^3)+(j*m^4))/5;
var x:= 1;
while x<= min([size(Age_cat),size(Period_lt)]) do ( var b:= a*n;
b:= b[@Age_cat=@Age_cat-1, @Period_lt=@Period_lt-1];
a:= if b=null then a else b;
x:= x+1);
a
Indices and function used in the code above:
Follow-up year (fu_year)
sequence(Start_year,Start_year+(Followup_time-1),1)
Year in life table (year_lt)
sequence(Start_year,Start_year+Followup_time+99,1)
Follow-up period in 5-year time steps (fu_period)
sequence(Start_year,Start_year+(Followup_time-1),5)
5-year period in life table (period_lt)
sequence(Start_year,Start_year+Followup_time+99,5)
Age of child (age_child)
sequence(0,4,1)
Si_pi function (si_pi)
Parameters: (data, kerroin;karkea,tarkka:indextype;indtieto)
Description
- Data = data to be divided into more detailed parts
- Kerroin = relative weight inside a cluster
- Karkea = index for the clustered data
- Tarkka = index for the detailed data
- Indtieto = Data about which detailed item belongs to which cluster
Analytica code:
var a:=sum((if indtieto=karkea then kerroin else 0), tarkka);
a:= sum((if indtieto=karkea then a else 0), karkea);
a:= kerroin/a;
sum((if indtieto=karkea then data*a else 0), karkea)