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)

----#: . Changed to a more versatile ovariable implementation, fixed loop so that no unwanted years show, also more efficient loop (rbind vs merge and melt at end), "pop" renamed "Population" due to some unknown bug, see comparison for details. --Teemu R 22:35, 22 July 2013 (EEST) (type: truth; paradigms: science: comment)

+ Show code- Hide code library(OpasnetUtils) library(ggplot2) ######################### # Lifetable function # ######################### lifetable <- function( pop = data.frame(Birth = 2000, Result = 1000), # population by birth year. mort_frac = data.frame(Age = 0:100, Result = 0.02), # mortality by age. followup = 100, # follow-up time. keep = (0:20) * 5, # Which years to save? starttime = 2000, # What year is the starting time or reference? Ignored if present in data. ... # Ovariable evaluation extra parameters ) { # Check for ovariable input; Instead of having 4 different methods for all combinations we can check # individually in one function. if(class(mort_frac)=="ovariable"){ if(nrow(mort_frac@output) == 0) mort_frac <- EvalOutput(mort_frac, ...) # Proceed to remove non-marginal columns temp <- result(mort_frac) mort_frac <- mort_frac@output[mort_frac@marginal] mort_frac$Time <- NULL mort_frac$Mort <- temp } if(class(pop)=="ovariable"){ if(nrow(pop@output) == 0) pop <- EvalOutput(pop, ...) # Proceed to remove non-marginal columns temp <- result(pop) pop <- pop@output[pop@marginal] pop$Result <- temp } # Prepare input data for yearly loops of follow-up. if(!"Time" %in% colnames(pop)) { pop$Time <- starttime # create Time index if it does not exist } else { starttime <- min(pop$Time) # else take lowest value as starting point } temp_pop <- data.frame() # init out <- pop # Follow the population year by year. for(i in 0:followup) { # Take any newborn population. temp_pop_add <- pop[pop$Time == starttime + i,] # Rbind newborn with leftover population from last iteration temp_pop <- rbind(temp_pop_add, temp_pop) temp_pop$Age <- temp_pop$Time - temp_pop$Birth # update age temp_pop <- merge(temp_pop, mort_frac, all.x = TRUE) # combine population and age-specific mortality temp_pop$Mort[is.na(temp_pop$Mort)] <- 1 # All die beyond the mortality table temp_pop$Result <- temp_pop$Result * (1 - temp_pop$Mort) # update the population size temp_pop <- temp_pop[colnames(pop)] temp_pop$Time <- temp_pop$Time + 1 # update time # Rbind wanted years with output if(i %in% keep) { out <- rbind(out, temp_pop) } } return(out) } ################################# # Ovariable Definitions # ################################# Population <- Ovariable(name = "Population", ddata = "Op_en5994.population") Population@data$Age <- as.integer(levels(Population@data$Age)[Population@data$Age]) Population@data$Time <- as.integer(levels(Population@data$Time)[Population@data$Time]) Population@data$Birth <- Population@data$Time - Population@data$Age mort <- Ovariable(name = "mort", ddata = "Op_en5994.mortality") mort@data$Age <- as.integer(levels(mort@data$Age)[mort@data$Age]) mort_frac <- Ovariable("mort_frac", dependencies = data.frame(Name = c("mort", "Population")), formula = function(...){ temp <- mort / Population temp <- temp@output[temp@output$Time == "2011" , !colnames(temp@output) %in% c("Birth")] return(temp) } ) ################################# # Evaluation and presentation # ################################# out <- lifetable(pop = Population, mort_frac = mort_frac, N = 100) cat("Mortality, Population and mortality rate by age in 2011.\n") oprint(mort_frac) #cat("Population data for the start years of subcohorts.\n") #oprint(Population) cat("Lifetable\n") oprint(out[1:100,]) #oprint(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	Time	Age	Result	Description
1	2011	0	56683
2	2011	1	56683
3	2011	2	56683
4	2011	3	56683
5	2011	4	56683
6	2011	5	60615
7	2011	6	60615
8	2011	7	60615
9	2011	8	60615
10	2011	9	60615
11	2011	10	66167
12	2011	11	66167
13	2011	12	66167
14	2011	13	66167
15	2011	14	66167
16	2011	15	63786
17	2011	16	63786
18	2011	17	63786
19	2011	18	63786
20	2011	19	63786
21	2011	20	66423
22	2011	21	66423
23	2011	22	66423
24	2011	23	66424
25	2011	24	66423
26	2011	25	65882
27	2011	26	65882
28	2011	27	65882
29	2011	28	65882
30	2011	29	65882
31	2011	30	61495
32	2011	31	61495
33	2011	32	61495
34	2011	33	61495
35	2011	34	61495
36	2011	35	72474
37	2011	36	72474
38	2011	37	72474
39	2011	38	72474
40	2011	39	72474
41	2011	40	75917
42	2011	41	75917
43	2011	42	75917
44	2011	43	75917
45	2011	44	75917
46	2011	45	76977
47	2011	46	76977
48	2011	47	76977
49	2011	48	76977
50	2011	49	76977
51	2011	50	80206
52	2011	51	80206
53	2011	52	80206
54	2011	53	80206
55	2011	54	80206
56	2011	55	80291
57	2011	56	80291
58	2011	57	80291
59	2011	58	80291
60	2011	59	80291
61	2011	60	54300
62	2011	61	54300
63	2011	62	54300
64	2011	63	54300
65	2011	64	54300
66	2011	65	48077
67	2011	66	48077
68	2011	67	48077
69	2011	68	48077
70	2011	69	48077
71	2011	70	41475
72	2011	71	41475
73	2011	72	41475
74	2011	73	41475
75	2011	74	41475
76	2011	75	34987
77	2011	76	34987
78	2011	77	34987
79	2011	78	34987
80	2011	79	34987
81	2011	80	23300
82	2011	81	23300
83	2011	82	23300
84	2011	83	23300
85	2011	84	23300
86	2011	85	11292
87	2011	86	11292
88	2011	87	11292
89	2011	88	11292
90	2011	89	11292
91	2011	90	4394
92	2011	91	4394
93	2011	92	4394
94	2011	93	4394
95	2011	94	4394
96	2011	95	886
97	2011	96	886
98	2011	97	886
99	2011	98	886
100	2011	99	886
101	2012	0	57000
102	2013	0	57000
103	2014	0	57000
104	2015	0	57000
105	2016	0	57000
106	2017	0	57000
107	2018	0	57000
108	2019	0	57000
109	2020	0	57000
110	2021	0	57000
111	2022	0	57000
112	2023	0	57000
113	2024	0	57000
114	2025	0	57000
115	2026	0	57000
116	2027	0	57000
117	2028	0	57000
118	2029	0	57000

Annual birth rate (birth_rate)

⇤--#: . Not needed because it can be given as a part of the population. --Jouni 06:46, 10 July 2013 (EEST) (type: truth; paradigms: science: attack)

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)