Difference between revisions of "Consumption of Baltic herring in Finland"

# This is code Op_en7403/answer on page [[Consumption of Baltic herring in Finland]]
library(OpasnetUtils)
library(ggplot2)
objects.latest("Op_en6007", code_name = "answer") # OpasnetUtils/Drafts
objects.latest("Op_en7403", code_name = "amount") # [[Consumption of Baltic herring in Finland]]

amount <- EvalOutput(amount)
amount$Age_group <- cut(as.numeric(as.character(amount$Ikä)), c(0, 20, 65, 80), include.lowest = TRUE)

ggplot(amount@output, aes(x = amountResult, colour = Age_group)) + stat_ecdf(size = 2) + 
  facet_wrap( ~ Sukupuoli) + scale_x_log10() + theme_gray(base_size = 24) +
  labs(
	title = "Baltic herring consumption by sex and age", 
	colour = "Age group", 
	x = "Consumption (g/d)", 
	y = "Cumulative probability distribution"
  )

#This code is Op_en7403/amount on page [[Consumption of Baltic herring in Finland]]
library(OpasnetUtils)

##########################################################################
# Questionnaire data about Baltic herring

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++++
silakka <- opbase.data("Op_fi3831", subset = "Silakka") # [[:op_fi:Silakan hyöty-riskiarvio]]
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++++

########## PREPROSESSING

silakka$Paino[silakka$Paino < 4 & silakka$Ikäryhmä == "Aikuinen"] <- 75 #Tehdään karkea inputointi aikuisten painoon
silakka$Paino <- ifelse(silakka$Paino < 4 & silakka$Ikäryhmä == "Lapsi", 5+((60-5)/(15)*silakka$Ikä), silakka$Paino) 
# Lineaarinen ekstrapolaatio 5-60 kg

colnames(silakka)[colnames(silakka) == "Nro"] <- "Rivi"
levels(silakka$Result)[levels(silakka$Result) == "Ei syö silakkaa ollenkaan"] <- "En syö silakkaa ollenkaan"

rannikko <- c(
	"Uusimaa",
	"Pohjanmaa", 
	"Kymenlaakso",
	"Etelä-Pohjanmaa",
	"Satakunta",
	"Keski-Pohjanmaa",
	"Pohjois-Pohjanmaa", 
	"Varsinais-Suomi"
)

sisämaa <- c(
	"Kanta-Häme",
	"Pirkanmaa",
	"Etelä-Karjala",
	"Pohjois-Savo",
	"Pohjois-Karjala",
	"Etelä-Savo",
	"Keski-Suomi",
	"Päijät-Häme",
	"Lappi",
	"Kainuu"
)

silakka$Rannikko <- ifelse(silakka$Maakunta %in% rannikko, "Rannikko", "Sisämaa")

ages <- factor(c(
	"0", "1-4", "5-9", "10-14", "15-19", "20-24",
	"25-29", "30-34", "35-39", "40-44", "45-49", "50-54", 
	"55-59", "60-64", "65-69", "70-74", "75-79", "80-84", 
	"85-89", "90-94", "95-"), ordered = TRUE
)

silakka$Age <- cut(silakka$Ikä, breaks = c(0,1,(1:20)*5), labels = ages, right = FALSE)
silakka$Hedelm <- silakka$Ikä >= 20 & silakka$Ikä < 45 # Onko henkilö hedelmällisessä iässä?

# lyhyt = lyhyt lista yksilökohtaisia määrittelyjä eli vain välttämättömät.
lyhyt <- silakka[c("Age", "Ikä", "Hedelm", "Sukupuoli", "Rivi", "Maakunta", "Rannikko")] 

kokoaik	<- Ovariable("kokoaik", data = data.frame(
	lyhyt, 
	Arvo = silakka$Result, 
	Result = silakka$Kokosilakka
)) # Kokonaisten silakoiden syönti kertaa/3 kk

osanaik	<- Ovariable("osanaik", data = data.frame(
	lyhyt,
	Result = silakka$Silakkaruoka
)) # Silakkaruokien syönti kertaa/3 kk

lisuaik	<- Ovariable("lisuaik", data = data.frame(
	lyhyt,
	Result = silakka$Silakkalisuke
)) # Silakkalisukkeiden syönti kertaa/3 kk

## Arvo yksi Iter jokaiselle kyselyn ihmiselle.

iterit <- Ovariable("iterit", data = data.frame(
	Iter = sample(get("N", envir = openv), nrow(silakka), replace = TRUE),
	Rivi = rownames(silakka),
	Result = 1
))

#### Arvo N riviä ehdot täyttävästä kyselyn osaryhmästä (tässä tapauksessa Hedelm+Sukupuoli-ryhmästä)

ehto <- unique(silakka[c("Sukupuoli", "Hedelm")]) # , "Age")]) Säästetään muistia #, "Maakunta")]) Eiköhän painokerroin huolehdi maakunnan

rivit <- data.frame()
for(i in 1:nrow(ehto)) {
	temp <- silakka[
		silakka$Sukupuoli == ehto$Sukupuoli[i] & silakka$Hedelm == ehto$Hedelm[i] , # Jätetään tästäkin Age pois ja Hedelm tilalle
		c("Sukupuoli", "Hedelm", "Painokerroin")
	] # Eikö tässä voisi yksinkertaistaa ja käyttää vain rivinumeroa eikä koko tempiä?
	temp <- temp[sample(1:nrow(temp), get("N", envir = openv), replace = TRUE, prob = temp$Painokerroin) , ]
	rivit <- rbind(rivit, data.frame(
		ehto[i , ],
		Iter = 1:get("N", envir = openv),
		Rivi = as.character(floor(as.numeric(rownames(temp)))),
		Result = 1
	))
}
rivit <- Ovariable(output = rivit, marginal = c(TRUE, TRUE, TRUE, TRUE, FALSE))

BW		<- rivit * Ovariable("BW", data = data.frame(lyhyt, Result = silakka$Paino)) # Ruumiinpaino
BW <- unkeep(BW, cols = c("Rivi", "Ikä"), prevresults = TRUE, sources = TRUE)

######################################################################################
# SILAKOIDEN ANNOSKOOT

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++
solet <- opbase.data("Op_fi3831.saantioletukset") # Silakkaoletukset sivulta [[:op_fi:Silakan hyöty-riskiarvio]]
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++

silakoita	<- Ovariable("silakoita", data = solet[solet$Muuttuja == "V95" , c("Arvo", "Result")]) # Silakoita per silakka-annos

silakanpaino	<- Ovariable("silakanpaino", data = solet[solet$Muuttuja == "Koko silakan paino" , ]["Result"])

raakaainepaino	<- Ovariable("raakaainepaino", data = data.frame(
	Age = ages,
	Result = c(
		rep(solet$Result[solet$Muuttuja == "Raaka-ainesilakan paino, lapset"], 4), 
		rep(solet$Result[solet$Muuttuja == "Raaka-ainesilakan paino"], 17)
	)
))

lisukepaino	<- Ovariable("lisukepaino", data = solet[solet$Muuttuja == "Lisukesilakan paino" , "Result", drop = FALSE])

#!!++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Silakkaruokia (kpl/3 kk)
kouluruokailut <- Ovariable("kouluruokailut", data = data.frame(Result = "0;0;0.8;1.6"))

# Silakka-annoksen koko (g)
annos <- Ovariable("koulut", data = data.frame(Age = ages, Result = c(
	0,
	"20; 39; 40; 50; 60; 70; 80; 80; 100; 100; 100",
	"20; 39; 50; 65; 70; 70; 70; 80; 100; 120",
	"20; 39; 50; 65; 70; 70; 70; 80; 100; 120",
	"30; 39; 50; 80; 100; 110; 120; 120; 130; 160",
	rep(0, 16)))
)

# Syökö lapsi silakka-annoksensa vai ei?
into <- Ovariable("into", data = data.frame(Age = ages, Result = c(
	0,
	"1;1;1;1;1;1;1;1;1;1;0",
	"1;1;1;1;0",
	"1;1;1;1;0",
	"1;1;1;1;1;1;1;0;0;0",
	rep(0, 16)))
)
#ii++++++++++++++++++++++++++++++++++++++++++++++++++++

amount <- Ovariable("amount", 
	dependencies = data.frame(Name = c(
		"kokoaik",
		"silakoita",
		"silakanpaino",
		"osanaik",
		"raakaainepaino",
		"lisuaik",
		"lisukepaino",
		"kouluruokailut",
		"annos",
		"into"
	)),
	formula = function(...) {

		out <- (kokoaik * silakoita * silakanpaino + osanaik * raakaainepaino + lisuaik * lisukepaino) / 91 
		# Per 3 kk -> per d

		out <- out + kouluruokailut * annos * into / 91

		# Muutetaan Age ja Maakunta epävarmaksi eli ei-marginaaliksi
		colnames(out@output)[colnames(out@output) == "Ikä"] <- "Ikä"
		out@marginal[colnames(out@output) %in% c("Age", "Ikä", "Maakunta", "Rannikko")] <- FALSE 
		# Sukupuoli ja Hedelm pidetään indekseinä koska niiden mukaan arvottiin
		out <- unkeep(out, cols = c("Rivi", "Arvo"), prevresults = TRUE, sources = TRUE)

		result(out)[result(out) == 0] <- 0.01 # Ei jätetä nollia saantiin
		return(out)
	}
)

########## Luodaan vaikutusarviointimallia varten ovariable, jossa arvotut yksilöt

riv <- rivit@output
riv$Result <- NULL

kokoaik@data <- merge(kokoaik@data, riv)
osanaik@data <- merge(osanaik@data, riv)
lisuaik@data <- merge(lisuaik@data, riv)

objects.store(
		"kokoaik",
		"silakoita",
		"silakanpaino",
		"osanaik",
		"raakaainepaino",
		"lisuaik",
		"lisukepaino",
		"kouluruokailut",
		"annos",
		"into",
		"amount"
)
cat("Objects
 		kokoaik,
		silakoita,
		silakanpaino,
		osanaik,
		raakaainepaino,
		lisuaik,
		lisukepaino,
		kouluruokailut,
		annos,
		into,
		amount
stored,\n")