+ Show code- Hide code
# This is code Op_en7749/ on page [[Goherr: Fish consumption study]]
library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(MASS)
#library(extrafont) # Needed to save Arial fonts to PDF or EPS
#library(thlVerse)
#library(car)
#library(vegan)
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) {
o$Group <- paste(o$Gender, o$Ages)
o$Group <- factor(o$Group, levels = c(
"Female 18-45",
"Male 18-45",
"Female >45",
"Male >45"
))
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(
dat,
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(
dat[setdiff(colnames(dat),"Weighting")],
id.vars = c("Country","Row"),
variable.name ="Reason",
value.name="Result"
)
tmp$Result <- (tmp$Result %in% c("Yes","1"))*1
tmp <- weight * (Ovariable("tmp",data=tmp))
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 <- (oapply(tmp, c("Reason","Response"), sum)/
oapply(tmp, c("Reason"),sum))
oprint(tmp@output)
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)
}
#' @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
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(
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,
aes(x=Eat.fish, weight=Weighting, group=Country), stat="count")+
geom_bar(aes(y=..prop.., fill=Country), position="dodge")+
theme_gray(base_size=BS)+
scale_fill_manual(values=colors)+
scale_y_continuous(labels=scales::percent_format())+
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
reasons(
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
"Reasons not to eat among fish non-consumers"
)
# Figure 5
cat("Figure 5 (to be converted to text)\n")
reasons(
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"
)
reasons(
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"
)
if(localcomp) ggsave("Figure 5.pdf", width=10,height=5)
if(localcomp) ggsave("Figure 5.png", width=10,height=5)
reasons(
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, Baltic salmon non-consumers"
)
if(localcomp) ggsave("Figure 6.pdf", width=10,height=5)
if(localcomp) ggsave("Figure 6.png", width=10,height=5)
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
ggplot(survey1[survey1$Eat.fish=="Yes",],
aes(x=Country, weight=Weighting, fill=Eat.herring))+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="Any herring consumption",
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 <- 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",],
aes(x=Country, weight=Weighting, fill=Eat.BH))+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(colors))+
labs(
title="Baltic herring consumption",
y="Fraction of herring consumers"
)
####### 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(
survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,99:108,154,157)], # 1 Country, 27 why Baltic herring
"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(
survey1[
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 <- rbind(
data.frame(
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
amount <- groups(amount)
#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)+
scale_fill_manual(values=colors)+
guides(fill=FALSE)+
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)
weight <- EvalOutput(Ovariable("weight",data=data.frame(
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)
############################
#### 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?")
| |