Impact Calculation Tool for R

[show] This page is a knowledge crystal of subtype method. The page identifier is Op_en5994
Moderator:Virpi Kollanus (see all)
Upload data {{#opasnet_base_link:Op_en5994}}

Life table calculations

R calculations

←--#: . Lifetable function works with data.frames! --Jouni 23:38, 9 July 2013 (EEST) (type: truth; paradigms: science: defence) ←--#: . And with ovariables! --Jouni 00:49, 10 July 2013 (EEST) (type: truth; paradigms: science: defence)

+ Show code- Hide code library(OpasnetUtils) library(reshape2) library(ggplot2) # 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", value.name = "Result") 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(mort@output)[colnames(mort@output) == rescolmort] <- "Result" callGeneric( pop = pop@output, mort = mort@output, followup = followup, keep = keep, starttime = starttime ) } ) pop <- opbase.data("Op_en5994", subset = "population") mort <- opbase.data("Op_en5994", subset = "mortality") 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[colnames(output) %in% c("Age", "Result", "Iter")] cat("Mortality rate by age.\n") oprint(mort) pop@output$Birth <- 2012 - pop@output$Age pop@output <- pop@output[colnames(output) %in% c("Birth", "popResult", "Iter")] cat("Population data for year 2012.\n") oprint(pop) out <- lifetable(pop = pop, mort = mort, starttime = 2012, followup = 1:50, keep = c(0, 5, 10, 20, 30, 40, 50)) oprint(head(out)) ggplot(out, aes(x = Time, y = Result, group = Birth)) + geom_line(aes(colour = Birth)) + theme_grey(base_size = 24)

Input data

Population structure in the beginning of the assessment follow-up period (pop_data)

population(-)

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)

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)

mortality(-)

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)