Difference between revisions of "Goherr: Fish consumption study"

  
library(OpasnetUtils)
library(ggplot2)
library(reshape2)

objects.latest("Op_en7749", "preprocess")

for(i in c(5:6, 16, 29:30, 46:49, 85:86, 95:98, 135)) {
  temp <- survey[!is.na(survey[[i]]),]
    p <- ggplot(temp, aes(x = temp[[i]])) + geom_bar() +
    theme_gray(base_size = 18) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
    labs(title = colnames(temp[i])) + xlab("") + facet_wrap(~ Country)
  print(p)
}

temp <- melt(survey, measure.vars = 145:153, variable.name = "Changing.factor", value.name = "Impact")
levs <- c("-5 strongly disagree", "-4", "-3", "-2", "-1", "0 Neutral", "1", "2", "3", "4", "5 strongly agree", "I don't know")
temp$Impact <- factor(temp$Impact, levels = levs, labels = levs, ordered = TRUE)

ggplot(temp, aes(x = Impact, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Changing.factor)+
  theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4))

surveytemp <- subset(survey, survey[[31]] == "Yes")

temp <- melt(surveytemp, measure.vars = 38:43, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + labs(title = "Baltic salmon sources")

temp <- melt(surveytemp, measure.vars = 50:59, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + labs(title = "Reasons to eat Baltic salmon")

temp <- melt(surveytemp, measure.vars = 73:82, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + labs(title = "Baltic salmon actions")

surveytemp <- subset(survey, survey[[31]] == "No")

temp <- melt(surveytemp, measure.vars = 62:70, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + labs(title = "Reasons not to eat Baltic salmon")


surveytemp <- subset(survey, survey[[86]] == "Yes")

temp <- melt(surveytemp, measure.vars = 87:92, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + labs(title = "Baltic herring sources")

temp <- melt(surveytemp, measure.vars = 99:108, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + labs(title = "Reasons to eat Baltic herring")

temp <- melt(surveytemp, measure.vars = 123:132, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + labs(title = "Baltic herring actions")

surveytemp <- subset(survey, survey[[86]] == "No")

temp <- melt(surveytemp, measure.vars = 111:120, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
   theme_gray(base_size = 24) + labs(title = "Reasons not to eat Baltic herring")

library(OpasnetUtils)

############## Generic functions and objects are defined first.

### webropol.convert converts a csv file from Webropol into a useful data.frame.

webropol.convert <- function(
  data, # Data.frame created from a Webropol csv file. The first row should contain headings.
  rowfact, # Row number where the factor levels start (in practice, last row + 3)
  textmark = "Other open" # The text that is shown in the heading if there is an open sub-question.
) {
  out <- dropall(data[2:(rowfact - 3) , ])
  subquestion <- t(data[1 , ])
  subquestion <- gsub("\xa0", " ", subquestion)
  subquestion <- gsub("\xb4", " ", subquestion)
  subquestion <- gsub("\n", " ", subquestion)
  #  subquestion <- gsub("\\(", " ", subquestion)
  #  subquestion <- gsub("\\)", " ", subquestion)
  textfield <- regexpr(textmark, subquestion) != -1
  subquestion <- strsplit(subquestion, ":") # Divide the heading into a main question and a subquestion.
  subqtest <- 0 # The previous question name.
  for(i in 1:ncol(out)) {
    #print(i)
    if(subquestion[[i]][1] != subqtest) { # If part of previous question, use previous fact.
      fact <- as.character(data[rowfact:nrow(data) , i]) # Create factor levels from the end of Webropol file.
      fact <- fact[fact != ""] # Remove empty rows
      fact <- gsub("\xa0", " ", fact)
      fact <- gsub("\xb4", " ", fact)
      fact <- gsub("\n", " ", fact)
      fact <- strsplit(fact, " = ") # Separate value (level) and interpretation (label)
    }
    if(length(fact) != 0 & !textfield[i]) { # Do this only if the column is not a text type column.
      out[[i]] <- factor(
        out[[i]], 
        levels = unlist(lapply(fact, function(x) x[1])), 
        labels = unlist(lapply(fact, function(x) x[2])), 
        ordered = TRUE
      )
    }
    subqtest <- subquestion[[i]][1]
  }
  return(out)
}

# merge.questions takes a multiple checkbox question and merges that into a single factor.
# First levels in levs have priority over others, if several levels apply to a row.

merge.questions <- function(
  dat, # data.frame with questionnaire data
  cols, # list of vectors of column names or numbers to be merged into one level in the factor
  levs, # vector (with the same length as cols) of levels of factors into which questions are merged.
  name # text string for the name of the new factor column in the data.
) {
  for(i in length(cols):1) {
    temp <- FALSE
    for(j in rev(cols[[i]])) {
      temp <- temp | !is.na(dat[[j]])
    }
    dat[[name]][temp] <- levs[i]
  }
  dat[[name]] <- factor(dat[[name]], levels = levs, ordered = TRUE)
  return(dat)
}

############# Data preprocessing

# Get the data either from Opasnet or your own hard drive.

#Survey original file: N:/Ymal/Projects/Goherr/WP5/Goherr_fish_consumption.csv

survey <- opasnet.csv("5/57/Goherr_fish_consumption.csv", sep = ";", fill = TRUE, quote = "\"")
surcol <- t(survey[1,])

survey <- webropol.convert(survey, 2121, textmark = ":Other open") # Data file is converted to data.frame using levels at row 1269.

# Take the relevant columnames from the table on the page.
colframe <- opbase.data("Op_en7749", subset = "Questions in the Goherr questionnaire")
colnames(survey) <- colframe$Result

survey$Row <- 1:nrow(survey)
survey$Weighting <- as.double(levels(survey$Weighting))[survey$Weighting]
survey$Ages <- ifelse(as.numeric(as.character(survey$Age)) < 46, "18-45",">45")

freqlist <- c(
  "less than once a year",
  "A few times a year",
  "1 - 3 times per month",
  "once a week",
  "2 - 4 times per week",
  "5 or more times per week"
)

amlist <- c(
  "1/6 plate or below (50 grams)",
  "1/3 plate (100 grams)",
  "1/2 plate (150 grams)",
  "2/3 plate (200 grams)",
  "5/6 plate (250 grams)",
  "full plate (300 grams)",
  "overly full plate (500 grams)"
#  "Not able to estimate"
)

sidel <- c(
  "1/6 plate or below (50 grams)",
  "1/4 plate (100 grams)",
  "1/2 plate (150 grams)",
  "2/3 plate (200 grams)",
  "5/6 plate (250 grams)"
#  "Not able to estimate"
)

ansl <- list(
  `How often eat fish` = freqlist,
  `How often Baltic salmon` = c("Never", freqlist),
  `How much Baltic salmon` = amlist,
  `How often side Baltic salmon` = c("Never", freqlist),
  `How much side Baltic salmon` = sidel,
  `How often Baltic herring` = c("Never", freqlist),
  `How much Baltic herring` = amlist,
  `How often side Baltic herring` = c("Never", freqlist),
  `How much side Baltic herring` = sidel
)

for (i in names(ansl)) {
  survey[[i]] <- factor(survey[[i]], levels = ansl[[i]])
}

objects.store(survey, surcol)
cat("objects survey, surcol were stored.\n")

# This is code Op_en7749/bayes on page [[Goherr: Fish consumption study]]

library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(rjags)
library(car)

# Fish intake in humans
# Data from data.frame survey from page [[Goherr: Fish consumption study]]
# Start with salmon questions 46:49 (amounts eaten)
# Some preprocessing should be moved away from this code.
# Model assumes that questionnaire data follows binomial distribution with parameter p and max value levels.
# p depends on country (4 groups), sex (2), age (2) that are known individually.
# We should test whether sex makes a difference or whether it can be omitted.
# Predicted p for each question, determined by the group, are produced from the bayes model.
# p is used to sample answers, which are then combined with estimates of amounts related to each answer.
# Total amount eaten of each fish per modelled individual is finally calculated.

objects.latest("Op_en7749", "preprocess") # [[Goherr: Fish consumption study]]: survey, surcol

agel <- as.character(unique(survey$Ages))
countryl <- sort(as.character(unique(survey$Country)))
genderl <- sort(as.character(unique(survey$Gender)))
fisl <- c("Salmon", "Herring")

# Interesting fish eating questions
surv <- survey[c(1,3,158,16,29,30,31,46:49,86,95:98)]
colnames(surv)
#[1] "Country"                       "Gender"                       
#[3] "Ages"                          "Fish eating"                  
#[5] "How often eat fish"            "Salmon eating"                
#[7] "Baltic salmon"                 "How often Baltic salmon"      
#[9] "How much Baltic salmon"        "How often side Baltic salmon" 
#[11] "How much side Baltic salmon"   "Eat Baltic herring"           
#[13] "How often Baltic herring"      "How much Baltic herring"      
#[15] "How often side Baltic herring" "How much side Baltic herring" 
test <- sapply(unique(surv[c(4,7,12)]), unique)#function(x) sum(is.na(x)))

oprint(table(surv[c(12,7,4)], useNA = "ifany"))
oprint(table(surv[c(13,12)], useNA = "ifany"))
# For estimating distributions, we should 
#1 remove people with Fish eating = No (142)
#2 merge Eat Baltic herring = I don't know with No (How often BH = NA always)
#3 merge Baltic salmon = NA with No (because they usually have answered BH questions)

oprint(table(is.na(rowSums(sapply(surv[4:16], as.numeric)))))
# BUT: there are so many missing values, that we just model BH and BS separately now.

surv$Ages <- match(surv$Ages, agel) # Not a factor, coerce to integer
surv <- as.data.frame(lapply(surv, FUN = function(x) as.integer(x))) # Coerce to integers
surv[is.na(surv[[12]]) | surv[[12]] == 3 , 12] <- 1 # Eat Baltic herring: I don't know --> No
surv[4:16] <- lapply(surv[4:16], FUN = function(x) x-1) # Scale to start from zero

# Row numbers for respondents that have eaten fish, Baltic salmon, and Baltic herring
eatfish <- surv[[4]] %in% 1
eatsalm <- !(surv[[7]] %in% 0 | is.na(rowSums(surv[7:11])))
eatherr <- surv[[12]] %in% 1 & !is.na(rowSums(surv[13:16]))

# Oletetaan, että covarianssimatriisi on vakio kaikille maille ja sukupuolille yms
# mutta keskiarvo on spesifi näille ja kysymykselle.

#qlen <-  c(4,2,2,2,6,2,2,7,7,7,5,2,7,7,7,5) # Number of options in each question
# qlen not needed when dbinom is not used.

agel
countryl
genderl
fisl

####################### Descriptive statistics

oprint(cor(surv, use = "pairwise.complete.obs"))
# --> Baltic salmon and herring eating are correlated, so they should be estimated together

############################# Plot original data

## Eating frequencies of fish and Baltic salmon and herring with random noise, all
pl <- surv[eatfish,c(5,8,10,13,15)]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

ggplot(data.frame(
  X = rep(1,5), Y = 5:1,
  legend = c("All", "Finland", "Sweden", "Denmark", "Estonia")
), aes(x=X, y=Y, label=legend))+
  geom_text()+
  labs(title="Fish eating questions with random noise")


## Fish eating questions with some random noise, all
pl <- surv[eatfish,c(5,6,7,12)]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

## Fish eating questions with some random noise, FI
pl <- surv[surv$Country == 1 & eatfish,c(5,6,7,12)]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

## Fish eating questions with some random noise, SWE
pl <- surv[surv$Country == 2 & eatfish,c(5,6,7,12)]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

## Fish eating questions with some random noise, Dk
pl <- surv[surv$Country == 3 & eatfish,c(5,6,7,12)]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

## Fish eating questions with some random noise, EST
pl <- surv[surv$Country == 4 & eatfish,c(5,6,7,12)]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

## Baltic herring questions with some random noise
pl <- surv[eatherr,13:16]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

## Baltic salmon questions with some random noise
pl <- surv[eatsalm,8:11]
scatterplotMatrix(
  data.matrix(pl) + runif(nrow(pl)*ncol(pl), -0.5, 0.5)
)

############################################################
cat("Version with multivariate normal.\n")
ggplot(data.frame(A=1:2, B=1:2), aes(A, y=B))+geom_point()+
  labs(title="Version with multivariate normal.")

mod <- textConnection("
  model{
    for(i in 1:S) {
      survs[i,1:4] ~ dmnorm(mus[], Omegas[,])
    }
    for(j in 1:H) {
      survh[j,1:4] ~ dmnorm(muh[], Omegah[,])
    }
    mus[1:4] ~ dmnorm(mu0[1:4], S2[1:4,1:4])
    Omegas[1:4,1:4] ~ dwish(S3[1:4,1:4],S)
    anss.pred ~ dmnorm(mus[], Omegas[,])
    muh[1:4] ~ dmnorm(mu0[1:4], S2[1:4,1:4])
    Omegah[1:4,1:4] ~ dwish(S3[1:4,1:4],H)
    ansh.pred ~ dmnorm(muh[], Omegah[,])
  }
")

jags <- jags.model(
  mod,
  data = list(
    survs = surv[eatsalm,c(8:11)],
    S = sum(eatsalm),
    survh = surv[eatherr,c(13:16)],
    H = sum(eatherr),
    mu0 = rep(2,4),
    S2 = diag(4)/100000,
    S3 = diag(4)/10000
  ),
  n.chains = 4,
  n.adapt = 100
)

update(jags, 100)
samps.j <- jags.samples(
  jags, 
  c('muh','Omegah', 'ansh.pred', 'mus', 'Omegas', 'anss.pred'), 
  1000
)
jh <- array(
  samps.j$ansh.pred,
  dim = c(4,1000,4),
  dimnames = list(Question = 1:4, Iter = 1:1000, Seed = 1:4)
)
scatterplotMatrix(t(jh[,,1]))

js <- array(
  samps.j$anss.pred,
  dim = c(4,1000,4),
  dimnames = list(Question = 1:4, Iter = 1:1000, Seed = 1:4)
)
scatterplotMatrix(t(js[,,1]))

samps.c <- coda.samples(
  jags, 
  c('muh','Omegah', 'ansh.pred', 'mus', 'Omegas', 'anss.pred'), 
  1000
)
plot(samps.c)

jhm <- melt(jh[,,1], value.name = "Result")
jsm <- melt(js[,,1], value.name = "Result")

oftenh  <- Ovariable("oftenh",  data = jsm[jsm$Question == 1, 2:3])
muchh   <- Ovariable("muchh" ,  data = jsm[jsm$Question == 2, 2:3])
ofsideh <- Ovariable("ofsideh", data = jsm[jsm$Question == 3, 2:3])
musideh <- Ovariable("musideh", data = jsm[jsm$Question == 4, 2:3])
oftens  <- Ovariable("oftens",  data = jhm[jsm$Question == 1, 2:3])
muchs   <- Ovariable("muchs" ,  data = jhm[jsm$Question == 2, 2:3])
ofsides <- Ovariable("ofsides", data = jhm[jsm$Question == 3, 2:3])
musides <- Ovariable("musides", data = jhm[jsm$Question == 4, 2:3])

objects.store(oftenh, muchh, ofsideh, musideh, oftens, muchs, ofsides, musides)
cat("Ovariables oftenh, muchh, ofsideh, musideh, oftens, muchs, ofsides, musides stored.\n")

# This is code Op_en7749/calculate_amount on page [[Goherr: Fish consumption study]]

library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(rjags)

objects.latest("Op_en7749", code_name = "bayes") #: ans.pred, p.pred, p2.pred

ql <- melt(ans.pred)

##Assumptions for amount and frequencies

assumptions <- opbase.data("Op_en7749.assumptions_for_calculations")

salmon_dishtimesassump <- Ovariable("salmon_dishtimesassump", data = assumptions[assumptions$Variable == "freq" , c("value", "Result")])

salmon_dishamountassump <- Ovariable("salmon_dishamountassump", data = assumptions[assumptions$Variable == "amdish" , c("value", "Result")])

salmon_ingridient <- Ovariable("salmon_ingridient", data = assumptions[assumptions$Variable == "ingridient", ]["Result"])

salmon_sidetimesassump <- Ovariable("salmon_wholetimesassump", data = assumptions[assumptions$Variable == "freq" , c("value", "Result")])

salmon_sideamountassump <- Ovariable("salmon_wholeamountassump", data = assumptions[assumptions$Variable == "amside" , c("value", "Result")])

##Build up the variables to calculate amount of Baltic salmon eating

salmon_dishtimes <-  Ovariable("salmon_dishtimes", data = data.frame(
  subset(ql, ql$Fish == "Salmon" & ql$Question == "How often" & ql$Seed == "S1"),
  Result = 1
))

salmon_dishtimes@data$Question <- NULL

salmon_dishamount <- Ovariable("salmon_dishamount", data = data.frame(
  subset(ql, ql$Fish == "Salmon" & ql$Question == "How much" & ql$Seed == "S1"),
  Result = 1
)) #salmon with dish amount per one serving

salmon_dishamount@data$Question <- NULL  

salmon_sidetimes <- Ovariable("salmon_sidetimes", data = data.frame(
  subset(ql, ql$Fish == "Salmon" & ql$Question == "How often side" & ql$Seed == "S1"),
  Result = 1
)) #salmon as such times / year

salmon_sidetimes@data$Question <- NULL

salmon_sideamount <- Ovariable("salmon_sideamount", data = data.frame(
  subset(ql, ql$Fish == "Salmon" & ql$Question == "How much side" & ql$Seed == "S1"),
  Result = 1
)) #salmon as whole amount per one serving

salmon_sideamount@data$Question <- NULL

amount <- Ovariable("amount", 
                    dependencies = data.frame(Name = c(
                      "salmon_dishamount",
                      "salmon_dishamountassump",
                      "salmon_ingridient",
                      "salmon_dishtimes",
                      "salmon_dishtimesassump",
                      "salmon_sideamount",
                      "salmon_sideamountassump",
                      "salmon_sidetimes",
                      "salmon_sidetimesassump"
                    )),
                    formula = function(...) {
                      
                      out1 <- (salmon_dishamount * salmon_dishamountassump * salmon_ingridient) 
                      out1@output$value <- NULL
                      out1 <- ((salmon_dishtimes * salmon_dishtimesassump) * out1) / 365
                      out1@output$value <- NULL
                      out2 <- (salmon_sideamount * salmon_sideamountassump)
                      out2@output$value <- NULL
                      out2 <- ((salmon_sidetimes * salmon_sidetimesassump) * out2) / 365
                      out2@output$value <- NULL
                      out <- out1 + out2
                      
                      # Per vuosi -> per d
                      
                      #result(out)[result(out) == 0] <- 0.01 # Ei jätetä nollia saantiin
                      return(out)
                    }
)

amountSalmon <- EvalOutput(amount, forceEval = TRUE, N = 1000)

temp <- amountSalmon@output 
temp <- aggregate(temp$amountResult, by=list(temp$Country, temp$Age, temp$Gender), FUN = mean)
colnames(temp) <- c("Country","Age","Gender","Salmon amount")
oprint(temp)

##Build up the variables to calculate amount of Baltic herring eating

herring_dishtimes <-  Ovariable("herring_dishtimes", data = data.frame(
  subset(ql, ql$Fish == "Herring" & ql$Question == "How often" & ql$Seed == "S1"),
  Result = 1
))

herring_dishtimes@data$Question <- NULL

herring_dishamount <- Ovariable("herring_dishamount", data = data.frame(
  subset(ql, ql$Fish == "Herring" & ql$Question == "How much" & ql$Seed == "S1"),
  Result = 1
)) #herring with dish amount per one serving

herring_dishamount@data$Question <- NULL  

herring_sidetimes <- Ovariable("herring_sidetimes", data = data.frame(
  subset(ql, ql$Fish == "Herring" & ql$Question == "How often side" & ql$Seed == "S1"),
  Result = 1
)) #herring as such times / year

herring_sidetimes@data$Question <- NULL

herring_sideamount <- Ovariable("herring_sideamount", data = data.frame(
  subset(ql, ql$Fish == "Herring" & ql$Question == "How much side" & ql$Seed == "S1"),
  Result = 1
)) #herring as whole amount per one serving

herring_sideamount@data$Question <- NULL

herring_dishtimesassump <- Ovariable("herring_dishtimesassump", data = assumptions[assumptions$Variable == "freq" , c("value", "Result")])

herring_dishamountassump <- Ovariable("herring_dishamountassump", data = assumptions[assumptions$Variable == "amdish" , c("value", "Result")])

herring_ingridient <- Ovariable("herring_ingridient", data = assumptions[assumptions$Variable == "ingridient", ]["Result"])

herring_sidetimesassump <- Ovariable("herring_wholetimesassump", data = assumptions[assumptions$Variable == "freq" , c("value", "Result")])

herring_sideamountassump <- Ovariable("herring_wholeamountassump", data = assumptions[assumptions$Variable == "amside" , c("value", "Result")])

amount <- Ovariable("amount", 
                    dependencies = data.frame(Name = c(
                      "herring_dishamount",
                      "herring_dishamountassump",
                      "herring_ingridient",
                      "herring_dishtimes",
                      "herring_dishtimesassump",
                      "herring_sideamount",
                      "herring_sideamountassump",
                      "herring_sidetimes",
                      "herring_sidetimesassump"
                    )),
                    formula = function(...) {
                      
                      out1 <- (herring_dishamount * herring_dishamountassump * herring_ingridient) 
                      out1@output$value <- NULL
                      out1 <- ((herring_dishtimes * herring_dishtimesassump) * out1) / 365
                      out1@output$value <- NULL
                      out2 <- (herring_sideamount * herring_sideamountassump)
                      out2@output$value <- NULL
                      out2 <- ((herring_sidetimes * herring_sidetimesassump) * out2) / 365
                      out2@output$value <- NULL
                      out <- out1 + out2
                      
                      # Per vuosi -> per d
                      
                      #result(out)[result(out) == 0] <- 0.01 # Ei jätetä nollia saantiin
                      return(out)
                    }
)

amountHerring <- EvalOutput(amount, forceEval = TRUE, N = 1000)

temp <- amountHerring@output 
temp <- aggregate(temp$amountResult, by=list(temp$Country, temp$Age, temp$Gender), FUN = mean)
colnames(temp) <- c("Country","Age","Gender","Herring amount")
oprint(temp)

objects.store(amountHerring, amountSalmon)
cat("objects amountHerring, amountSalmon were stored.\n")