Goherr: Fish consumption study: Difference between revisions

From Opasnet
Jump to navigation Jump to search
(added references)
 
(21 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 288: Line 247:
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 11.7.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=Fa8TMg6h8DbEmhjp]
* 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 452: Line 417:
==== Figures, tables and stat analyses 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]
* Model run 20.6.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=bWeOP4teZT5c4T9Z]
* 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].
* Model run 11.7.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=yakkgFPuWKkXu3qF]
 
* Model run 9.8.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=r3UiLe8x2sRz8DpR]
Methodological concerns:
* Model run 12.8.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=zZWea5znWTliAEHo]
* [https://stackoverflow.com/questions/12953045/warning-non-integer-successes-in-a-binomial-glm-survey-packages Is the warning in logistic regression important?]
* Model run 12.9.2018 with country-specific regression analyses, open-ended answers and fig about eating any herring [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=Z4hgFmuaqm7mNulZ]
* [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 467: Line 435:
library(reshape2)
library(reshape2)
library(MASS)
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)


groups <- function(o) {
colors <- c("#74AF59FF","#2F62ADFF","#D692BDFF","#29A0C1FF","#BE3F72FF","#FBB848FF")
   o$Group <- paste(o$Gender, o$Ages)
 
########################## 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) {
   o$Group <- paste(o$Gender, o$Ages)
   o$Group <- factor(o$Group, levels = c(
   o$Group <- factor(o$Group, levels = c(
     "Female 18-45",
     "Female 18-45",
Line 486: Line 574:
#' @param dat a data.frame containing survey answers for the target subgroup.
#' @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
#' @param title an atomic string with the name to be shown as the title of graph
#' @return ggplot object


reasons <- function(
reasons <- function(
Line 491: Line 580:
   title
   title
) {
) {
  require(reshape2)
   weight <- Ovariable("weight",data=data.frame(
   weight <- Ovariable("weight",data=data.frame(
     Row=dat$Row,
     Row=dat$Row,
Line 518: Line 608:
   oprint(tmp@output)
   oprint(tmp@output)
    
    
   tmp <- ggplot(tmp@output, aes(x=Reason, y=Result,colour=Country, group=Country))+
   if(thl) {
    geom_point(shape=21, size=6, fill="Grey", stroke=2)+
    tmp <- thlLinePlot(tmp@output, xvar=Reason, yvar=Result,groupvar=Country,
    geom_line()+
                      colors= c("#519B2FFF", "#2F62ADFF", "#BE3F72FF","#88D0E6FF"), # #29A0C1FF"),
    coord_flip()+
                      # THL colors but fourth is brigter
    theme_gray(base_size=BS)+
                      legend.position = c(0.85,0.2), base.size = BS, title=title,
    scale_y_continuous(labels=scales::percent_format())+
                      subtitle="Fraction of population")+
    labs(
      coord_flip()+
      title=title,
      scale_y_continuous(labels=scales::percent_format(accuracy=1))
      x="Answer",
  } else {
      y="Fraction of population")
    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)
   return(tmp)
}
}


impacts <- function(dat,title) {
#' @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(
   tmp <- melt(
     dat,
     dat,
     id.vars = "Country",
     id.vars = c("Country","Weighting"),
     variable.name ="Reason",
     variable.name ="Reason",
     value.name="Response"
     value.name="Response"
Line 563: Line 670:
   )] <- "Increase"
   )] <- "Increase"
    
    
   popu <- aggregate(as.numeric(tmp$Response), tmp["Reason"],sum)
  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)]
   popu <- popu$Reason[order(popu$x)]
   tmp$Reason <- factor(tmp$Reason, levels=popu)
   tmp$Reason <- factor(tmp$Reason, levels=popu)
   levels(tmp$Reason) <- gsub("( BH| BS)", "", gsub("\\.", " ", levels(tmp$Reason)))
   levels(tmp$Reason) <- gsub("( BH| BS)", "", gsub("\\.", " ", levels(tmp$Reason)))
    
    
  tmp$Result <- 1
   tmp <- (oapply(tmp, c("Reason","Response"), sum)/
  tmp <- EvalOutput(Ovariable("tmp",data=tmp))
             oapply(tmp, c("Reason"),sum))
   tmp <- (oapply(tmp, c("Country","Reason","Response"), sum)/
             oapply(tmp, c("Country","Reason"),sum))
    
    
  oprint(tmp@output)
    
    
   tmp <- ggplot(tmp@output,
   print(ggplot(tmp@output,
                aes(x=Reason, weight=Result, fill=Response))+geom_bar()+
              aes(x=Reason, weight=Result, fill=Response))+geom_bar()+
    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")+
    scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
          scale_fill_manual(values=c("Grey",colors[c(1,6,5)]))+
    scale_y_continuous(breaks=c(0.5,1),labels=scales::percent_format())+
          scale_y_continuous(breaks=c(0.5,1),labels=scales::percent_format())+
    labs(
          labs(
      title=title,
            title=title,
      x="Cause",
            x="Cause",
      y="Fraction of population"
            y="Fraction of population"
    )
          ))
    
    
   return(tmp)
   return(tmp@output)
}
}


################## Get data
#' @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


objects.latest("Op_en7749", "preprocess2") # [[Goherr: Fish consumption study]]: survey1, surv
or95 <- function(fit) {
objects.latest("Op_en7749", "nonsamplejsp") # [[Goherr: Fish consumption study]]: jsp every respondent from data exactly once.
  cat("Odds ratios with 95 % confidence intervals\n")
objects.latest("Op_en7749", "initiate") # [[Goherr: Fish consumption study]]: amount etc.
  tmp <- exp(cbind(coef(fit), confint(fit)))
openv.setN(50)
  tmp <- data.frame(
survey1 <- groups(survey1)
    rownames(tmp),
effinfo <- 0 # No policies implemented.
    paste0(sprintf("%.2f",tmp[,1]), " (",sprintf("%.2f", tmp[,2]),"-",sprintf("%.2f",tmp[,3]),")")
effrecomm <- 0
  )
amount <- EvalOutput(amount)
  colnames(tmp) <- c("Parameter","OR (95% CI)")
#result(amount)[result(amount)==0] <- 0.1 Not needed if not log-transformed.
  return(tmp)
}


#### General fish eating
#' @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


su1 <- EvalOutput(Ovariable(
printregr <- function(dat, y, title, countries = TRUE) {
   "su1",
    
   data = data.frame(
   dat$y <- dat[[y]]
    survey1[c("Country","Gender","Education","Purchasing.power")],
  cat("\n############### Logistic regression analysis:",title, "(", y, ")\n")
     Result=1
  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)


#table 2: statistics of the survey population in each country (n, female %, education, purchasing power)
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
 
su1 <- EvalOutput(Ovariable(
  "su1",
  data = data.frame(
    survey1[c("Country","Gender","Education","Purchasing.power")],
    Result=1
  )
))
 
# Figure 1. Origin of consumed fish in Finland between 1999 and 2016.
# Data not on this page, drawn separately.
 
# Table 1. Dimensions of embeddedness, modified from Hass (2007, p. 16)
# Written directly on the manuscript.
 
# Table 2: statistics of the survey population in each country (n, female %, education, purchasing power)


tmp <- round(t(data.frame(
tmp <- round(t(data.frame(
Line 654: Line 856:
   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 population")
   labs(title="Fraction of fish eaters within population")


############## Answers to open questions
############## Answers to open questions
cat("Number of open-ended answers per question and country\n")
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 <- 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,
  tmp <- reshape(
  varying=list(
    tmp,
    Condition = colnames(tmp)[grepl("other",tolower(colnames(tmp)))],
    varying=list(
    Text = colnames(tmp)[grepl("open",colnames(tmp))]
      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"
    times = colnames(tmp)[grepl("other",tolower(colnames(tmp)))],
)
    direction="long"
tmp <- tmp[tmp$Why.not.fish.open!="" , ]
  )
oprint(table(tmp$time,tmp$Country))
  tmp <- tmp[tmp$Why.not.fish.open!="" , ]
oprint(tmp)
  oprint(table(tmp$time,tmp$Country))
  oprint(tmp)
} # End if(FALSE)


#### Baltic salmon
#### Baltic salmon
Line 680: Line 885:
   "Reasons not to eat among fish non-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,154,157)], # 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 used among salmon consumers"
   "Species used among salmon consumers"
)
)
reasons(
reasons(
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,50:59,154,157)], # 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 to eat 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:70,154,157)], # 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 not to eat among Baltic salmon non-consumers"
   "Reasons not to eat, Baltic salmon non-consumers"
)
)


impacts(
if(localcomp) ggsave("Figure 6.pdf", width=10,height=5)
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,73:82)],
if(localcomp) ggsave("Figure 6.png", width=10,height=5)
   "Effect on Baltic salmon consumption"
 
impacts.sal <- impacts(
   survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,73:82,154)],
   "Effect on Baltic salmon consumption",
  population
)
)
####### Any herring
####### Any herring


Line 706: Line 925:
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
   scale_fill_manual(values=colors)+
   labs(
   labs(
     title="Any herring consumption",
     title="Any herring consumption",
Line 714: Line 933:
####### Baltic herring
####### Baltic herring


ggplot(survey1[survey1$Eat.fish=="Yes",],
### Figure 2. Comparison of Baltic herring consumption habits.
      aes(x=Country, weight=Weighting, fill=Eat.BH))+geom_bar(position="fill")+
 
  coord_flip()+#facet_grid(.~Country)+
survey1$What <- ifelse(is.na(survey1$Eat.BH), "Other fish", as.character(survey1$Eat.BH))
  theme_gray(base_size=BS)+
survey1$What <- factor(
  theme(legend.position = "bottom")+
  survey1$What,
  scale_y_continuous(labels=scales::percent_format())+
  levels=(c("Yes","I don't know","No","Other fish")),
  #  scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
  labels=(c("Baltic herring","Some herring","Other herring","Other fish"))
  labs(
)
    title="Baltic herring consumption",
 
    y="Fraction of fish consumers"
tmp <- survey1[survey1$Eat.fish=="Yes",]
   )
colnames(tmp)[colnames(tmp)=="Weighting"] <- "Result"
tmp <- EvalOutput(Ovariable("tmp", data=tmp))
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(survey1[survey1$Eat.herring=="Yes",],
ggplot(survey1[survey1$Eat.herring=="Yes",],
Line 732: Line 989:
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
   scale_fill_manual(values=c(colors))+
   labs(
   labs(
     title="Baltic herring consumption",
     title="Baltic herring consumption",
Line 746: Line 1,003:
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
   scale_fill_manual(values=colors)+
   labs(
   labs(
     title="Any salmon consumption",
     title="Any salmon consumption",
Line 760: Line 1,017:
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
   scale_fill_manual(values=colors)+
   labs(
   labs(
     title="Baltic salmon consumption",
     title="Baltic salmon consumption",
Line 772: Line 1,029:
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
   scale_fill_manual(values=colors)+
   labs(
   labs(
     title="Baltic salmon consumption",
     title="Baltic salmon consumption",
Line 787: Line 1,044:
   theme(legend.position = "bottom")+
   theme(legend.position = "bottom")+
   scale_y_continuous(labels=scales::percent_format())+
   scale_y_continuous(labels=scales::percent_format())+
   scale_fill_manual(values=c("Grey","Green","Yellow","Red"))+
   scale_fill_manual(values=colors)+
   labs(title="Awareness of food recommendations about Baltic fish")
   labs(title="Awareness of food recommendations about Baltic fish")
# Figure 3. Percentages of reasons to eat Baltic herring


reasons(
reasons(
Line 794: Line 1,053:
   "Reasons to eat 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:120,154,157)], # 1 Country, 28 Why not Baltic herring
   survey1[
   "Reasons not to eat 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
)
)


impacts(
tmp <- rbind(
   survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,123:132)], # 1 Country, 29 Influence of Baltic herring policies
   data.frame(
  title="Effect on Baltic herring consumption"
    Fish="Herring",
    impacts.herr
  ),
  data.frame(
    Fish="Salmon",
    impacts.sal
  )
)
)
oprint(aggregate(tmp$Result,tmp[c("Reason","Response","Fish")],sum))
if(thl) {
  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
######### Amounts estimated for each respondent
Line 815: Line 1,133:
#  ]
#  ]
tmp <- cbind(
tmp <- cbind(
   oapply(amount, c("Country","Group","Fish"), FUN=mean)@output,
   oapply(amount, c("Country","Group","Gender","Ages","Fish"), FUN=mean)@output,
   SD = oapply(amount, c("Country","Gender","Ages","Fish"), FUN=sd)$amountResult
   SD=oapply(amount, c("Country","Group","Gender","Ages","Fish"), FUN=sd)$amountResult
)
)  
ggplot(tmp, aes(x=Group, weight=amountResult, fill=Fish))+
 
   geom_bar()+facet_grid(Country~Fish)+
# 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)+
  scale_fill_manual(values=colors)+
  guides(fill=FALSE)+
   theme(axis.text.x = element_text(angle = -90))+
   theme(axis.text.x = element_text(angle = -90))+
   labs(
   labs(
     title="Fish consumption in subgroups",
     title="Baltic herring consumption in subgroups",
     y="Average consumption (g/day)")
     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$Result <- paste0(sprintf("%.1f",tmp$amountResult), " (",sprintf("%.1f", tmp$SD,1),")")
Line 832: Line 1,159:
oprint(tmp)
oprint(tmp)


result(amount)[result(amount)==0] <- 0.1
weight <- EvalOutput(Ovariable("weight",data=data.frame(
ggplot(amount@output, aes(x=amountResult, colour=Group))+stat_ecdf()+scale_x_log10()
  amount@output[c("Row","Iter","Fish","Country")],
  Result=amount$Weighting
)))
 
tmp <- (oapply(amount * weight, c("Fish","Country"), sum) / oapply(weight,c("Fish","Country"), sum))@output
 
cat("Average fish consumption per country (kg/year)\n")
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)


############################
############################
Line 839: Line 1,177:
#### Statistical analyses
#### Statistical analyses


dat <- merge(survey1,amount@output)
dat <- merge(survey1,amount@output) # [amount$Fish=="Herring",]) # This results in 50 * 2 times of rows


# These should be corrected to preprocess2. NOT ordered.
# These should be corrected to preprocess2. NOT ordered.
Line 897: Line 1,235:
###################### Logistic regression
###################### Logistic regression


cat("###################### What explains whether people eat fish at all?\n")
dat2 <- dat[dat$Iter==1 & dat$Fish=="Herring",]
fit <- glm(Eat.fish ~ Ages + Gender + Country + Education + Purchasing.power, family=binomial(),data=dat[dat$Iter==1,])
 
oprint(summary(fit))
tmp <- printregr(dat2, "Eat.fish", "What explains whether people eat fish at all?")
step <- stepAIC(fit, direction = "both")
oprint(step$anova)


for(i in unique(dat$Country)) {
tmp <- printregr(dat2, "Eat.herring", "What explains whether people eat herring compared with other fish?")
  cat("\n#### Country-specific regression analysis:", i, "\n")
  fit <- glm(Eat.fish ~ Ages + Gender + Education + Purchasing.power, family=binomial(),
            data=dat[dat$Iter==1 & dat$Country==i,])
  oprint(summary(fit))
  step <- stepAIC(fit, direction = "both")
  oprint(step$anova)
}


cat("###################### What explains whether people eat any herring?\n")
dat2$Eat.BH2 <- ifelse(dat2$Eat.BH!="Yes" | is.na(dat2$Eat.BH),0,1)
fit <- glm(Eat.herring ~ Ages + Gender + Country + Education + Purchasing.power, family=binomial(),data=dat[dat$Iter==1,])
tmp <- printregr(dat2, "Eat.BH2", "What explains whether people eat Baltic herring compared with everyone else?")
oprint(summary(fit))
step <- stepAIC(fit, direction = "both")
oprint(step$anova)


for(i in unique(dat$Country)) {
dat2$What2 <- ifelse(dat2$What=="Baltic herring",1,ifelse(dat2$What=="Other herring",0,NA))
  cat("\n#### Country-specific regression analysis:", i, "\n")
tmp <- printregr(dat2, "What2", "What explains whether people eat Baltic herring compared with other herring?")
  fit <- glm(Eat.herring ~ Ages + Gender + Education + Purchasing.power, family=binomial(),
            data=dat[dat$Iter==1 & dat$Country==i,])
  oprint(summary(fit))
  step <- stepAIC(fit, direction = "both")
  oprint(step$anova)
}


cat("###################### What explains whether people eat Baltic herring?\n")
fit <- glm(Eat.BH ~ Ages + Gender + Country + Education + Purchasing.power, family=binomial(),data=dat[dat$Iter==1,])
oprint(summary(fit))
step <- stepAIC(fit, direction = "both")
oprint(step$anova)


for(i in unique(dat$Country)) {
cat("###################### How fish-specific causes are explained by
  cat("\n#### Country-specific regression analysis:", i, "\n")
    population determinants?\n")
  fit <- glm(Eat.BH ~ Ages + Gender + Education + Purchasing.power, family=binomial(),
            data=dat[dat$Iter==1 & dat$Country==i,])
  oprint(summary(fit))
  step <- stepAIC(fit, direction = "both")
  oprint(step$anova)
}


cat("###################### What explains whether people eat any salmon?\n")
fit <- glm(Eat.salmon ~ Ages + Gender + Country + Education + Purchasing.power, family=binomial(),data=dat[dat$Iter==1,])
oprint(summary(fit))
step <- stepAIC(fit, direction = "both")
oprint(step$anova)


for(i in unique(dat$Country)) {
for(j in c(
   cat("\n#### Country-specific regression analysis:", i, "\n")
   "Tastes.good.BH",
   fit <- glm(Eat.salmon ~ Ages + Gender + Education + Purchasing.power, family=binomial(),
  "Self.caught.BH",
            data=dat[dat$Iter==1 & dat$Country==i,])
  "Easy.to.cook.BH",
   oprint(summary(fit))
  "Quick.to.cook.BH",
   step <- stepAIC(fit, direction = "both")
  "Readily.available.BH",
  oprint(step$anova)
   "Healthy.BH",
  "Inexpensive.BH",
  "Family.likes.it.BH",
  "Environmental.BH",
   "Traditional.BH"
)) {
   tmp <- printregr(dat2, j, "What explains the reason to eat Baltic herring?")
}
}


cat("###################### What explains whether people eat Baltic salmon?\n")
for(j in c(
fit <- glm(Baltic.salmon ~ Ages + Gender + Country + Education + Purchasing.power, family=binomial(),data=dat[dat$Iter==1,])
  "Bad.taste.BH",
oprint(summary(fit))
  "Not.used.to.BH",
step <- stepAIC(fit, direction = "both")
  "Cannot.cook.BH",
oprint(step$anova)
  "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?")
}


for(i in unique(dat$Country)) {
# tmp <- printregr(dat2, "Eat.salmon", "What explains whether people eat any salmon compared with other fish?")
  cat("\n#### Country-specific regression analysis:", i, "\n")
# tmp <- printregr(dat2, "Baltic.salmon", "What explains whether people eat Baltic salmon compared with other salmon?")
  fit <- glm(Baltic.salmon ~ Ages + Gender + Education + Purchasing.power, family=binomial(),
            data=dat[dat$Iter==1 & dat$Country==i,])
  oprint(summary(fit))
  step <- stepAIC(fit, direction = "both")
  oprint(step$anova)
}
</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,350: Line 1,684:
==== Initiate ovariables ====
==== Initiate ovariables ====


'''jsp taken directly from data WITHOUT salmpling.
=====jsp taken directly from data WITHOUT salmpling=====


<rcode name="nonsamplejsp" label="Initiate jsp from data without sampling (for developers only)" embed=1>
<rcode name="nonsamplejsp" label="Initiate jsp from data without sampling (for developers only)" embed=1>
Line 1,360: Line 1,694:
jsp <- Ovariable(
jsp <- Ovariable(
   "jsp",  
   "jsp",  
   dependencies = data.frame(Name="survey1", Ident="Op_en7749/preprocess2"), # [[Goherr: Food cunsumption study]]
   dependencies = data.frame(Name="survey1", Ident="Op_en7749/preprocess2"), # [[Goherr: Fish consumption study]]
   formula = function(...) {
   formula = function(...) {
     require(reshape2)
     require(reshape2)
      
      
     sur <- survey1[c(157,1,3,158,162:178,81,82,131,132)] # Removed, not needed:16,29,30?
     sur <- survey1[c(157,1,3,158,162:178,81,82,131,132,154)] # Removed, not needed:16,29,30?
      
      
     #colnames(sur)
     #colnames(sur)
Line 1,375: Line 1,709:
     #[19] "Eatfish"                  "Eatsalm"                  "Eatherr"                 
     #[19] "Eatfish"                  "Eatsalm"                  "Eatherr"                 
     #[22] "Recommended.BS"          "Not.recommended.BS"      "Recommended.BH"           
     #[22] "Recommended.BS"          "Not.recommended.BS"      "Recommended.BH"           
     #[25] "Not.recommended.BH"                        
     #[25] "Not.recommended.BH"       "Weighting"             
      
      
     # Make sure that Row is kept separate from Iter because in the sampling version they are different.
     # Make sure that Row is kept separate from Iter because in the sampling version they are different.
Line 1,386: Line 1,720:
     sur <- melt(
     sur <- melt(
       sur,
       sur,
       id.vars=c(1:4),
       id.vars=c(1:4,26),
       variable.name="Question",
       variable.name="Question",
       value.name="Result"
       value.name="Result"
Line 1,417: Line 1,751:
</rcode>
</rcode>


'''Amount estimated from a bayesian model.
===== 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,461: 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,532: 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,661: 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,693: 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,744: 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.