+ Show code- Hide code # This is code Op_en7749/ on page [[Goherr: Fish consumption study]]
library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(MASS)
#library(car)
#library(vegan)
BS <- 24
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
reasons <- function(
dat,
title
) {
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)
tmp <- ggplot(tmp@output, aes(x=Reason, y=Result,colour=Country, group=Country))+
geom_point(shape=21, size=6, fill="Grey", stroke=2)+
geom_line()+
coord_flip()+
theme_gray(base_size=BS)+
scale_y_continuous(labels=scales::percent_format())+
labs(
title=title,
x="Answer",
y="Fraction of population")
return(tmp)
}
impacts <- function(dat,title) {
tmp <- melt(
dat,
id.vars = "Country",
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"
popu <- aggregate(as.numeric(tmp$Response), tmp["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$Result <- 1
tmp <- EvalOutput(Ovariable("tmp",data=tmp))
tmp <- (oapply(tmp, c("Country","Reason","Response"), sum)/
oapply(tmp, c("Country","Reason"),sum))
tmp <- 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","Green","Yellow","Red"))+
scale_y_continuous(breaks=c(0.5,1),labels=scales::percent_format())+
labs(
title=title,
x="Cause",
y="Fraction of population"
)
return(tmp)
}
################## 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)
effinfo <- 0 # No policies implemented.
effrecomm <- 0
amount <- EvalOutput(amount)
#result(amount)[result(amount)==0] <- 0.1 Not needed if not log-transformed.
#### General fish eating
su1 <- EvalOutput(Ovariable(
"su1",
data = data.frame(
survey1[c("Country","Gender","Education","Purchasing.power")],
Result=1
)
))
#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_y_continuous(labels=scales::percent_format())+
labs(title="Fraction of fish eaters within population")
############## Answers to open questions
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)
#### 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"
)
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"
)
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 among Baltic salmon non-consumers"
)
impacts(
survey1[survey1$Baltic.salmon=="Yes"& !is.na(survey1$Baltic.salmon),c(1,73:82)],
"Effect on Baltic salmon consumption"
)
####### 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=c("Grey","Green","Yellow","Red"))+
labs(
title="Any herring consumption",
y="Fraction of fish consumers"
)
####### Baltic herring
ggplot(survey1[survey1$Eat.fish=="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("Grey","Green","Yellow","Red"))+
labs(
title="Baltic herring consumption",
y="Fraction of fish consumers"
)
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("Grey","Green","Yellow","Red"))+
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=c("Grey","Green","Yellow","Red"))+
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=c("Grey","Green","Yellow","Red"))+
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=c("Grey","Green","Yellow","Red"))+
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=c("Grey","Green","Yellow","Red"))+
labs(title="Awareness of food recommendations about Baltic fish")
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"
)
reasons(
survey1[survey1$Eat.BH=="No"& !is.na(survey1$Eat.BH),c(1,111:120,154,157)], # 1 Country, 28 Why not Baltic herring
"Reasons not to eat among Baltic herring non-consumers"
)
impacts(
survey1[survey1$Eat.BH=="Yes"& !is.na(survey1$Eat.BH),c(1,123:132)], # 1 Country, 29 Influence of Baltic herring policies
title="Effect on Baltic herring consumption"
)
######### 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","Fish"), FUN=mean)@output,
SD = oapply(amount, c("Country","Gender","Ages","Fish"), FUN=sd)$amountResult
)
ggplot(tmp, aes(x=Group, weight=amountResult, fill=Fish))+
geom_bar()+facet_grid(Country~Fish)+
theme_gray(base_size=BS)+
theme(axis.text.x = element_text(angle = -90))+
labs(
title="Fish consumption in subgroups",
y="Average consumption (g/day)")
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)
result(amount)[result(amount)==0] <- 0.1
ggplot(amount@output, aes(x=amountResult, colour=Group))+stat_ecdf()+scale_x_log10()
############################
#### Statistical analyses
dat <- merge(survey1,amount@output)
# 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
cat("###################### What explains whether people eat fish at all?\n")
fit <- glm(Eat.fish ~ 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("\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")
fit <- glm(Eat.herring ~ 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("\n#### Country-specific regression analysis:", i, "\n")
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("\n#### Country-specific regression analysis:", i, "\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)) {
cat("\n#### Country-specific regression analysis:", i, "\n")
fit <- glm(Eat.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)
}
cat("###################### What explains whether people eat Baltic salmon?\n")
fit <- glm(Baltic.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)) {
cat("\n#### Country-specific regression analysis:", i, "\n")
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)
}
| |