Goherr: Fish consumption study: Difference between revisions

From Opasnet
Jump to navigation Jump to search
(→‎Figures for the first manuscript: all (?) graphs created)
(added references)
 
(30 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{study|moderator=Arja|stub=Yes}}
{{progression class
|progression=Reviewed
|curator=THL
|date=2019-08-26
}}
{{study|moderator=Arja}}
 
:''This page contains a detailed description about data management and analysis of an international survey related to scientific article ''Forage Fish as Food: Consumer Perceptions on Baltic Herring'' by Mia Pihlajamäki, Arja Asikainen, Suvi Ignatius, Päivi Haapasaari, and Jouni T. Tuomisto.<ref name="pihlajamaki2019">Mia Pihlajamäki, Arja Asikainen, Suvi Ignatius, Päivi Haapasaari and Jouni T. Tuomisto. Forage Fish as Food: Consumer Perceptions on Baltic Herring. Sustainability 2019, 11(16), 4298; https://doi.org/10.3390/su11164298</ref> The results of this survey where also used in another article ''Health effects of nutrients and environmental pollutants in Baltic herring and salmon: a quantitative benefit-risk assessment'' by the same group.<ref name="tuomisto2020">Tuomisto, J.T., Asikainen, A., Meriläinen, P., Haapasaari, P. Health effects of nutrients and environmental pollutants in Baltic herring and salmon: a quantitative benefit-risk assessment. BMC Public Health 20, 64 (2020). https://doi.org/10.1186/s12889-019-8094-1</ref>


== Question ==
== Question ==
Line 6: Line 13:
== Answer ==
== Answer ==


Original questionnaire analysis results  
* Model run with all the results of the article<ref name="pihlajamaki2019"/> [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=B2jMOHfuUSmTjfrn 26.8.2018]
*[http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QaMJZqUX0cPaTfOF 13.3.2017] {{comment|# |These should be presented somewhere|--[[User:Arja|Arja]] ([[User talk:Arja|talk]]) 07:39, 26 April 2017 (UTC)}}
* Original questionnaire analysis results [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QaMJZqUX0cPaTfOF 13.3.2017]  
 
* Consumption amount estimates
Consumption amount estimates
** Model run 21.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=xc0kaCs8cgzpjwo9] first distribution
* Model run 21.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=xc0kaCs8cgzpjwo9] first distribution
** Model run 18.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=YnuQMDJTQgW1Se5a with modelled data]; [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=KXXFiP0aj0DYEPdx with direct survey data]
* Model run 18.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=YnuQMDJTQgW1Se5a with modelled data]; [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=KXXFiP0aj0DYEPdx with direct survey data]
 
<rcode graphics=1 variables="
name:usesurvey|description:Do you want to use directly the survey data rather than modelled data?|type:selection|options:
FALSE;No, use modelled data;TRUE;Yes, use survey data
">
# This is code Op_en7749/ on page [[Goherr: Fish consumption study#Answer]]
 
library(OpasnetUtils)
library(ggplot2)
 
objects.latest("Op_en7749", code_name="initiate") # [[Goherr: Fish consumption study]] ovariables
DecisionTableParser(opbase.data("Op_en7748", subset = "Decisions"))
 
if(usesurvey) {
  objects.latest("Op_en7749", code_name="surveyjsp") # jsp ovariable directly based on survey data (N=2217)
  openv.setN(nrow(jsp@data))
}
 
amount <- EvalOutput(amount)
 
if(usesurvey) {
  oprint(summary(amount, marginals=c("Gender", "Country", "Fish","Ages","Foodpolicy")))
 
  print(ggplot(amount@output, aes(x=amountResult+0.1, colour=Country))+stat_ecdf()+scale_x_log10()+facet_wrap(~ Fish)+
    labs(x="Fish consumption (g /d)", y="Cumulative frequency")+theme_gray(base_size=24))
  print(ggplot(amount@output, aes(x=amountResult+0.1, colour=Ages))+stat_ecdf()+scale_x_log10()+facet_grid(Country ~ Fish)+
    labs(x="Fish consumption (g /d)", y="Cumulative frequency")+theme_gray(base_size=24))
  print(ggplot(amount@output, aes(x=amountResult+0.1, colour=Gender))+stat_ecdf()+scale_x_log10()+facet_grid(Country ~ Fish)+
    labs(x="Fish consumption (g /d)", y="Cumulative frequency")+theme_gray(base_size=24))
 
  print(ggplot(often@output, aes(x=oftenResult+0.1, colour=Country))+stat_ecdf()+scale_x_log10()+facet_wrap(~ Fish))
  print(ggplot(oftenside@output, aes(x=oftensideResult+0.1, colour=Country))+stat_ecdf()+scale_x_log10()+facet_wrap(~ Fish))
  print(ggplot(much@output, aes(x=muchResult+0.1, colour=Country))+stat_ecdf()+scale_x_log10()+facet_wrap(~ Fish))
  print(ggplot(muchside@output, aes(x=muchsideResult+0.1, colour=Country))+stat_ecdf()+scale_x_log10()+facet_wrap(~ Fish))
 
} else {
  oprint(summary(amount, marginals=c("Fish")))
 
  print(ggplot(amount@output, aes(x=amountResult+0.1, colour=Fish))+stat_ecdf()+scale_x_log10()+
    labs(x="Fish consumption (g /d)", y="Cumulative frequency")+theme_gray(base_size=24))
 
  print(ggplot(often@output, aes(x=oftenResult+0.1, colour=Fish))+stat_ecdf()+scale_x_log10())
  print(ggplot(oftenside@output, aes(x=oftensideResult+0.1, colour=Fish))+stat_ecdf()+scale_x_log10())
  print(ggplot(much@output, aes(x=muchResult+0.1, colour=Fish))+stat_ecdf()+scale_x_log10())
  print(ggplot(muchside@output, aes(x=muchsideResult+0.1, colour=Fish))+stat_ecdf()+scale_x_log10())
}
</rcode>


== Rationale ==
== Rationale ==
Line 287: Line 246:


This code is used to preprocess the original questionnaire data from the above .csv file and to store the data as a usable variable to Opasnet base. The code stores a data.frame named survey.
This code is used to preprocess the original questionnaire data from the above .csv file and to store the data as a usable variable to Opasnet base. The code stores a data.frame named survey.
* Model run 26.2.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=fgMgRmZztjBKmdXN]
* Model run 11.7.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=Fa8TMg6h8DbEmhjp]
* Model run 27.3.2019 with country codes DK EE FI SE [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=6UJy6JKNE3EyNHEQ]


<rcode name="preprocess2" label="Preprocess (only for developers)">
<rcode name="preprocess2" label="Preprocess (only for developers)">
Line 300: Line 260:
# Get the data either from Opasnet or your own hard drive.
# Get the data either from Opasnet or your own hard drive.


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


survey1 <- opasnet.csv(
survey1 <- opasnet.csv(
Line 306: Line 266:
   wiki = "opasnet_en", sep = ";", fill = TRUE, quote = "\""
   wiki = "opasnet_en", sep = ";", fill = TRUE, quote = "\""
)
)
#survey1 <- re#ad.csv(file = "N:/Ymal/Projects/Goherr/WP5/Goherr_fish_consumption.csv",
#survey1 <- re#ad.csv(fi_le = "N:/Ymal/Projects/Goherr/WP5/Goherr_fish_consumption.csv",
#                  header=FALSE, sep=";", fill = TRUE, quote="\"")
#                  header=FALSE, sep=";", fill = TRUE, quote="\"")


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


Line 329: Line 289:
   ifelse(as.numeric(as.character(survey1$Age)) < 46, "18-45",">45"),
   ifelse(as.numeric(as.character(survey1$Age)) < 46, "18-45",">45"),
   levels = c("18-45", ">45"), ordered = TRUE
   levels = c("18-45", ">45"), ordered = TRUE
)
survey1$Country <- factor(
  survey1$Country,
  levels=c("DK","EST","FI","SWE"),
  labels=c("DK","EE","FI","SE")
)
)
# Anonymize data
# Anonymize data
Line 450: Line 415:
* Sketches about modelling determinants of eating (spring 2018) [http://en.opasnet.org/en-opwiki/index.php?title=Benefit-risk_assessment_of_Baltic_herring_and_salmon_intake&oldid=41947#Sketches_about_modelling_determinants_of_eating]
* Sketches about modelling determinants of eating (spring 2018) [http://en.opasnet.org/en-opwiki/index.php?title=Benefit-risk_assessment_of_Baltic_herring_and_salmon_intake&oldid=41947#Sketches_about_modelling_determinants_of_eating]


==== Figures for the first manuscript ====
==== Figures, tables and stat analyses for the first manuscript ====


* Model run 15.6.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=eiVRACT5B9AUJMUF]
* Model run 8.5.2019 with thlVerse code (not run on Opasnet) [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=46qKuh08I7PHqXRX] [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=lnq3yNAPXAJFxwaG]
* Model run 16.6.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=RfEKhGyyVlMG6frc]
* Model run 26.8.2018 with fig 6 as table [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=B2jMOHfuUSmTjfrn]
* Previous model runs are [http://en.opasnet.org/en-opwiki/index.php?title=Goherr:_Fish_consumption_study&oldid=42828#Figures.2C_tables_and_stat_analyses_for_the_first_manuscript archived].
 
Methodological concerns:
* [https://stackoverflow.com/questions/12953045/warning-non-integer-successes-in-a-binomial-glm-survey-packages Is the warning in logistic regression important?]
* [https://stats.stackexchange.com/questions/46345/how-to-calculate-goodness-of-fit-in-glm-r Goodness of fit in logistic regression]
* [https://stats.stackexchange.com/questions/8661/logistic-regression-in-r-odds-ratio Calculate odds ratios in logistic regression]
* [https://www3.nd.edu/~rwilliam/stats3/ordinalindependent.pdf You might treat independent ordinal variables as continuous]
* [https://www.color-blindness.com/coblis-color-blindness-simulator/ Color blindness simulator to adjust colors for the color blind and for black and white printing]


<rcode graphics=1>
<rcode graphics=1>
Line 461: Line 434:
library(ggplot2)
library(ggplot2)
library(reshape2)
library(reshape2)
library(MASS)
#library(extrafont) # Needed to save Arial fonts to PDF or EPS
#library(thlVerse)
#library(car)
#library(car)
#library(vegan)
#library(vegan)


BS <- 24
BS <- 24
localcomp <- FALSE
thl <- FALSE
# Set colours
#library(thlVerse)
#library(tidyverse)
#colors <- c(
#  thlColors(6,"quali","bar",thin=0.8),
#  thlColors(6,"quali","bar",thin=1)
#)[c(1,8,3,10,11,6)]
#ggplot(tibble(A=as.character(1:6),B=1),aes(x=A,weight=B,fill=A))+
#  geom_bar()+
#  scale_fill_manual(values=colors)
colors <- c("#74AF59FF","#2F62ADFF","#D692BDFF","#29A0C1FF","#BE3F72FF","#FBB848FF")
########################## Functions
#### thlLinePlot was adjusted to enable different point shapes
thlLinePlot <- function (data, xvar, yvar, groupvar = NULL, ylabel = yvar,
                        xlabel = NULL, colors = thlColors(n = 12, type = "quali",
                                                          name = "line"), title = NULL, subtitle = NULL, caption = NULL,
                        legend.position = "none", base.size = 16, linewidth = 3,
                        show.grid.x = FALSE, show.grid.y = TRUE, lang = "fi", ylimits = NULL,
                        marked.treshold = 10, plot.missing = FALSE, xaxis.breaks = waiver(),
                        yaxis.breaks = waiver(), panels = FALSE, nrow.panels = 1,
                        labels.end = FALSE)
{
  lwd <- thlPtsConvert(linewidth)
  gg <- ggplot(data, aes_(x = substitute(xvar), y = substitute(yvar),
                          group = ifelse(!is.null(substitute(groupvar)), substitute(groupvar),
                                        NA), colour = ifelse(!is.null(substitute(groupvar)),
                                                              substitute(groupvar), ""))) + geom_line(size = lwd)
  if (isTRUE(plot.missing)) {
    df <- thlNaLines(data = data, xvar = deparse(substitute(xvar)),
                    yvar = deparse(substitute(yvar)), groupvar = unlist(ifelse(deparse(substitute(groupvar)) !=
                                                                                  "NULL", deparse(substitute(groupvar)), list(NULL))))
    if (!is.null(df)) {
      gg <- gg + geom_line(data = df, aes_(x = substitute(xvar),
                                          y = substitute(yvar), group = ifelse(!is.null(substitute(groupvar)),
                                                                                substitute(groupvar), NA), colour = ifelse(!is.null(substitute(groupvar)),
                                                                                                                          substitute(groupvar), "")), linetype = 2,
                          size = lwd)
    }
  }
  if (!is.null(marked.treshold)) {
    if (length(unique(data[, deparse(substitute(xvar))])) >
        marked.treshold) {
      if (is.factor(data[, deparse(substitute(xvar))]) ||
          is.character(data[, deparse(substitute(xvar))]) ||
          is.logical(data[, deparse(substitute(xvar))])) {
        levs <- levels(factor(data[, deparse(substitute(xvar))]))
        min <- levs[1]
        max <- levs[length(levs)]
      }
      else {
        min <- min(data[, deparse(substitute(xvar))])
        max <- max(data[, deparse(substitute(xvar))])
      }
      subdata <- data[c(data[, deparse(substitute(xvar))] %in%
                          c(min, max)), ]
      gg <- gg + geom_point(data = subdata, aes_(
        x = substitute(xvar),
        y = substitute(yvar),
        group = ifelse(!is.null(substitute(groupvar)), substitute(groupvar), NA),
        colour = ifelse(!is.null(substitute(groupvar)), substitute(groupvar), ""),
        shape = substitute(groupvar)),
        fill = "white",
        stroke = 1.35 * lwd,
        size = 10/3 * lwd)+scale_shape_manual(values=21:25)
    }
    else {
      gg <- gg + geom_point(
        aes_(
          shape = substitute(groupvar)
        ),
        fill = "white",
        stroke = 1.35 * lwd,
        size = 10/3 * lwd)+scale_shape_manual(values=21:25)
    }
  }
  if (isTRUE(labels.end)) {
    if (is.factor(data[, deparse(substitute(xvar))]) ||
        is.character(data[, deparse(substitute(xvar))]) ||
        is.logical(data[, deparse(substitute(xvar))])) {
      levs <- levels(factor(data[, deparse(substitute(xvar))]))
      maxd <- data[data[, deparse(substitute(xvar))] ==
                    levs[length(levs)], ]
    }
    else {
      maxd <- data[data[, deparse(substitute(xvar))] ==
                    max(data[, deparse(substitute(xvar))]), ]
    }
    brks <- maxd[, deparse(substitute(yvar))]
    labsut <- maxd[, deparse(substitute(groupvar))]
  }
  else (brks <- labsut <- waiver())
  gg <- gg + ylab(ifelse(deparse(substitute(ylabel)) == "yvar",
                        deparse(substitute(yvar)), ylabel)) + labs(title = title,
                                                                    subtitle = subtitle, caption = caption) + thlTheme(show.grid.y = show.grid.y,
                                                                                                                      show.grid.x = show.grid.x, base.size = base.size, legend.position = legend.position,
                                                                                                                      x.axis.title = ifelse(!is.null(xlabel), TRUE, FALSE)) +
    xlab(ifelse(!is.null(xlabel), xlabel, "")) + scale_color_manual(values = colors) +
    thlYaxisControl(lang = lang, limits = ylimits, breaks = yaxis.breaks,
                    sec.axis = labels.end, sec.axis.breaks = brks, sec.axis.labels = labsut)
  if (is.factor(data[, deparse(substitute(xvar))]) || is.character(data[,
                                                                        deparse(substitute(xvar))]) || is.logical(data[, deparse(substitute(xvar))])) {
    gg <- gg + scale_x_discrete(breaks = xaxis.breaks, expand = expand_scale(mult = c(0.05)))
  }
  else (gg <- gg + scale_x_continuous(breaks = xaxis.breaks))
  if (isTRUE(panels)) {
    fmla <- as.formula(paste0("~", substitute(groupvar)))
    gg <- gg + facet_wrap(fmla, scales = "free", nrow = nrow.panels)
  }
  gg
}


groups <- function(o) {
groups <- function(o) {
Line 476: Line 570:
   return(o)
   return(o)
}
}
#' @title reasons produces a ggplot object for a graph that shows fractions of reasons mentioned in the data.
#' @param dat a data.frame containing survey answers for the target subgroup.
#' @param title an atomic string with the name to be shown as the title of graph
#' @return ggplot object


reasons <- function(
reasons <- function(
   dat, # data.frame to be used
   dat,
   title # title on graph
   title
) {
) {
  require(reshape2)
  weight <- Ovariable("weight",data=data.frame(
    Row=dat$Row,
    Result=dat$Weight
  ))
  levels(dat$Country) <- paste0(
    levels(dat$Country), " (n=",
    aggregate(dat$Country, dat["Country"], length)$x, ")"
  )
   tmp <- melt(
   tmp <- melt(
     dat,
     dat[setdiff(colnames(dat),"Weighting")],
     id.vars = "Country",
     id.vars = c("Country","Row"),
     variable.name ="Reason",
     variable.name ="Reason",
     value.name="Result"
     value.name="Result"
   )
   )
   tmp$Result <- (tmp$Result %in% c("Yes","1"))*1
   tmp$Result <- (tmp$Result %in% c("Yes","1"))*1
   tmp <- EvalOutput(Ovariable("tmp",data=tmp))
   tmp <- weight * (Ovariable("tmp",data=tmp))
   tmp <- (oapply(tmp, c("Country","Reason"), sum)/oapply(tmp, "Country",sum))@output
  cat(title, ": number of individual answers.\n")
  oprint(oapply(tmp, c("tmpResult","Country"),length)@output)
 
   tmp <- (oapply(tmp, c("Country","Reason"), sum)/oapply(tmp, "Country",sum))
  popu <- oapply(tmp, "Reason",sum)@output
  popu <- popu$Reason[order(popu$Result)]
  tmp$Reason <- factor(tmp$Reason, levels=popu)
  levels(tmp$Reason) <- gsub("( BH| BS)", "", gsub("\\.", " ", levels(tmp$Reason)))
  cat(title, ": fractions shown on graph.\n")
  oprint(tmp@output)
 
  if(thl) {
    tmp <- thlLinePlot(tmp@output, xvar=Reason, yvar=Result,groupvar=Country,
                      colors= c("#519B2FFF", "#2F62ADFF", "#BE3F72FF","#88D0E6FF"), # #29A0C1FF"),
                      # THL colors but fourth is brigter
                      legend.position = c(0.85,0.2), base.size = BS, title=title,
                      subtitle="Fraction of population")+
      coord_flip()+
      scale_y_continuous(labels=scales::percent_format(accuracy=1))
  } else {
    tmp <- ggplot(tmp@output, aes(x=Reason, y=Result,colour=Country, group=Country))+
      geom_point(shape=21, size=5, fill="Grey", stroke=2)+
      geom_line(size=1.2)+
      coord_flip()+
      theme_gray(base_size=BS)+
      scale_y_continuous(labels=scales::percent_format())+#accuracy=1))+
      scale_colour_manual(values=colors)+
      labs(
        title=title,
        x="Answer",
        y="Fraction of population")
  }
  return(tmp)
}
 
#' @title Function impacts produces a graph showing the magnitude of impact if some things change
#' @param dat data.frame containing answers
#' @param title string to be used as the title of the graph
#' @param population ovariable of country population weights
#' @return ggplot object
 
impacts <- function(dat,title, population) {
  tmp <- melt(
    dat,
    id.vars = c("Country","Weighting"),
    variable.name ="Reason",
    value.name="Response"
  )
  tmp$Response <- factor(
    tmp$Response,
    levels=c(
      "I don't know",
      "Increase it over to double",
      "Increase it over by half",
      "Increase it a bit",
      "No effect",
      "Decrease it a bit",
      "Decrease it to less than half",
      "Decrease it to zero"
    ),
    ordered=TRUE
  )
  levels(tmp$Response)[levels(tmp$Response) %in% c(
    "Decrease it to zero",
    "Decrease it to less than half",
    "Decrease it a bit"
  )] <- "Decrease"
  levels(tmp$Response)[levels(tmp$Response) %in% c(
    "Increase it a bit",
    "Increase it over by half",
    "Increase it over to double"
  )] <- "Increase"
 
  colnames(tmp)[colnames(tmp)=="Weighting"] <- "Result"
  tmp <- Ovariable("tmp",data=tmp) * population
 
  popu <- aggregate(as.numeric(tmp$Response), tmp@output["Reason"],sum)
  popu <- popu$Reason[order(popu$x)]
  tmp$Reason <- factor(tmp$Reason, levels=popu)
  levels(tmp$Reason) <- gsub("( BH| BS)", "", gsub("\\.", " ", levels(tmp$Reason)))
    
    
   tmp <- ggplot(tmp, aes(x=Reason, weight=Result,fill=Country))+
   tmp <- (oapply(tmp, c("Reason","Response"), sum)/
    geom_bar(position="dodge")+
            oapply(tmp, c("Reason"),sum))
    coord_flip()+
 
    theme_gray(base_size=BS)+
  oprint(tmp@output)
    scale_y_continuous(labels=scales::percent_format())+
 
    labs(title=title)
  print(ggplot(tmp@output,
              aes(x=Reason, weight=Result, fill=Response))+geom_bar()+
          coord_flip()+#facet_grid(.~Country)+
          theme_gray(base_size=BS)+
          theme(legend.position = "bottom")+
          scale_fill_manual(values=c("Grey",colors[c(1,6,5)]))+
          scale_y_continuous(breaks=c(0.5,1),labels=scales::percent_format())+
          labs(
            title=title,
            x="Cause",
            y="Fraction of population"
          ))
 
  return(tmp@output)
}
 
#' @title or95 prints the odds ratio and 95 % confidence interval of a glm fit
#' @param fit a glm fit
#' @return data.frame with OR and 95%CI
 
or95 <- function(fit) {
  cat("Odds ratios with 95 % confidence intervals\n")
  tmp <- exp(cbind(coef(fit), confint(fit)))
  tmp <- data.frame(
    rownames(tmp),
    paste0(sprintf("%.2f",tmp[,1]), " (",sprintf("%.2f", tmp[,2]),"-",sprintf("%.2f",tmp[,3]),")")
  )
  colnames(tmp) <- c("Parameter","OR (95% CI)")
   return(tmp)
   return(tmp)
}
}


objects.latest("Op_en7749", "preprocess2") # [[Goherr: Fish consumption study]]: survey, surv
#' @title Perform and print regression analyses
#' @param dat data.frame with data
#' @param y name of column to be used as dependent variable
#' @param countries logical atom to tell whether to run country-specific analyses
#' @return returns the model fit object
 
printregr <- function(dat, y, title, countries = TRUE) {
 
  dat$y <- dat[[y]]
  cat("\n############### Logistic regression analysis:",title, "(", y, ")\n")
  fit <- glm(
    y ~ Ages + Gender + Country + as.numeric(Education) + as.numeric(Purchasing.power),
    family=binomial(),na.action="na.omit",weights=Weighting,
    data=dat
  )
  step <- stepAIC(fit, direction = "both")
  oprint(step$anova)
  oprint(summary(fit))
  oprint(or95(fit))
 
  if(countries) {
    for(i in unique(dat$Country)) {
      cat("\n#### Country-specific logistic regression analysis:", i, y, "\n")
      fit <- glm(
        Eat.fish ~ Ages + Gender + as.numeric(Education) + as.numeric(Purchasing.power),
        family=binomial(),na.action="na.omit",weights=Weighting,
        data=dat[dat$Country==i,]
      )
      step <- stepAIC(fit, direction = "both")
      oprint(step$anova)
      oprint(summary(fit))
      oprint(or95(fit))
    }
  }
  return(fit)
}
 
#################### Data for Figure 1.
 
if(thl) { # Commented out because needs package thlVerse.
 
  dat <- opbase.data("Op_en7749", subset="Fish consumption as food in Finland") # [[Goherr: Fish consumption study]]
  dat$Year <- as.numeric(as.character(dat$Year))
 
  gr <- thlLinePlot(
    dat[dat$Species=="Total",],
    xvar=Year,
    yvar=Result,
    ylimits=c(0,12),
    groupvar=Origin,
    base.size=24,
    title="Fish consumption in Finland",
    ylabel="",
    subtitle="(kg/a per person)",
    legend.position = "bottom"
  )
 
  if(localcomp) ggsave("Figure 1.pdf", width=8,height=6)
  if(localcomp) ggsave("Figure 1.png", width=8,height=6)
 
} # End if
################## Get data
 
objects.latest("Op_en7749", "preprocess2") # [[Goherr: Fish consumption study]]: survey1, surv
objects.latest("Op_en7749", "nonsamplejsp") # [[Goherr: Fish consumption study]]: jsp every respondent from data exactly once.
objects.latest("Op_en7749", "initiate") # [[Goherr: Fish consumption study]]: amount etc.
openv.setN(50)
survey1 <- groups(survey1)
#levels(survey1$Country) <- c("FI","SE","DK","EE") # WHY was this here? It is a potential hazard.
effinfo <- 0 # No policies implemented.
effrecomm <- 0
 
population <- EvalOutput(Ovariable(
  "population",
  ddata = "Op_en7748",
  subset = "Population"
))
levels(population$Country) <- c("DK","EE","FI","SE")
population <- oapply(population, c("Country"),sum)
 
amount <- EvalOutput(amount)
#result(amount)[result(amount)==0] <- 0.1 Not needed if not log-transformed.
amount <- amount * 365.25 / 1000 # g/day --> kg/year


#### General fish eating
#### General fish eating
Line 512: Line 808:
))
))


oapply(su1, "Country", sum)@output
# Figure 1. Origin of consumed fish in Finland between 1999 and 2016.
tmp <- (oapply(su1, c("Country","Gender"), sum)/oapply(su1,"Country",sum))@output
# Data not on this page, drawn separately.
cat("Fraction of genders\n")
 
oprint(tmp[order(tmp$Country),])
# Table 1. Dimensions of embeddedness, modified from Hass (2007, p. 16)
tmp <- (oapply(su1, c("Country","Education"), sum)/oapply(su1,"Country",sum))@output
# Written directly on the manuscript.
cat("Fraction of education levels\n")
 
oprint(tmp[order(tmp$Country),])
# Table 2: statistics of the survey population in each country (n, female %, education, purchasing power)
tmp <- (oapply(su1, c("Country","Purchasing.power"), sum)/oapply(su1,"Country",sum))@output
 
cat("Fraction of purchasing power\n")
tmp <- round(t(data.frame(
oprint(tmp[order(tmp$Country),])
  n = tapply(rep(1,nrow(survey1)), survey1$Country,sum),
  Females = tapply(survey1$Gender=="Female", survey1$Country,mean)*100,
  Old = tapply(survey1$Ages==">45", survey1$Country,mean)*100,
  PurcVlow = tapply(survey1$Purchasing.power=="Very low", survey1$Country,mean)*100,
  PurcLow = tapply(survey1$Purchasing.power=="Low", survey1$Country,mean)*100,
  PurcSuf = tapply(survey1$Purchasing.power=="Sufficient", survey1$Country,mean)*100,
  PurcGood = tapply(survey1$Purchasing.power=="Good", survey1$Country,mean)*100,
  PurcVgood = tapply(survey1$Purchasing.power=="Very good", survey1$Country,mean)*100,
  PurcExc = tapply(survey1$Purchasing.power=="Excellent", survey1$Country,mean)*100,
  EducPri = tapply(survey1$Education=="Primary education", survey1$Country,mean)*100,
  EducSec = tapply(survey1$Education=="Secondary education (gymnasium, vocational school or similar)", survey1$Country,mean)*100,
  EducCol = tapply(survey1$Education=="Lower level college education or similar", survey1$Country,mean)*100,
  EducHig = tapply(survey1$Education=="Higher level college education or similar", survey1$Country,mean)*100
)))
rownames(tmp) <- c(
  "Number of respondents",
  "Females (%)",
  ">45 years (%)",
  "Purchasing power: Very low (%)",
  "Low (%)",
  "Sufficient (%)",
  "Good (%)",
  "Very good (%)",
  "Excellent (%)",
  "Education: Primary education (%)",
  "Secondary education (%)",
  "Lower level college (%)",
  "Higher level college (%)"
)
oprint(tmp, digits=0)
 
cat("Weighted fraction of fish eaters\n")
tmp <- aggregate(survey1$Weighting, survey1[c("Eat.fish","Country")],FUN=sum)
oprint(aggregate(tmp$x, tmp["Country"],FUN=function(x) x[2]/sum(x)))


ggplot(survey1,  
ggplot(survey1,  
       aes(x=Eat.fish, group=Country), stat="count")+
       aes(x=Eat.fish, weight=Weighting, group=Country), stat="count")+
   geom_bar(aes(y=..prop.., fill=Country), position="dodge")+
   geom_bar(aes(y=..prop.., fill=Country), position="dodge")+
   theme_gray(base_size=BS)+
   theme_gray(base_size=BS)+
  scale_fill_manual(values=colors)+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   labs(title="Fraction of fish eaters within the study population")
   labs(title="Fraction of fish eaters within population")
 
############## Answers to open questions
if(FALSE){
  cat("Number of open-ended answers per question and country\n")
 
  tmp <- survey1[c(1,3,27:28,36:37,44:45,60:61,71:72,83:84,93:94,109:110,121:122,133:134,154,157:158)]
  tmp <- reshape(
    tmp,
    varying=list(
      Condition = colnames(tmp)[grepl("other",tolower(colnames(tmp)))],
      Text = colnames(tmp)[grepl("open",colnames(tmp))]
    ),
    times = colnames(tmp)[grepl("other",tolower(colnames(tmp)))],
    direction="long"
  )
  tmp <- tmp[tmp$Why.not.fish.open!="" , ]
  oprint(table(tmp$time,tmp$Country))
  oprint(tmp)
} # End if(FALSE)


#### Baltic salmon
#### Baltic salmon


reasons(
reasons(
   survey1[survey1$Eat.fish=="No"& !is.na(survey1$Eat.fish),c(1,17:27)], # 1 Country, 8 Why don't you eat fish
   survey1[survey1$Eat.fish=="No"& !is.na(survey1$Eat.fish),c(1,17:26,154,157)], # 1 Country, 8 Why don't you eat fish
   "Reason among non-fish-consumers"
   "Reasons not to eat among fish non-consumers"
)
)
# Figure 5
cat("Figure 5 (to be converted to text)\n")
reasons(
reasons(
   survey1[survey1$Eat.salmon=="Yes"& !is.na(survey1$Eat.salmon),c(1,31:36)], # 1 Country, 11 Which salmon species
   survey1[survey1$Eat.salmon=="Yes"& !is.na(survey1$Eat.salmon),c(1,31:36,154,157)], # 1 Country, 11 Which salmon species
   "Species among salmon consumers"
   "Species used among salmon consumers"
)
)
reasons(
reasons(
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,50:60)], # 1 Country, 17 Why Baltic salmon
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,50:59,154,157)], # 1 Country, 17 Why Baltic salmon
   "Reasons among Baltic salmon consumers"
   "Reasons to eat among Baltic salmon consumers"
)
)
if(localcomp) ggsave("Figure 5.pdf", width=10,height=5)
if(localcomp) ggsave("Figure 5.png", width=10,height=5)
reasons(
reasons(
   survey1[survey1$Baltic.salmon=="No"& !is.na(survey1$Baltic.salmon),c(1,62:71)], # 1 Country, 18 Why not Baltic salmon
   survey1[survey1$Baltic.salmon=="No"& !is.na(survey1$Baltic.salmon),c(1,62:70,154,157)], # 1 Country, 18 Why not Baltic salmon
   "Reasons among Baltic salmon non-consumers"
   "Reasons not to eat, Baltic salmon non-consumers"
)
)


tmp <- melt(
if(localcomp) ggsave("Figure 6.pdf", width=10,height=5)
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,73:83)],
if(localcomp) ggsave("Figure 6.png", width=10,height=5)
   id.vars = "Country",
 
   variable.name ="Reason",
impacts.sal <- impacts(
  value.name="Response"
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,73:82,154)],
   "Effect on Baltic salmon consumption",
   population
)
)
tmp$Response <- factor(
 
   tmp$Response,
####### Any herring
   levels=c(
 
    "Decrease it to zero",
ggplot(survey1[survey1$Eat.fish=="Yes",],
     "Decrease it to less than half",
      aes(x=Country, weight=Weighting, fill=Eat.herring))+geom_bar(position="fill")+
     "Decrease it a bit",
   coord_flip()+#facet_grid(.~Country)+
    "No effect",
   theme_gray(base_size=BS)+
    "Increase it a bit",
  theme(legend.position = "bottom")+
    "Increase it over by half",
  scale_y_continuous(labels=scales::percent_format())+
    "Increase it over to double",
  scale_fill_manual(values=colors)+
    "I don't know"
  labs(
  ),
     title="Any herring consumption",
   ordered=TRUE
     y="Fraction of fish consumers"
  )
 
####### Baltic herring
 
### Figure 2. Comparison of Baltic herring consumption habits.
 
survey1$What <- ifelse(is.na(survey1$Eat.BH), "Other fish", as.character(survey1$Eat.BH))
survey1$What <- factor(
  survey1$What,
  levels=(c("Yes","I don't know","No","Other fish")),
   labels=(c("Baltic herring","Some herring","Other herring","Other fish"))
)
)


tmp$Result <- 1
tmp <- survey1[survey1$Eat.fish=="Yes",]
tmp <- EvalOutput(Ovariable("tmp",data=tmp))
colnames(tmp)[colnames(tmp)=="Weighting"] <- "Result"
tmp <- (oapply(tmp, c("Country","Reason","Response"), sum)/
tmp <- EvalOutput(Ovariable("tmp", data=tmp))
          oapply(tmp, c("Country","Reason"),sum))@output
tmp <- oapply(tmp, c("What","Country"), sum) / oapply(tmp,"Country",sum)
tmp$Country <- factor(tmp$Country, levels=rev(levels(tmp$Country)))
 
if(thl) {
  thlBarPlot(tmp@output, xvar=Country, yvar=Result, groupvar=What, horizontal = TRUE, stacked = TRUE,
            legend.position = "bottom",,
            title="Baltic herring consumption",
            subtitle="Fraction of fish consumers",
  )+
    scale_y_continuous(labels=scales::percent_format())
 
} else {
  ggplot(tmp@output, aes(x=Country, weight=Result, fill=What))+geom_bar(position="fill")+
    coord_flip()+#facet_grid(.~Country)+
    theme_gray(base_size=BS)+
    theme(legend.position = "bottom")+
    scale_y_continuous(labels=scales::percent_format())+
    scale_fill_manual(values=c("gray",colors))+
    guides(fill=guide_legend(reverse=TRUE))+
    labs(
      title="Baltic herring consumption",
      y="Fraction of fish consumers"
    )
 
}
 
if(localcomp) ggsave("Figure 2.pdf",width=10, height=5)
if(localcomp) ggsave("Figure 2.png",width=6, height=3)
 
tmp <- survey1[survey1$Eat.fish=="Yes",]
tmp <- aggregate(tmp$Weighting, tmp[c("What","Country")], sum)
 
rown <- tmp$What[tmp$Country=="FI"]
tmp <- aggregate(tmp$x, tmp["Country"],FUN=function(x) x/sum(x))
colnames(tmp[[2]]) <- rown
oprint(as.matrix(tmp))
 


ggplot(tmp,
ggplot(survey1[survey1$Eat.herring=="Yes",],
       aes(x=Reason, weight=Result, fill=Response))+geom_bar()+
       aes(x=Country, weight=Weighting, fill=Eat.BH))+geom_bar(position="fill")+
   coord_flip()+facet_grid(.~Country)+
   coord_flip()+#facet_grid(.~Country)+
   theme_gray(base_size=BS)+
   theme_gray(base_size=BS)+
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   labs(title="Effect on Baltic salmon consumption")
  scale_y_continuous(labels=scales::percent_format())+
  scale_fill_manual(values=c(colors))+
   labs(
    title="Baltic herring consumption",
    y="Fraction of herring consumers"
  )


####### Baltic herring
####### Any salmon
 
ggplot(survey1[survey1$Eat.fish=="Yes",],
      aes(x=Country, weight=Weighting, fill=Eat.salmon))+geom_bar(position="fill")+
  coord_flip()+
  theme_gray(base_size=BS)+
  theme(legend.position = "bottom")+
  scale_y_continuous(labels=scales::percent_format())+
  scale_fill_manual(values=colors)+
  labs(
    title="Any salmon consumption",
    y="Fraction of fish consumers"
  )
 
####### Baltic salmon
 
ggplot(survey1[survey1$Eat.fish=="Yes",],
      aes(x=Country, weight=Weighting, fill=Baltic.salmon))+geom_bar(position="fill")+
  coord_flip()+
  theme_gray(base_size=BS)+
  theme(legend.position = "bottom")+
  scale_y_continuous(labels=scales::percent_format())+
  scale_fill_manual(values=colors)+
  labs(
    title="Baltic salmon consumption",
    y="Fraction of fish consumers"
  )
 
ggplot(survey1[survey1$Eat.salmon=="Yes",],
      aes(x=Country, weight=Weighting, fill=Baltic.salmon))+geom_bar(position="fill")+
  coord_flip()+
  theme_gray(base_size=BS)+
  theme(legend.position = "bottom")+
  scale_y_continuous(labels=scales::percent_format())+
  scale_fill_manual(values=colors)+
  labs(
    title="Baltic salmon consumption",
    y="Fraction of salmon consumers"
  )
 
############## Recommendation awareness
 
levels(survey1$Recommendation.awareness) <- c("Not aware","Aware","Know content")
ggplot(survey1,
      aes(x=Country, weight=Weighting, fill=Recommendation.awareness))+geom_bar(position="fill")+
  coord_flip()+#facet_grid(.~Country)+
  theme_gray(base_size=BS)+
  theme(legend.position = "bottom")+
  scale_y_continuous(labels=scales::percent_format())+
  scale_fill_manual(values=colors)+
  labs(title="Awareness of food recommendations about Baltic fish")
 
# Figure 3. Percentages of reasons to eat Baltic herring


reasons(
reasons(
   survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,99:109)], # 1 Country, 27 why Baltic herring
   survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,99:108,154,157)], # 1 Country, 27 why Baltic herring
   "Reasons among Baltic herring consumers"
   "Reasons to eat among Baltic herring consumers"
)
)
if(localcomp) ggsave("Figure 3.pdf", width=10, height=5)
if(localcomp) ggsave("Figure 3.png", width=10, height=5)
# Figure 4. Reasons for not to eat Baltic herring.
# Denmark and Estonia were omitted because they had less than 20 observations.
reasons(
reasons(
   survey1[survey1$Eat.BH=="No"& !is.na(survey1$Eat.BH),c(1,111:121)], # 1 Country, 28 Why not Baltic herring
   survey1[
   "Reasons among Baltic herring non-consumers"
    survey1$Eat.BH=="No"&
      survey1$Country %in% c("FI","SE") &
      !is.na(survey1$Eat.BH),c(1,111:120,154,157)], # 1 Country, 28 Why not Baltic herring
   "Reasons not to eat, Baltic herring non-consumers"
)+scale_color_manual(values=c("#BE3F72FF","#88D0E6FF"))+
  scale_shape_manual(values=23:24)
 
if(localcomp) ggsave("Figure 4.pdf", width=10, height=5)
if(localcomp) ggsave("Figure 4.png", width=10, height=5)
 
# Figure 7. The role of different determinants on Baltic salmon consumption
 
impacts.herr <- impacts(
  survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,123:132,154)], # 1 Country, 29 Influence of Baltic herring policies
  title="Effect on Baltic herring consumption",
  population
)
)


tmp <- melt(
tmp <- rbind(
   survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,123:133)], # 1 Country, 29 Influence of Baltic herring policies
   data.frame(
  id.vars = "Country",
    Fish="Herring",
  variable.name ="Reason",
    impacts.herr
  value.name="Response"
)
tmp$Response <- factor(
  tmp$Response,
  levels=c(
    "Decrease it to zero",
    "Decrease it to less than half",
    "Decrease it a bit",
    "No effect",
    "Increase it a bit",
    "Increase it over by half",
    "Increase it over to double",
    "I don't know"
   ),
   ),
   ordered=TRUE
   data.frame(
    Fish="Salmon",
    impacts.sal
  )
)
)


tmp$Result <- 1
oprint(aggregate(tmp$Result,tmp[c("Reason","Response","Fish")],sum))
tmp <- EvalOutput(Ovariable("tmp",data=tmp))
 
tmp <- (oapply(tmp, c("Country","Reason","Response"), sum)/
if(thl) {
          oapply(tmp, c("Country","Reason"),sum))@output
  tmp$Response <- factor(tmp$Response, levels=rev(levels(tmp$Response)))
  thlBarPlot(aggregate(tmp$Result,tmp[c("Reason","Response","Fish")],sum),
            xvar=Reason, yvar=x, groupvar=Response, legend.position = "bottom",
            horizontal = TRUE, stacked = TRUE,
            base.size = BS, title="Effect on Baltic fish consumption",
            subtitle="Fraction of population")+
    facet_grid(.~Fish)+
    scale_fill_manual(values=rev(c("gray",colors[c(1,6,5)])))+
    scale_y_continuous(breaks=c(0.5,1),labels=scales::percent_format(accuracy=1))
 
} else {
  ggplot(tmp, aes(x=Reason, weight=Result, fill=Response))+geom_bar()+
    coord_flip()+facet_grid(.~Fish)+
    theme_gray(base_size=BS)+
    theme(legend.position = "bottom")+
    scale_fill_manual(values=c("gray",colors[c(1,6,5)]))+
    scale_y_continuous(breaks=c(0.5,1),labels=scales::percent_format())+#accuracy=1))+
    guides(fill=guide_legend(reverse=TRUE))+
    labs(
      title="Effect on Baltic fish consumption",
      x="Cause",
      y="Fraction of population"
    )
}
 
if(localcomp) ggsave("Figure 7.pdf", width=12, height=7)
if(localcomp) ggsave("Figure 7.png", width=12, height=7)
 
######### Amounts estimated for each respondent


ggplot(tmp,
amount <- groups(amount)
      aes(x=Reason, weight=Result, fill=Response))+geom_bar()+
 
   coord_flip()+facet_grid(.~Country)+
#tmp <- amount*info
#tmp <- tmp[
#  tmp$Info.improvements=="BAU" &
#    tmp$Recomm.herring=="BAU" &
#    tmp$Recomm.salmon=="BAU" ,
#  ]
tmp <- cbind(
  oapply(amount, c("Country","Group","Gender","Ages","Fish"), FUN=mean)@output,
  SD=oapply(amount, c("Country","Group","Gender","Ages","Fish"), FUN=sd)$amountResult
)
 
# Figure 3 (old, unused). Average consumption (kg/year) of Baltic herring in four countries,
# calculated with Monte Carlo simulation (1000 iterations).
 
ggplot(tmp[tmp$Fish=="Herring",], aes(x=Group, weight=amountResult, fill=Group))+
   geom_bar()+facet_wrap(~Country)+
   theme_gray(base_size=BS)+
   theme_gray(base_size=BS)+
   theme_gray(base_size=BS)+
   scale_fill_manual(values=colors)+
   theme(legend.position = "bottom")+
  guides(fill=FALSE)+
   labs(title="Effect on Baltic salmon consumption")
   theme(axis.text.x = element_text(angle = -90))+
   labs(
    title="Baltic herring consumption in subgroups",
    y="Average consumption (kg/year)")
 
# if(localcomp) ggsave("Figure3.pdf", width=9, height=10)
# if(localcomp) ggsave("Figure3.png", width=9, height=10)
 
tmp$Result <- paste0(sprintf("%.1f",tmp$amountResult), " (",sprintf("%.1f", tmp$SD,1),")")
tmp <- reshape(tmp, v.names="Result", timevar="Country", idvar=c("Fish","Group"),drop=c("Gender","Ages","amountResult","SD"), direction="wide")
colnames(tmp) <- gsub("Result\\.", "", colnames(tmp))
cat("Average fish consumption in subgroups, mean (sd)\n")
oprint(tmp)


objects.latest("Op_en7748", code_name="hia")
weight <- EvalOutput(Ovariable("weight",data=data.frame(
  amount@output[c("Row","Iter","Fish","Country")],
  Result=amount$Weighting
)))


tmp <- amount*info
tmp <- (oapply(amount * weight, c("Fish","Country"), sum) / oapply(weight,c("Fish","Country"), sum))@output
tmp <- tmp[
  tmp$Info.improvements=="BAU" &
    tmp$Recomm.herring=="BAU" &
    tmp$Recomm.salmon=="BAU" ,
  ]
tmp <- groups(oapply(tmp, c("Country","Gender","Ages","Fish"), FUN=mean)@output)
ggplot(tmp, aes(x=Group, weight=Result, fill=Fish))+
  geom_bar()+facet_grid(Country~Fish)+
  theme_gray(base_size=BS)+
  labs(title="Fish consumption in subgroups")


cat("Average fish consumption in subgroups\n")
cat("Average fish consumption per country (kg/year)\n")
oprint(tmp)
oprint(tmp)
result(amount)[result(amount)==0] <- 0.1
ggplot(amount@output, aes(x=amountResult, colour=Group))+stat_ecdf()+scale_x_log10()+
  facet_wrap(~Fish)
############################
#### Statistical analyses
dat <- merge(survey1,amount@output) # [amount$Fish=="Herring",]) # This results in 50 * 2 times of rows
# These should be corrected to preprocess2. NOT ordered.
dat$Country <- factor(dat$Country, ordered=FALSE)
dat$Ages <- factor(dat$Ages, ordered=FALSE)
dat$Gender <- factor(dat$Gender, ordered=FALSE)
ggplot(dat, aes(x=amountResult, weight=Weighting, color=Country))+stat_ecdf()+facet_wrap(~Fish)+
  scale_x_log10()
gro <- unique(dat$Group)
out <- list(data.frame(), data.frame())
for(i in unique(dat$Fish)) {
  for(j in unique(dat$Country)){
    for(k in unique(dat$Iter)) {
      res <- kruskal.test(amountResult ~ Group, data=dat[
        dat$Iter==k  & dat$Country==j & dat$Fish==i,])
      out[[1]] <- rbind(out[[1]], data.frame(
        test="Kruskal-Wallis",
        Fish=i,
        Country=j,
        Iter=k,
        p.value=res[[3]])
      )
      for(l in 1:(length(gro)-1)) {
        for(m in (l+1):length(gro)) {
          res2 <- wilcox.test(
            x=dat$amountResult[dat$Iter==k  & dat$Country==j & dat$Fish==i & dat$Group==gro[l]],
            y=dat$amountResult[dat$Iter==k  & dat$Country==j & dat$Fish==i & dat$Group==gro[m]], conf.int = FALSE)
          out[[2]] <- rbind(out[[2]], data.frame(
            Fish=i,
            Country=j,
            Iter=k,
            Pair1=gro[l],
            Pair2=gro[m],
            p.value = res2$p.value
          ))
        }
      }
    }
  }
}
ggplot(out[[1]], aes(x=paste(Fish, Country), y = p.value, colour = Country))+geom_jitter()+
  geom_hline(yintercept=0.05)+coord_flip()+theme_gray(base_size=BS)+
  labs(title="Kruskal-Wallis non-parametric test across all groups")
ggplot(out[[2]], aes(x=paste(Fish, Country), y=p.value, colour=Pair1, shape=Pair2))+geom_jitter()+
  geom_hline(yintercept=0.05)+coord_flip()+theme_gray(base_size=BS)+
  labs(title="Mann-Whitney U test across all pairs of groups")
cat("Kurskal-Wallis non-parametric test for 50 iterations of amount\n")
oprint(aggregate(out[[1]]["p.value"], out[[1]][c("Country","Fish")], mean))
cat("Mann-Whitney U non-parametric test for 50 iterations of amount\n")
oprint(aggregate(out[[2]]["p.value"], out[[2]][c("Pair1","Pair2","Country","Fish")], mean))
###################### Logistic regression
dat2 <- dat[dat$Iter==1 & dat$Fish=="Herring",]
tmp <- printregr(dat2, "Eat.fish", "What explains whether people eat fish at all?")
tmp <- printregr(dat2, "Eat.herring", "What explains whether people eat herring compared with other fish?")
dat2$Eat.BH2 <- ifelse(dat2$Eat.BH!="Yes" | is.na(dat2$Eat.BH),0,1)
tmp <- printregr(dat2, "Eat.BH2", "What explains whether people eat Baltic herring compared with everyone else?")
dat2$What2 <- ifelse(dat2$What=="Baltic herring",1,ifelse(dat2$What=="Other herring",0,NA))
tmp <- printregr(dat2, "What2", "What explains whether people eat Baltic herring compared with other herring?")
cat("###################### How fish-specific causes are explained by
    population determinants?\n")
for(j in c(
  "Tastes.good.BH",
  "Self.caught.BH",
  "Easy.to.cook.BH",
  "Quick.to.cook.BH",
  "Readily.available.BH",
  "Healthy.BH",
  "Inexpensive.BH",
  "Family.likes.it.BH",
  "Environmental.BH",
  "Traditional.BH"
)) {
  tmp <- printregr(dat2, j, "What explains the reason to eat Baltic herring?")
}
for(j in c(
  "Bad.taste.BH",
  "Not.used.to.BH",
  "Cannot.cook.BH",
  "Difficult.to.cook.BH",
  "Not.available.BH",
  "Health.risks.BH",
  "Better.for.feed.BH",
  "Quality.issues.BH",
  "Sustainability.BH",
  "Not.traditional.BH"
)) {
  tmp <- printregr(dat2, j, "What explains the reason to not eat Baltic herring?")
}
# tmp <- printregr(dat2, "Eat.salmon", "What explains whether people eat any salmon compared with other fish?")
# tmp <- printregr(dat2, "Baltic.salmon", "What explains whether people eat Baltic salmon compared with other salmon?")
</rcode>
</rcode>
Luke data about fish consumption in Finland [https://stat.luke.fi/en/fish-consumption-2017_en][http://statdb.luke.fi/PXWeb/pxweb/en/LUKE/LUKE__06%20Kala%20ja%20riista__06%20Muut__02%20Kalan%20kulutus/2_Kalankulutus.px/table/tableViewLayout1/?rxid=dc711a9e-de6d-454b-82c2-74ff79a3a5e0]
<t2b name="Fish consumption as food in Finland" index="Origin,Species,Year" locations="1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017" unit="kg/a per person">
domestic fish|Total|6.1|6.1|5.9|6.2|5.8|5.3|5.2|5.0|5.0|4.4|4.5|4.3|3.8|3.8|3.8|4.0|4.1|4.1|4.1
domestic fish|Farmed rainbow trout|1.6|1.6|1.6|1.6|1.3|1.3|1.4|1.1|1.1|1.2|1.2|1.1|1.0|1.0|1.1|1.1|1.3|1.2|1.2
domestic fish|Baltic herring|0.8|1.2|1.1|1.1|0.9|0.8|0.7|0.5|0.4|0.4|0.4|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3
domestic fish|Pike|0.8|0.7|0.7|0.7|0.6|0.7|0.7|0.8|0.7|0.6|0.6|0.6|0.5|0.4|0.4|0.5|0.5|0.4|0.4
domestic fish|Perch|0.7|0.7|0.7|0.7|0.6|0.6|0.6|0.7|0.7|0.5|0.5|0.5|0.4|0.4|0.4|0.5|0.5|0.4|0.4
domestic fish|Vendace|0.7|0.7|0.7|0.7|0.8|0.8|0.7|0.6|0.6|0.6|0.6|0.7|0.6|0.6|0.6|0.6|0.6|0.5|0.6
domestic fish|European whitefish|0.4|0.4|0.4|0.4|0.3|0.3|0.3|0.3|0.5|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.2|0.3|0.3
domestic fish|Pike perch|0.3|0.2|0.2|0.2|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.3|0.4|0.4
domestic fish|Other domestic fish|0.8|0.6|0.5|0.8|1.0|0.5|0.5|0.7|0.7|0.5|0.6|0.5|0.4|0.5|0.4|0.5|0.4|0.5|0.5
imported fish|Total|6.0|6.2|7.0|7.1|8.0|8.6|7.9|8.6|9.7|9.7|9.3|10.2|11.1|10.9|10.8|10.9|10.2|9.1|9.8
imported fish|Farmed rainbow trout|0.2|0.3|0.4|0.6|0.9|0.6|0.6|0.7|1.0|1.0|0.8|0.8|0.9|1.0|0.9|0.9|0.8|0.9|0.8
imported fish|Farmed salmon|1.0|0.9|1.2|1.3|1.6|2.2|1.9|2.0|2.7|2.6|2.9|3.1|3.9|4.2|4.0|4.4|4.1|3.5|4.0
imported fish|Tuna (prepared and preserved)|1.0|1.2|1.4|1.4|1.5|1.6|1.6|1.5|1.7|1.7|1.6|1.7|1.7|1.6|1.9|1.7|1.6|1.4|1.5
imported fish|Saithe (frozen fillet)|0.7|0.6|0.7|0.7|0.7|0.6|0.4|0.5|0.5|0.5|0.5|0.5|0.6|0.5|0.5|0.5|0.5|0.4|0.4
imported fish|Shrimps|0.5|0.4|0.5|0.5|0.5|0.5|0.5|0.5|0.6|0.6|0.6|0.7|0.7|0.6|0.6|0.5|0.5|0.4|0.4
imported fish|Herring and Baltic herring (preserved)|0.6|0.5|0.6|0.6|0.5|0.4|0.5|0.6|0.5|0.6|0.3|0.5|0.4|0.3|0.4|0.5|0.5|0.5|0.5
imported fish|Other imported fish|2.0|2.3|2.2|2.0|2.3|2.7|2.4|2.8|2.7|2.7|2.6|2.9|2.9|2.7|2.5|2.3|2.3|1.9|2.2
</t2b>


==== Descriptive statistics ====
==== Descriptive statistics ====
Line 1,022: Line 1,684:
==== Initiate ovariables ====
==== Initiate ovariables ====


'''Amount estimated from a bayesian model.
=====jsp taken directly from data WITHOUT salmpling=====
 
<rcode name="nonsamplejsp" label="Initiate jsp from data without sampling (for developers only)" embed=1>
# This is code Op_en7749/nonsamplejsp on page [[Goherr: Fish consumption study]]
# The code produces amount esimates (jsp ovariable) directly from data rather than bayesian model or sampling.
 
library(OpasnetUtils)
 
jsp <- Ovariable(
  "jsp",
  dependencies = data.frame(Name="survey1", Ident="Op_en7749/preprocess2"), # [[Goherr: Fish consumption study]]
  formula = function(...) {
    require(reshape2)
   
    sur <- survey1[c(157,1,3,158,162:178,81,82,131,132,154)] # Removed, not needed:16,29,30?
   
    #colnames(sur)
    #[1] "Row"                      "Country"                  "Gender"                 
    #[4] "Ages"                    "Baltic.salmon.n"          "How.often.BS.n"         
    #[7] "How.much.BS.n"            "How.often.side.BS.n"      "How.much.side.BS.n"     
    #[10] "Better.availability.BS.n" "Less.chemicals.BS.n"      "Eat.BH.n"               
    #[13] "How.often.BH.n"          "How.much.BH.n"            "How.often.side.BH.n"   
    #[16] "How.much.side.BH.n"      "Better.availability.BH.n" "Less.chemicals.BH.n"   
    #[19] "Eatfish"                  "Eatsalm"                  "Eatherr"               
    #[22] "Recommended.BS"          "Not.recommended.BS"      "Recommended.BH"         
    #[25] "Not.recommended.BH"      "Weighting"             
   
    # Make sure that Row is kept separate from Iter because in the sampling version they are different.
    # sur contained columns Eat.fish, How.often.fish, Eat.salmon. Are these needed, as all other questions are melted? No
    # Columns 5-7 removed, so colnames list above does not match.
   
    colnames(sur)[5] <- "Eat.BS"
    colnames(sur) <- gsub("\\.n","",colnames(sur))
    sur[22:25] <- sapply(sur[22:25], as.numeric)
    sur <- melt(
      sur,
      id.vars=c(1:4,26),
      variable.name="Question",
      value.name="Result"
    )
    sur$Fish <- ifelse(grepl("BH",sur$Question),"Herring","Salmon")
    sur$Question <- gsub("\\.BS","",sur$Question)
    sur$Question <- gsub("\\.BH","",sur$Question)
   
    ########## If the missing values are not adjusted, they drop out in the next stage.   
    if(TRUE) {
      # The adjustments below probably should go to the preprocess2 code.
      sur$Result[sur$Question %in% c(
        "Better.availability",
        "Less.chemicals",
        "Recommended",
        "Not.recommended"
      ) & is.na(sur$Result)] <- 4 # If missing --> no change 
      sur$Result[is.na(sur$Result)] <- 1 # Replace missing values with 1. That will produce 0 g/d.
    }
   
    return(Ovariable(
      output = sur,
      marginal = colnames(sur) %in% c("Fish", "Iter", "Question")
    ))
  }
)
 
objects.store(jsp)
cat("Ovariable jsp with actual survey data: every respondent is kept in data without sampling.\n")
</rcode>
 
===== Amount estimated from a bayesian model =====


* Model run 24.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=3UORzPwospQxp82h]
* Model run 24.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=3UORzPwospQxp82h]
Line 1,066: Line 1,795:
</rcode>
</rcode>


'''Amount estimates directly from data rather than from a bayesian model.
===== Amount estimates directly from data rather than from a bayesian model =====


* Initiation run 18.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=crW1kboP72BN1JbK]
* Initiation run 18.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=crW1kboP72BN1JbK]
Line 1,137: Line 1,866:
</rcode>
</rcode>


'''Initiate other ovariables
===== Initiate other ovariables =====


* Code stores ovariables assump, often, much, oftenside, muchside, amount.
* Code stores ovariables assump, often, much, oftenside, muchside, amount.
Line 1,266: Line 1,995:
)
)


effrecomm <- Ovariable( # Effect of recommendations
effrecommRaw <- Ovariable( # Effect of recommendations
   "effrecomm",
   "effrecommRaw",
   dependencies = data.frame(
   dependencies = data.frame(
     Name=c("jsp","assump"),
     Name=c("jsp","assump"),
Line 1,298: Line 2,027:


     return(Ovariable(output=out, marginal=colnames(out) %in% c("Iter","Fish","Recommendation")))
     return(Ovariable(output=out, marginal=colnames(out) %in% c("Iter","Fish","Recommendation")))
  },
  unit = "change of amount in fractions"
)
# effrecommRaw and effrecomm are needed to enable decisions to affect the effect of recommendation before
# the ovariable is collapsed. So, target the decisions to effrecommRaw and collapses to effrecomm.
effrecomm <- Ovariable(
  "effrecomm",
  dependencies = data.frame(Name="effrecommRaw"),
  formula = function(...) {
    return(effrecommRaw)
   }
   }
)
)
Line 1,349: Line 2,090:
objects.store(list=ls())
objects.store(list=ls())
cat("Ovariables", ls(), "stored.\n")
cat("Ovariables", ls(), "stored.\n")
</rcode>
==== Other code ====
This is code for analysing EFSA food intake data about fish for BONUS GOHERR manuscript Ronkainen L, Lehikoinen A, Haapasaari P, Tuomisto JT. 2019.
<rcode>
################################
# This is code for analysing EFSA food intake data about fish for manuscript
# Ronkainen L, Lehikoinen A, Haapasaari P, Tuomisto JT. 2019
library(ggplot2)
# Change the encoding from UCS-2 to UTF-8 first.
dat <- read.csv(
#  "C:/_Eivarmisteta/EFSA fish chronic intake.csv",
  "C:/_Eivarmisteta/EFSA herring salmon chronic intake.csv",
  header=TRUE, sep=",",dec=".",quote='\"'
)
#> colnames(dat)
#[1] "ï..Survey.s.country"    "Survey.start.year"      "Survey"               
#[4] "Population.Group..L2."  "pop.order"              "Exposure.hierarchy..L1."
#[7] "Number.of.subjects"      "Number.of.consumers"    "Mean"                 
#[10] "Standard.Deviation"      "X5th.percentile"        "X10th.percentile"     
#[13] "Median"                  "X95th.percentile"        "X97.5th.percentile"   
#[16] "X99th.percentile"        "Comment"               
colnames(dat)[c(1,2,4,12)] <- c("Country","Year","Group","Fish")
#> levels(dat$Group)
#[1] "Adolescents"    "Adults"          "Elderly"        "Infants"        "Lactating women"
#[6] "Other children"  "Pregnant women"  "Toddlers"        "Very elderly" 
levels(dat$Group) <- c("Children","Adults","Elderly","Children","Pregnant women",
                      "Children","Pregnant women","Children","Elderly")
dat <- dat[dat$Fish !="Diadromous fish",] # Foodex2 level 7 entries
ggplot(dat,aes(x=Country, y=Mean, colour=Group))+
  geom_point()+coord_flip()+facet_wrap(~Fish, scales="free")+
  labs(title="Fish consumption in European countries",
      y="Average consumption (g/d)")
</rcode>
</rcode>



Latest revision as of 15:39, 12 August 2020

Progression class
In Opasnet many pages being worked on and are in different classes of progression. Thus the information on those pages should be regarded with consideration. The progression class of this page has been assessed:
This page is reviewed
It has been read with a critical eye and commented on by an outside source, and given impairment suggestions have been included in the page. An equivalent to a peer-reviewed article.
The content and quality of this page is/was being curated by the project that produced the page.

The quality was last checked: 2019-08-26.



This page contains a detailed description about data management and analysis of an international survey related to scientific article Forage Fish as Food: Consumer Perceptions on Baltic Herring by Mia Pihlajamäki, Arja Asikainen, Suvi Ignatius, Päivi Haapasaari, and Jouni T. Tuomisto.[1] The results of this survey where also used in another article Health effects of nutrients and environmental pollutants in Baltic herring and salmon: a quantitative benefit-risk assessment by the same group.[2]

Question

How Baltic herring and salmon are used as human food in Baltic sea countries? Which determinants affect on people’s eating habits of these fish species?

Answer

Rationale

Survey of eating habits of Baltic herring and salmon in Denmark, Estonia, Finland and Sweden has been done in September 2016 by Taloustutkimus oy. Content of the questionnaire can be accessed in Google drive. The actual data can be found from the link below (see Data).

Data

Questionnaire

Original datafile File:Goherr fish consumption.csv.



Assumptions

The following assumptions are used:

Assumptions for calculations(-)
ObsVariableValueUnitResultDescription
1freq1times /a0Never
2freq2times /a0.5 - 0.9less than once a year
3freq3times /a2 - 5A few times a year
4freq4times /a12 - 361 - 3 times per month
5freq5times /a52once a week
6freq6times /a104 - 2082 - 4 times per week
7freq7times /a260 - 3645 or more times per week
8amdish1g /serving20 - 701/6 plate or below (50 grams)
9amdish2g /serving70 - 1301/3 plate (100 grams)
10amdish3g /serving120 - 1801/2 plate (150 grams)
11amdish4g /serving170 - 2302/3 plate (200 grams)
12amdish5g /serving220 - 2805/6 plate (250 grams)
13amdish6g /serving270 - 400full plate (300 grams)
14amdish7g /serving400 - 550overly full plate (500 grams)
15ingredientfraction0.1 - 0.3Fraction of fish in the dish
16amside1g /serving20 - 701/6 plate or below (50 grams)
17amside2g /serving70 - 1301/4 plate (100 grams)
18amside3g /serving120 - 1801/2 plate (150 grams)
19amside4g /serving170 - 2302/3 plate (200 grams)
20amside5g /serving220 - 2805/6 plate (250 grams)
21change1fraction-1 - -0.8Decrease it to zero
22change2fraction-0.9 - -0.5Decrease it to less than half
23change3fraction-0.6 - -0.1Decrease it a bit
24change4fraction0No effect
25change5fraction0.1 - 0.6Increase it a bit
26change6fraction0.5 - 0.9Increase it over by half
27change7fraction0.8 - 1.3Increase it over to double
28change8fraction-0.3 - 0.3Don't know

Preprocessing

This code is used to preprocess the original questionnaire data from the above .csv file and to store the data as a usable variable to Opasnet base. The code stores a data.frame named survey.

  • Model run 11.7.2018 [2]
  • Model run 27.3.2019 with country codes DK EE FI SE [3]

+ Show code

Analyses

  • Sketches about modelling determinants of eating (spring 2018) [4]

Figures, tables and stat analyses for the first manuscript

  • Model run 8.5.2019 with thlVerse code (not run on Opasnet) [5] [6]
  • Model run 26.8.2018 with fig 6 as table [7]
  • Previous model runs are archived.

Methodological concerns:

+ Show code

Luke data about fish consumption in Finland [8][9]

Fish consumption as food in Finland(kg/a per person)
ObsOriginSpecies1999200020012002200320042005200620072008200920102011201220132014201520162017
1domestic fishTotal6.16.15.96.25.85.35.25.05.04.44.54.33.83.83.84.04.14.14.1
2domestic fishFarmed rainbow trout1.61.61.61.61.31.31.41.11.11.21.21.11.01.01.11.11.31.21.2
3domestic fishBaltic herring0.81.21.11.10.90.80.70.50.40.40.40.30.30.30.30.30.30.30.3
4domestic fishPike0.80.70.70.70.60.70.70.80.70.60.60.60.50.40.40.50.50.40.4
5domestic fishPerch0.70.70.70.70.60.60.60.70.70.50.50.50.40.40.40.50.50.40.4
6domestic fishVendace0.70.70.70.70.80.80.70.60.60.60.60.70.60.60.60.60.60.50.6
7domestic fishEuropean whitefish0.40.40.40.40.30.30.30.30.50.30.30.30.30.30.30.30.20.30.3
8domestic fishPike perch0.30.20.20.20.30.30.30.30.30.30.30.30.30.30.30.30.30.40.4
9domestic fishOther domestic fish0.80.60.50.81.00.50.50.70.70.50.60.50.40.50.40.50.40.50.5
10imported fishTotal6.06.27.07.18.08.67.98.69.79.79.310.211.110.910.810.910.29.19.8
11imported fishFarmed rainbow trout0.20.30.40.60.90.60.60.71.01.00.80.80.91.00.90.90.80.90.8
12imported fishFarmed salmon1.00.91.21.31.62.21.92.02.72.62.93.13.94.24.04.44.13.54.0
13imported fishTuna (prepared and preserved)1.01.21.41.41.51.61.61.51.71.71.61.71.71.61.91.71.61.41.5
14imported fishSaithe (frozen fillet)0.70.60.70.70.70.60.40.50.50.50.50.50.60.50.50.50.50.40.4
15imported fishShrimps0.50.40.50.50.50.50.50.50.60.60.60.70.70.60.60.50.50.40.4
16imported fishHerring and Baltic herring (preserved)0.60.50.60.60.50.40.50.60.50.60.30.50.40.30.40.50.50.50.5
17imported fishOther imported fish2.02.32.22.02.32.72.42.82.72.72.62.92.92.72.52.32.31.92.2

Descriptive statistics

Correlation matrix of all questions in the survey (answers converted to numbers).

Model must contain predictors such as country, gender, age etc. Maybe we should first study what determinants are important? Model must also contain determinants that would increase or decrease fish consumption. This should be conditional on the current consumption. How? Maybe we should look at principal coordinates analysis with all questions to see how they behave.

Also look at correlation table to see clusters.

Some obvious results:

  • If reports no fish eating, many subsequent answers are NA.
  • No vitamins correlates negatively with vitamin intake.
  • Unknown salmon correlates negatively with the types of salmon eaten.
  • Different age categories correlate with each other.

However, there are also meaningful negative correlations:

  • Country vs allergy
  • Country vs Norwegian salmon and Rainbow trout
  • Country vs not traditional.
  • Country vs recommendation awareness
  • Allergy vs economic wellbeing
  • Baltic salmon use (4 questions) vs Don't like taste and Not used to
  • All questions between Easy to cook ... Traditional dish

Meaningful positive correlations:

  • All questions between Baltic salmon ... Rainbow trout
  • How often Baltic salmon/herring/side salmon/side herring
  • How much Baltic salmon/herring/side salmon/side herring
  • Better availability ... Recommendation
  • All questions between Economic wellbeing...Personal aims
  • Omega3, Vitamin D, and Other vitamins

Model runs

+ Show code

Bayes model

  • Model run 3.3.2017. All variables assumed independent. [11]
  • Model run 3.3.2017. p has more dimensions. [12]
  • Model run 25.3.2017. Several model versions: strange binomial+multivarnormal, binomial, fractalised multivarnormal [13]
  • Model run 27.3.2017 [14]
  • Other models except multivariate normal were archived and removed from active code 29.3.2017.
  • Model run 29.3.2017 with raw data graphs [15]
  • Model run 29.3.2017 with salmon and herring ovariables stored [16]
  • Model run 13.4.2017 with first version of coordinate matrix and principal coordinate analysis [17]
  • Model run 20.4.2017 [18] code works but needs a safety check against outliers
  • Model run 21.4.2017 [19] some model results plotted
  • Model run 21.4.2017 [20] ovariables produced by the model stored.
  • Model run 18.5.2017 [21] small updates
  • 13.2.2018 old model run but with new Opasnet [22]

+ Show code

Initiate ovariables

jsp taken directly from data WITHOUT salmpling

+ Show code

Amount estimated from a bayesian model
  • Model run 24.5.2017 [23]

+ Show code

Amount estimates directly from data rather than from a bayesian model
  • Initiation run 18.5.2017 [24]
  • Initiation run 24.2.2018: sampling from survey rather than each respondent once [25]

+ Show code

Initiate other ovariables
  • Code stores ovariables assump, often, much, oftenside, muchside, amount.
  • Model run 19.5.2017 [26]
  • Initiation run 24.5.2017 without jsp [27]
  • Model run 8.6.2017 [28]

+ Show code

Other code

This is code for analysing EFSA food intake data about fish for BONUS GOHERR manuscript Ronkainen L, Lehikoinen A, Haapasaari P, Tuomisto JT. 2019.

+ Show code

Dependencies

The survey data will be used as input in the benefit-risk assessment of Baltic herring and salmon intake, which is part of the WP5 work in Goherr-project.

See also

Keywords

References

  1. Jump up to: 1.0 1.1 Mia Pihlajamäki, Arja Asikainen, Suvi Ignatius, Päivi Haapasaari and Jouni T. Tuomisto. Forage Fish as Food: Consumer Perceptions on Baltic Herring. Sustainability 2019, 11(16), 4298; https://doi.org/10.3390/su11164298
  2. Tuomisto, J.T., Asikainen, A., Meriläinen, P., Haapasaari, P. Health effects of nutrients and environmental pollutants in Baltic herring and salmon: a quantitative benefit-risk assessment. BMC Public Health 20, 64 (2020). https://doi.org/10.1186/s12889-019-8094-1

Related files

Goherr: Fish consumption study. Opasnet . [29]. Accessed 15 May 2025.