Difference between revisions of "Goherr: Fish consumption study"

From Testiwiki
Jump to: navigation, search
(technical edits)
(Analyses: debugged. Now both codes work)
Line 351: Line 351:
  
 
=== Analyses ===
 
=== Analyses ===
 +
 +
==== Descriptive statistics ====
  
 
[[File:Goherr survey correlations.png|thumb|300px|Correlation matrix of all questions in the survey (answers converted to numbers).]]
 
[[File:Goherr survey correlations.png|thumb|300px|Correlation matrix of all questions in the survey (answers converted to numbers).]]
Line 386: Line 388:
 
Model runs
 
Model runs
 
* [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QaMJZqUX0cPaTfOF Model run 13.3.2017]
 
* [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QaMJZqUX0cPaTfOF Model run 13.3.2017]
 +
* Model run 21.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=0Lmu6OVAjomgJwAN] old code from Answer merged to this code and debugged
  
<rcode name="descrstat" label="Descriptive statistics" graphics=1>
+
<rcode label="Descriptive statistics" graphics=1>
# This is code Op_en7749/descrstat on page [[Goherr: Fish consumption study]]
+
# This is code Op_en7749/ on page [[Goherr: Fish consumption study]]
  
 
library(OpasnetUtils)
 
library(OpasnetUtils)
Line 402: Line 405:
 
for(i in c(5:6, 16, 29:30, 46:49, 85:86, 95:98, 135)) {
 
for(i in c(5:6, 16, 29:30, 46:49, 85:86, 95:98, 135)) {
 
   temp <- survey[!is.na(survey[[i]]),]
 
   temp <- survey[!is.na(survey[[i]]),]
    p <- ggplot(temp, aes(x = temp[[i]])) + geom_bar() +
+
  p <- ggplot(temp, aes(x = temp[[i]])) + geom_bar() +
 
     theme_gray(base_size = 18) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
 
     theme_gray(base_size = 18) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
 
     labs(title = colnames(temp[i])) + xlab("") + facet_wrap(~ Country)
 
     labs(title = colnames(temp[i])) + xlab("") + facet_wrap(~ Country)
Line 412: Line 415:
 
temp$Impact <- factor(temp$Impact, levels = levs, labels = levs, ordered = TRUE)
 
temp$Impact <- factor(temp$Impact, levels = levs, labels = levs, ordered = TRUE)
  
ggplot(temp, aes(x = Impact, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Changing.factor)+
+
ggplot(temp, aes(x = Impact, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Changing.factor)+
 
   theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4))
 
   theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4))
  
Line 418: Line 421:
  
 
temp <- melt(surveytemp, measure.vars = 38:43, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 38:43, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + labs(title = "Baltic salmon sources")
+
  theme_gray(base_size = 24) + labs(title = "Baltic salmon sources")
  
 
temp <- melt(surveytemp, measure.vars = 50:59, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 50:59, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + labs(title = "Reasons to eat Baltic salmon")
+
  theme_gray(base_size = 24) + labs(title = "Reasons to eat Baltic salmon")
  
 
temp <- melt(surveytemp, measure.vars = 73:82, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 73:82, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + labs(title = "Baltic salmon actions")
+
  theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + labs(title = "Baltic salmon actions")
  
 
surveytemp <- subset(survey, survey[[31]] == "No")
 
surveytemp <- subset(survey, survey[[31]] == "No")
  
 
temp <- melt(surveytemp, measure.vars = 62:70, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 62:70, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + labs(title = "Reasons not to eat Baltic salmon")
+
  theme_gray(base_size = 24) + labs(title = "Reasons not to eat Baltic salmon")
  
  
Line 439: Line 442:
  
 
temp <- melt(surveytemp, measure.vars = 87:92, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 87:92, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + labs(title = "Baltic herring sources")
+
  theme_gray(base_size = 24) + labs(title = "Baltic herring sources")
  
 
temp <- melt(surveytemp, measure.vars = 99:108, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 99:108, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + labs(title = "Reasons to eat Baltic herring")
+
  theme_gray(base_size = 24) + labs(title = "Reasons to eat Baltic herring")
  
 
temp <- melt(surveytemp, measure.vars = 123:132, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 123:132, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + labs(title = "Baltic herring actions")
+
  theme_gray(base_size = 24) + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + labs(title = "Baltic herring actions")
  
 
surveytemp <- subset(survey, survey[[86]] == "No")
 
surveytemp <- subset(survey, survey[[86]] == "No")
  
 
temp <- melt(surveytemp, measure.vars = 111:120, variable.name = "Variable", value.name = "Value")
 
temp <- melt(surveytemp, measure.vars = 111:120, variable.name = "Variable", value.name = "Value")
ggplot(temp, aes(x = Value, fill = Country)) + geom_histogram(position = "dodge") + facet_wrap(~ Variable)+
+
ggplot(temp, aes(x = Value, fill = Country)) + geom_bar(position = "dodge") + facet_wrap(~ Variable)+
  theme_gray(base_size = 24) + labs(title = "Reasons not to eat Baltic herring")
+
  theme_gray(base_size = 24) + labs(title = "Reasons not to eat Baltic herring")
  
 
#####################################3
 
#####################################3
Line 485: Line 488:
 
############################# Plot original data
 
############################# Plot original data
  
## Eating frequencies of fish and Baltic salmon and herring with random noise, all
+
# Eating frequencies of fish and Baltic salmon and herring with random noise, all
 
pl <- surv[surv$Eatfish,c(5,8,10,13,15)]
 
pl <- surv[surv$Eatfish,c(5,8,10,13,15)]
 
scatterplotMatrix(
 
scatterplotMatrix(
Line 571: Line 574:
 
## Kuva koko hypoteeseista
 
## Kuva koko hypoteeseista
  
traits <- as.factor(c(rep("bipedalism", 10), rep("brain", 10), rep("hairlessness", 8),
 
                      rep("fat", 4), rep("larynx", 3), rep("speech", 7), rep("other", 9)))
 
 
colstr <- c("palevioletred1","royalblue1","seagreen1","violet","khaki2","skyblue", "orange")
 
colstr <- c("palevioletred1","royalblue1","seagreen1","violet","khaki2","skyblue", "orange")
trait.cols <- colstr[traits]
 
  
hypo_sizes <- (5 - colMeans(answ))
+
#hypo_sizes <- (5 - colMeans(answ))
leg_sizes <- c(4, 3, 2, 1, 0.01)
+
#leg_sizes <- c(4, 3, 2, 1, 0.01)
  
 
#pdf(file="pcoa_plot.pdf", height=6, width=7.5)
 
#pdf(file="pcoa_plot.pdf", height=6, width=7.5)
Line 602: Line 602:
 
* Model run 13.4.2017 with first version of coordinate matrix and principal coordinate analysis [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=2k2dKhYPc2UkOCY5]
 
* Model run 13.4.2017 with first version of coordinate matrix and principal coordinate analysis [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=2k2dKhYPc2UkOCY5]
 
* Model run 20.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QXRSZOliFNfJil39] code works but needs a safety check against outliers
 
* Model run 20.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QXRSZOliFNfJil39] code works but needs a safety check against outliers
 +
* Model run 21.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=xbt2OoXPG8fi5COL] some model results plotted
 +
* Model run 21.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=waDgScuIkUvEShLy] ovariables produced by the model stored.
  
<rcode name="bayes" label="Initiate Bayes model (for developers only)" graphics=1 variables="
+
<rcode name="bayes" label="Initiate Bayes model (for developers only)" graphics=1>
name:store|description:Do you want to plot model results or store ovariables derived?|type:selection|options:
 
FALSE;Plot results;
 
TRUE;Store ovariables derived
 
">
 
 
# This is code Op_en7749/bayes on page [[Goherr: Fish consumption study]]
 
# This is code Op_en7749/bayes on page [[Goherr: Fish consumption study]]
  
Line 625: Line 623:
  
 
cat("Version with multivariate normal.\n")
 
cat("Version with multivariate normal.\n")
 +
 +
# Development needs:
 +
## Correlation between salmon.often and herring.often needs to be estimated.
 +
## Gender, country and age-spesific values should be estimated.
  
 
mod <- textConnection("
 
mod <- textConnection("
Line 688: Line 690:
 
)
 
)
  
if(!store) {
+
# Mu for all questions about salmon
  # Mu for all questions about salmon
+
scatterplotMatrix(t(js[,,1,1]))
  scatterplotMatrix(t(js[,,1,1]))
+
# All parameters for question 1 about salmon
  # All parameters for question 1 about salmon
+
scatterplotMatrix(js[1,,,1])
  scatterplotMatrix(t(js[1,,,1]))
+
# Same for herring
 +
scatterplotMatrix(js[1,,,2])
  
  jsd <- melt(js)
+
jsd <- melt(js)
  ggplot(jsd, aes(x=value, colour=Question))+geom_density()+facet_grid(Parameter ~ Fish)
+
ggplot(jsd, aes(x=value, colour=Question))+geom_density()+facet_grid(Parameter ~ Fish)
  #ggplot(as.data.frame(js), aes(x = anss.pred, y = Sampled))+geom_point()+stat_ellipse()
+
#ggplot(as.data.frame(js), aes(x = anss.pred, y = Sampled))+geom_point()+stat_ellipse()
  
  coda.j <- coda.samples(
+
coda.j <- coda.samples(
    jags,  
+
  jags,  
    c('mus', 'Omegas', 'anss.pred', 'mu0'),  
+
  c('mus', 'Omegas', 'anss.pred'),  
    1000
+
  1000
  )
+
)
  
  plot(coda.j)
+
plot(coda.j)
  
} else {
+
######## fish.param contains expected values of the distribution parameters from the model
  
  ######## fish.param contains expected values of the distribution parameters from the model
+
fish.param <- list(
 +
  mu = apply(js[,,1,], MARGIN = c(1,3), FUN = mean),
 +
  Omega = lapply(
 +
    1:2,
 +
    FUN = function(x) {
 +
      solve(apply(js[,,2:5,], MARGIN = c(1,3,4), FUN = mean)[,,x])
 +
    } # solve matrix: precision->covariace
 +
  )
 +
)
  
  fish.param <- list(
+
jsp <- Ovariable(
    mu = apply(js[,,1,], MARGIN = c(1,3), FUN = mean),
+
  "jsp",
     Omega = lapply(
+
  dependencies = data.frame(Name = "fish.param"),
 +
  formula = function(...) {
 +
     jsp <- lapply(
 
       1:2,
 
       1:2,
 
       FUN = function(x) {
 
       FUN = function(x) {
         solve(apply(js[,,2:5,], MARGIN = c(1,3,4), FUN = mean)[,,x])
+
         mvrnorm(openv$N, fish.param$mu[,x], fish.param$Omega[[x]])
       } # solve matrix: precision->covariace
+
       }
 
     )
 
     )
  )
+
   
 
+
    jsp <- rbind(
  jsp <- Ovariable(
+
      cbind(
    "jsp",
+
        Fish = "Salmon",
    dependencies = data.frame(Name = "fish.param"),
+
        Iter = 1:nrow(jsp[[1]]),
    formula = function(...) {
+
        as.data.frame(jsp[[1]])
       jsp <- lapply(
+
      ),
         1:2,
+
       cbind(
         FUN = function(x) {
+
         Fish = "Herring",
          mvrnorm(openv$N, fish.param$mu[,x], fish.param$Omega[[x]])
+
         Iter = 1:nrow(jsp[[2]]),
        }
+
        as.data.frame(jsp[[2]])
 
       )
 
       )
     
+
    )
      jsp <- rbind(
+
    jsp <- melt(jsp, id.vars = c("Iter", "Fish"), variable.name = "Question", value.name = "Result")
        cbind(
+
    jsp <- Ovariable(output=jsp, marginal = colnames(jsp) %in% c("Iter", "Fish", "Question"))
          Fish = "Salmon",
+
    return(jsp)
          Iter = 1:nrow(jsp[[1]]),
+
  }
          as.data.frame(jsp[[1]])
+
)
        ),
 
        cbind(
 
          Fish = "Herring",
 
          Iter = 1:nrow(jsp[[2]]),
 
          as.data.frame(jsp[[2]])
 
        )
 
      )
 
      jsp <- melt(jsp, id.vars = c("Iter", "Fish"), variable.name = "Question", value.name = "Result")
 
      jsp <- Ovariable(output=jsp, marginal = colnames(jsp) %in% c("Iter", "Fish", "Question"))
 
      return(jsp)
 
    }
 
  )
 
  
  # Combine modelled survey answers with estimated amounts and frequencies by:
+
# Combine modelled survey answers with estimated amounts and frequencies by:
  # Rounding the modelled result and merging that with value in ovariable assump
+
# Rounding the modelled result and merging that with value in ovariable assump
  
 
often <- Ovariable(
 
often <- Ovariable(
Line 836: Line 837:
 
)
 
)
  
  objects.store(fish.param, jsp, assump, often, much, oftenside, muchside, amount)
+
objects.store(fish.param, jsp, assump, often, much, oftenside, muchside, amount)
  cat("List fish.param and ovariables jsp, assump, often, much, oftenside, muchside, amount stored.\n")
+
cat("List fish.param and ovariables jsp, assump, often, much, oftenside, muchside, amount stored.\n")
}
 
 
</rcode>
 
</rcode>
  

Revision as of 11:45, 21 April 2017


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)

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 13.4.2017 [1]
  • Model run 20.4.2017 [2] (contains surv and helping vectors)
  • Model run 21.4.2017 [3] surv contains Eatfish, Eatherr and Eatsalm as columns.

+ Show code

Analyses

Descriptive statistics

Error creating thumbnail: Unable to save thumbnail to destination
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. [5]
  • Model run 3.3.2017. p has more dimensions. [6]
  • Model run 25.3.2017. Several model versions: strange binomial+multivarnormal, binomial, fractalised multivarnormal [7]
  • Model run 27.3.2017 [8]
  • Other models except multivariate normal were archived and removed from active code 29.3.2017.
  • Model run 29.3.2017 with raw data graphs [9]
  • Model run 29.3.2017 with salmon and herring ovariables stored [10]
  • Model run 13.4.2017 with first version of coordinate matrix and principal coordinate analysis [11]
  • Model run 20.4.2017 [12] code works but needs a safety check against outliers
  • Model run 21.4.2017 [13] some model results plotted
  • Model run 21.4.2017 [14] ovariables produced by the model stored.

+ 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


Related files

<mfanonymousfilelist></mfanonymousfilelist>

Goherr: Fish consumption study. Opasnet . [15]. Accessed 17 May 2024.