Goherr: Fish consumption study: Difference between revisions

From Opasnet
Jump to navigation Jump to search
(→‎Calculations: updated to use new Bayesian model outputs, but results are really high comapred to real questionnaire data)
Line 228: Line 228:
library(reshape2)
library(reshape2)
library(rjags)
library(rjags)
objects.latest("Op_en7748", code_name = "bayes") #: pcd.pred, ans.pred, mu.pred
pl <- melt(pcd.pred)
mul <- melt(mu.pred)
ql <- melt(ans.pred)
##Assumptions for amount and frequencies


objects.latest("Op_en7749", "preprocess")
assumptions <- opbase.data("Op_en7749.assumptions_for_calculations")
assumptions <- opbase.data("Op_en7749.assumptions_for_calculations")


for (i in 46:49) {
salmon_dishtimesassump <- Ovariable("salmon_dishtimesassump", data = assumptions[assumptions$Variable == "freq" , c("value", "Result")])
  survey[[i]] <- as.character(survey[[i]])
}


survey[,46:49][is.na(survey[,46:49])] <- ""
salmon_dishamountassump <- Ovariable("salmon_dishamountassump", data = assumptions[assumptions$Variable == "amdish" , c("value", "Result")])


for (i in 95:98) {
salmon_ingridient <- Ovariable("salmon_ingridient", data = assumptions[assumptions$Variable == "ingridient", ]["Result"])
  survey[[i]] <- as.character(survey[[i]])
}


survey[,95:98][is.na(survey[,95:98])] <- ""
salmon_sidetimesassump <- Ovariable("salmon_wholetimesassump", data = assumptions[assumptions$Variable == "freq" , c("value", "Result")])


survey$Row <- 1:nrow(survey)
salmon_sideamountassump <- Ovariable("salmon_wholeamountassump", data = assumptions[assumptions$Variable == "amside" , c("value", "Result")])
survey$Weighting <- as.double(levels(survey$Weighting))[survey$Weighting]
survey$Ages <- ifelse(as.numeric(as.character(survey$Age)) < 46, "18-45",">45")
short <- survey[c("Country", "Gender", "Row", "Ages")]


salmon_dishtimes <- Ovariable(
##Build up the variables to calculate amount of Baltic salmon eating
  "salmon_dishtimes",
  data = data.frame(
    short,
    Value = survey[[46]],
    Result = 1
  )
) #salmon with dish times / year


salmon_dishamount <- Ovariable(
salmon_dishtimes <- Ovariable("salmon_dishtimes", data = data.frame(
  "salmon_dishamount",  
                      subset(ql, ql$Question == "How often Baltic salmon" & ql$Seed == "S1"),
  data = data.frame(
                      Result = 1
    short,
))
    Value = survey[[47]],
    Result = 1
  )
) #salmon with dish amount per one serving


salmon_wholetimes <- Ovariable("salmon_wholetimes", data = data.frame(
salmon_dishtimes@data$Question <- NULL
  short,
  Value = survey[[48]],
  Result = 1
)) #salmon as such times / year


salmon_wholeamount <- Ovariable("salmon_wholeamount", data = data.frame(
salmon_dishamount <- Ovariable("salmon_dishamount", data = data.frame(
  short,
                    subset(ql, ql$Question == "How much Baltic salmon" & ql$Seed == "S1"),
  Value = survey[[49]],
                    Result = 1
  Result = 1
)) #salmon with dish amount per one serving
)) #salmon as whole amount per one serving


salmon_dishtimesassump <- Ovariable("salmon_dishtimesassump", data = assumptions[assumptions$Variable == "V46" , c("Value", "Result")])
salmon_dishamount@data$Question <- NULL 


salmon_dishamountassump <- Ovariable("salmon_dishamountassump", data = assumptions[assumptions$Variable == "V47" , c("Value", "Result")])
salmon_sidetimes <- Ovariable("salmon_sidetimes", data = data.frame(
                    subset(ql, ql$Question == "How often side Baltic salmon" & ql$Seed == "S1"),
                    Result = 1
)) #salmon as such times / year


salmon_ingridient <- Ovariable("salmon_ingridient", data = assumptions[assumptions$Variable == "salmon_ingridient", ]["Result"])
salmon_sidetimes@data$Question <- NULL
 
 
salmon_wholetimesassump <- Ovariable("salmon_wholetimesassump", data = assumptions[assumptions$Variable == "V48" , c("Value", "Result")])
salmon_sideamount <- Ovariable("salmon_sideamount", data = data.frame(
                      subset(ql, ql$Question == "How much side Baltic salmon" & ql$Seed == "S1"),
                      Result = 1
)) #salmon as whole amount per one serving


salmon_wholeamountassump <- Ovariable("salmon_wholeamountassump", data = assumptions[assumptions$Variable == "V49" , c("Value", "Result")])
salmon_sideamount@data$Question <- NULL


amount <- Ovariable("amount",  
amount <- Ovariable("amount",  
Line 296: Line 285:
                       "salmon_dishtimes",
                       "salmon_dishtimes",
                       "salmon_dishtimesassump",
                       "salmon_dishtimesassump",
                       "salmon_wholeamount",
                       "salmon_sideamount",
                       "salmon_wholeamountassump",
                       "salmon_sideamountassump",
                       "salmon_wholetimes",
                       "salmon_sidetimes",
                       "salmon_wholetimesassump"
                       "salmon_sidetimesassump"
                     )),
                     )),
                     formula = function(...) {
                     formula = function(...) {
                        
                        
                       out1 <- (salmon_dishamount * salmon_dishamountassump * salmon_ingridient)  
                       out1 <- (salmon_dishamount * salmon_dishamountassump * salmon_ingridient)  
                       out1@output$Value <- NULL
                       out1@output$value <- NULL
                       out1 <- ((salmon_dishtimes * salmon_dishtimesassump) * out1) / 365
                       out1 <- ((salmon_dishtimes * salmon_dishtimesassump) * out1) / 365
                       out1@output$Value <- NULL
                       out1@output$value <- NULL
                       out2 <- (salmon_wholeamount * salmon_wholeamountassump)
                       out2 <- (salmon_sideamount * salmon_sideamountassump)
                       out2@output$Value <- NULL
                       out2@output$value <- NULL
                       out2 <- ((salmon_wholetimes * salmon_wholetimesassump) * out2) / 365
                       out2 <- ((salmon_sidetimes * salmon_sidetimesassump) * out2) / 365
                       out2@output$Value <- NULL
                       out2@output$value <- NULL
                       out <- out1 + out2
                       out <- out1 + out2
                        
                        
Line 319: Line 308:
                     }
                     }
)
)
#amount <- EvalOutput(amount)
#### Draw N rows from the questionnaire subset filling the condition (in this case from the Country+Gender-group)
condition <- unique(survey[c("Country", "Gender")])
rows <- data.frame()
for(i in 1:nrow(condition)) {
  temp <- survey[
    survey$Gender == condition$Gender[i] & survey$Country == condition$Country[i] ,
    c("Gender", "Country", "Weighting")
    ] # Eikö tässä voisi yksinkertaistaa ja käyttää vain rivinumeroa eikä koko tempiä?
  temp <- temp[sample(1:nrow(temp), get("N", envir = openv), replace = TRUE, prob = temp$Weighting) , ]
  rows <- rbind(rows, data.frame(
    condition[i , ],
    Iter = 1:get("N", envir = openv),
    Row = as.character(floor(as.numeric(rownames(temp)))),
    Result = 1
  ))
}
# New version with same functionality
# However, the whole approach is inefficient. With N=1000 only ca 200 iterations per country and sex results from here.
# Instead, some probability distributions should be fitted based on the survey.
rows <- data.frame()
for(i in 1:nrow(condition)) {
  temp <- (1:nrow(survey))[survey$Gender == condition$Gender[i] &
                            survey$Country == condition$Country[i]]
  temp <- sample(temp, openv$N, replace = TRUE, prob = survey$Weighting[temp])
  rows <- rbind(rows, data.frame(
    Iter = 1:openv$N,
    Row = temp
  ))
}


 
amountSalmon <- EvalOutput(amount, forceEval = TRUE, N = 1000)
#rows <- Ovariable(output = rows, marginal = c(TRUE, TRUE, TRUE, TRUE, FALSE))
 
########## create ovariable for impact assessment model with the drawn individuals
 
#row <- rows@output
#row$Result <- NULL
 
salmon_dishamount@data <- merge(salmon_dishamount@data, rows) # One merge is enough to choose the rows wanted
#salmon_dishtimes@data <- merge(salmon_dishtimes@data, row)
#salmon_wholeamount@data <- merge(salmon_wholeamount@data, row)
#salmon_wholetimes@data <- merge(salmon_wholetimes@data, row)
 
amountSalmon <- EvalOutput(amount, forceEval = TRUE)


temp <- amountSalmon@output  
temp <- amountSalmon@output  
temp <- aggregate(temp$amountResult, by=list(temp$Country, temp$Ages, temp$Gender), FUN = mean)
temp <- aggregate(temp$amountResult, by=list(temp$Country, temp$Age, temp$Gender), FUN = mean)
colnames(temp) <- c("Country","Age","Gender","Salmon amount")
colnames(temp) <- c("Country","Age","Gender","Salmon amount")
oprint(temp)
oprint(temp)


#################Bayesian approach
mod <- textConnection("
  model{
    for(j in 1:4) {
      for(i in 1:N) {
        sal[i , j] ~ dbin(p[j], 9)
      }
      p[j] ~ dunif(0,1)
    }
  }
")
# This row takes data and converts it to numbers. But what should be done is to take
# each variable, remove NA, order answers meaningfully and then convert to numbers.
# A binomial distribution is assumed for bins of answer choices.
salm <- as.data.frame(lapply(survey[46:49], FUN = function(x) as.numeric(as.factor(x))))
jags <- jags.model(mod,
  data = list(N = nrow(survey),
    sal = salm),
  n.chains = 4,
  n.adapt = 100
)
update(jags, 1000)
samps <- jags.samples(jags, 'p', 1000)
plot(samps[[1]])
#### Calculate herring amount
herring_dishtimes <- Ovariable("herring_dishtimes", data = data.frame(
  short,
  Value = survey[[95]],
  Result = 1
)) #herring with dish times / year
herring_dishamount <- Ovariable("herring_dishamount", data = data.frame(
  short,
  Value = survey[[96]],
  Result = 1
)) #herring with dish amount per one serving
herring_wholetimes <- Ovariable("herring_wholetimes", data = data.frame(
  short,
  Value = survey[[97]],
  Result = 1
)) #herring as such times / year
herring_wholeamount <- Ovariable("herring_wholeamount", data = data.frame(
  short,
  Value = survey[[98]],
  Result = 1
)) #herring as whole amount per one serving
herring_dishtimesassump <- Ovariable("herring_dishtimesassump", data = assumptions[assumptions$Variable == "V95" , c("Value", "Result")])
herring_dishamountassump <- Ovariable("herring_dishamountassump", data = assumptions[assumptions$Variable == "V96" , c("Value", "Result")])
herring_ingridient <- Ovariable("herring_ingridient", data = assumptions[assumptions$Variable == "herring_ingridient", ]["Result"])
herring_wholetimesassump <- Ovariable("herring_wholetimesassump", data = assumptions[assumptions$Variable == "V97" , c("Value", "Result")])
herring_wholeamountassump <- Ovariable("herring_wholeamountassump", data = assumptions[assumptions$Variable == "V98" , c("Value", "Result")])
amount <- Ovariable("amount",
                    dependencies = data.frame(Name = c(
                      "herring_dishamount",
                      "herring_dishamountassump",
                      "herring_ingridient",
                      "herring_dishtimes",
                      "herring_dishtimesassump",
                      "herring_wholeamount",
                      "herring_wholeamountassump",
                      "herring_wholetimes",
                      "herring_wholetimesassump"
                    )),
                    formula = function(...) {
                     
                      out1 <- (herring_dishamount * herring_dishamountassump * herring_ingridient)
                      out1@output$Value <- NULL
                      out1 <- ((herring_dishtimes * herring_dishtimesassump) * out1) / 365
                      out1@output$Value <- NULL
                      out2 <- (herring_wholeamount * herring_wholeamountassump)
                      out2@output$Value <- NULL
                      out2 <- ((herring_wholetimes * herring_wholetimesassump) * out2) / 365
                      out2@output$Value <- NULL
                      out <- out1 + out2
                     
                      # Per vuosi -> per d
                     
                      return(out)
                    }
)
#### Draw N rows from the questionnaire subset filling the condition (in this case from the Country+Gender-group)
rows <- data.frame()
for(i in 1:nrow(condition)) {
  temp <- survey[
    survey$Gender == condition$Gender[i] & survey$Country == condition$Country[i] ,
    c("Gender", "Country", "Weighting")
    ] # Eikö tässä voisi yksinkertaistaa ja käyttää vain rivinumeroa eikä koko tempiä?
  temp <- temp[sample(1:nrow(temp), get("N", envir = openv), replace = TRUE, prob = temp$Weighting) , ]
  rows <- rbind(rows, data.frame(
    condition[i , ],
    Iter = 1:get("N", envir = openv),
    Row = as.character(floor(as.numeric(rownames(temp)))),
    Result = 1
  ))
}
rows <- Ovariable(output = rows, marginal = c(TRUE, TRUE, TRUE, TRUE, FALSE))
########## create ovariable for impact assessment model with the drawn individuals
row <- rows@output
row$Result <- NULL
herring_dishamount@data <- merge(herring_dishamount@data, row)
herring_dishtimes@data <- merge(herring_dishtimes@data, row)
herring_wholeamount@data <- merge(herring_wholeamount@data, row)
herring_wholetimes@data <- merge(herring_wholetimes@data, row)
amountHerring <- EvalOutput(amount, forceEval = TRUE)
temp <- amountHerring@output
temp <- aggregate(temp$amountResult, by=list(temp$Country,temp$Ages,temp$Gender), FUN = mean)
colnames(temp) <- c("Country","Age","Gender","Herring amount")
oprint(temp)


objects.store(amountHerring, amountSalmon)
#objects.store(amountHerring, amountSalmon)
cat("objects amountHerring, amountSalmon were stored.\n")
#cat("objects amountHerring, amountSalmon were stored.\n")
</rcode>
</rcode>


Revision as of 12:25, 2 March 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

Survey data will be analysed during winter 2016-2017 and results will be updated here.

+ Show code

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 will be uploaded to Opasnet base on Octobere 2016.

The R-code to analyse the survey data will be provided on this page later on.

Data

Original datafile File:Goherr fish consumption.csv

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.

+ Show code

Calculations

This code calculates how much (g/day) Baltic herring and salmon are eaten based on the questionnaire data. Calculation takes several minutes!

+ Show code

Assumptions

The following assumptions are used:

Assumptions for calculations(-)
ObsVariablevalueExplanationResult
1freq6times per year260 - 364
2freq5times per year104 - 208
3freq4times per year52
4freq3times per year12 - 36
5freq2times per year2 - 5
6freq1times per year0.5 - 0.9
7freq0times per year0
8amdish0grams / serving20 - 50
9amdish1grams / serving70 - 100
10amdish2grams / serving120 - 150
11amdish3grams / serving170 - 200
12amdish4grams / serving220 - 250
13amdish5grams / serving270 - 300
14amdish6grams / serving450 - 500
15ingridientfraction0.1 - 0.3
16amside0grams / serving20 - 50
17amside1grams / serving70 - 100
18amside2grams / serving120 - 150
19amside3grams / serving170 - 200
20amside4grams / serving220 - 250

Questionnaire


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.

Formula

See also

Keywords

References


Related files

<mfanonymousfilelist></mfanonymousfilelist>

Goherr: Fish consumption study. Opasnet . [1]. Accessed 27 Nov 2024.


Retrieved from "https://en.opasnet.org/index.php?title=Goherr:_Fish_consumption_study&oldid=40124"