setwd("D://PhD Project/3rd Analysis/Data") getwd() #Install Packages install.packages("support.BWS") install.packages("crossdes") install.packages("mlogit") install.packages("gmnl") install.packages("ggplot2") #Load packages library("support.BWS") library("crossdes") library("mlogit") library("gmnl") #to estimate more kinds of model including random parameter model library("ggplot2") #Set No of attributes to appear in each question options(digit=4) ### Designing a choice situation # Setting variables for Case 1 BWS ### Generating a BIBD # Setting the seed in order to generate the same BIBD in the future set.seed(8020) des <- find.BIB( trt = 13, b = 13, k = 4) des # Checking whether the resultant design is a BIBD isGYD(des) ### Creating Questions # Storing the names of the items (attributes) items.meat <- c("Origin","Price","Breed","Brand","Color", "Taste","Tenderness","Halal_Label","Freshness", "Cleanliness","Animal_Welfare","Food_Safety_Certification", "Quality_Appearance") items.meat # Obtaining a series of BWS questions bws.questionnaire(choice.sets = des, design.type = 2, item.names = items.meat) ########## Conducting the Survey # Creating a Data set # Importing dataset meatbws0 <- read.csv("meatbws2.csv") View(meatbws0) # Have a look print(meatbws0) ########### Check best worst answers for bad answers ############ # to check the best worst answer. All of the combination of best-worst answers should be different. table(meatbws0$b1 == meatbws0$w1) table(meatbws0$b2 == meatbws0$w2) table(meatbws0$b3 == meatbws0$w3) table(meatbws0$b4 == meatbws0$w4) # include false answers to BWS choice set 4 // perhaps selecting one attribute as both, best and worst table(meatbws0$b5 == meatbws0$w5) table(meatbws0$b6 == meatbws0$w6) table(meatbws0$b7 == meatbws0$w7) table(meatbws0$b8 == meatbws0$w8) table(meatbws0$b9 == meatbws0$w9) table(meatbws0$b10 == meatbws0$w10) table(meatbws0$b11 == meatbws0$w11) table(meatbws0$b12 == meatbws0$w12) table(meatbws0$b13 == meatbws0$w13) # include false answer to BWS choice set 13 // perhaps selecting one attribute as both, best and worst meatbws1 <- meatbws0[(meatbws0$b4 != meatbws0$w4)&(meatbws0$b13 != meatbws0$w13),] # to remove bad answers View(meatbws1) ## New Dataset after removing respondents with incomplete answers or false answers to BWS choice sets ## checking if there are still bad answers table(meatbws1$b1 == meatbws1$w1) table(meatbws1$b2 == meatbws1$w2) table(meatbws1$b3 == meatbws1$w3) table(meatbws1$b4 == meatbws1$w4) table(meatbws1$b5 == meatbws1$w5) table(meatbws1$b6 == meatbws1$w6) table(meatbws1$b7 == meatbws1$w7) table(meatbws1$b8 == meatbws1$w8) table(meatbws1$b9 == meatbws1$w9) table(meatbws1$b10 == meatbws1$w10) table(meatbws1$b11 == meatbws1$w11) table(meatbws1$b12 == meatbws1$w12) table(meatbws1$b13 == meatbws1$w13) ## No bad answers in the new dataset called meatbws1 ########### Creating dummy variables for categorical Socio-demographic variables ######### ## GENDER meatbws1$female=ifelse(meatbws1$gend==2, 1, 0) # creating gender variable for female only table(meatbws1$female) ## AGE meatbws1$age1=ifelse(meatbws1$age==1, 1, 0) #young-aged ## REFERENCE table(meatbws1$age1) meatbws1$age2=ifelse(meatbws1$age>=2, 1, 0) #middle-aged table(meatbws1$age2) ## EDUCATION meatbws1$educ3=ifelse(meatbws1$educ==6, 1, 0) #College Level table(meatbws1$educ3) meatbws1$educ4=ifelse(meatbws1$educ==7, 1, 0) #Bachelor's degree table(meatbws1$educ4) meatbws1$educ5=ifelse(meatbws1$educ==8, 1, 0) #Graduate degree table(meatbws1$educ5) meatbws1$educ0=ifelse(meatbws1$educ>= 4 & meatbws1$educ<= 5, 1, 0) #No_degree ##REFERENCE table(meatbws1$educ0) ## Household INCOME meatbws1$Inc_low=ifelse(meatbws1$inc<4, 1, 0) # Low_income table(meatbws1$Inc_low) meatbws1$Inc_mid=ifelse(meatbws1$inc>=4 & meatbws1$inc<6, 1, 0) # Middle_income table(meatbws1$Inc_mid) meatbws1$Inc_high=ifelse(meatbws1$inc>=6 & meatbws1$inc<11, 1, 0) # High_income ##REFERENCE table(meatbws1$Inc_high) meatbws1$undisclosed=ifelse(meatbws1$inc==11, 1, 0) # undisclosed table(meatbws1$undisclosed) ## Meat consumption frequency #Habitual consumers: between 1 ~ 5 times a week meatbws1$Habitual=ifelse(meatbws1$meat_cons<3, 1, 0) ##REFERENCE table(meatbws1$Habitual) #Occasional consumers meatbws1$Occasional=ifelse(meatbws1$meat_cons>=3 & meatbws1$meat_cons<6, 1, 0) table(meatbws1$Occasional) #Never eat meat meatbws1$Never_eat=ifelse(meatbws1$meat_cons==6, 1, 0) table(meatbws1$Never_eat) ##Risk Preference ## Risk_Pref will be incorporated as continuous variable describe(meatbws1$Risk_Pref) ## In the data, Risk preference range is from 1 - 11, whereas in the questionnaire it is from 0 - 10. #The Questpro data collection software didn't accept 0 as an answer to the first value of the range. ########### Creating a dataset ######### #3.8.1:Preparing a dataset including responses to the questions ## My data is already stored as "meatbws1" dim(meatbws1) #Displays number of responses and columns names(meatbws1) # names of columns ## Combining the design and the responses # create the dataset for analysis using the command bws.dataset() # we store the names of the response variables used in the respondent dataset meatbws1 in the vector response.meat response.meat <- colnames(meatbws1)[2:27] response.meat ## Note: the respondents answer the BWS questions on the basis of the maxdiff model # Dataset for the maxdiff model meat.data <- bws.dataset( data = meatbws1, # data frame containing response dataset response.type = 1, # format of response variables: 1 = row number format choice.sets = des, # choice sets design.type = 2, # 2 if a BIBD is assigned to choice.sets item.names = items.meat, # the names of items id = "id", # the name of respondent id variable response = response.meat, # the names of response variables model = "maxdiff") # the type of dataset created by the function is the maxdiff model dim(meat.data) ##OUTput: 35412 rows (=227 repondents X 13 BWS questions X 13 possible pairs BWS question) and 123 columns (variables) meat.data[1:12, ] ## to see responses of the first 12 subjects ######## Measuring Preferences based on counting approach ############# ## The counting approach cs <- bws.count(meat.data, cl = 2) dim(cs) ## OUTPUT: 227 ROWS (respondents) 97 columns (variables) names(cs) ## names of the 109 columns (variables) ## To understand the relative importance of items among individuals, we should note the means of the BWS scores of items as well as the variances (standard deviations). ## The graphical presentation of individual BWS scores allows us to easily read the relationship between the means and standard deviations of BWS scores. plot( x = cs, score = "bw", # BW scores are used pos = 4, # Labels are added to the right of the points ylim = c(1, 2.5), # the y axis ranges from 1 to 2.5 xlim =c(-3, 3)) # the x axis ranges from -3 to 3 ## Heterogeneity in detail using bar plot of BWS scores par(mar = c(5, 4, 2, 1)) barplot( score = "bw", # BW scores are used height = cs, mfrow = c(4, 4)) # Bar plots are drawn in a 4-by-4 array ###############New Heterogeneity bar plot // Figure 2 ## The counting approach #Install tidyverse package install.packages("tidyverse") # Load libraries library(tidyverse) # loads ggplot2, dplyr, tidyr # 1. Extract only BW score columns bw_data <- cs[, attributes(cs)$bw.names] # 2. Add respondent ID bw_data$Respondent <- 1:nrow(bw_data) # 3. Reshape to long format bw_long <- bw_data %>% pivot_longer( cols = -Respondent, names_to = "Attribute", values_to = "Score" ) # 4. Clean attribute names bw_long <- bw_long %>% mutate(Attribute = gsub("^bw\\.", "", Attribute), Attribute = str_to_title(gsub("\\.", " ", Attribute))) # Capitalize + replace dots with spaces #5. Overlay Histogram plus Density with bold text ggplot(bw_long, aes(x = Score)) + geom_histogram(aes(y = ..density..), binwidth = 1, fill = "gray85", color = "black", boundary = 0) + geom_density(color = "blue", fill = "turquoise", alpha = 0.5) + facet_wrap(~ Attribute, scales = "fixed", ncol = 4) + theme_minimal() + labs(x = "Score", y = "Density") + theme( strip.text = element_text(size = 12, face = "bold"), # Facet labels axis.title.x = element_text(size = 14, face = "bold"), # X-axis title axis.title.y = element_text(size = 14, face = "bold"), # Y-axis title axis.text = element_text(size = 12), # Axis tick labels plot.title = element_text(size = 16, face = "bold") # (optional) Plot title ) ################ FIGURE 4 ## Drawing bar plots for mean standardized BWS scores par(mar = c(5,11,2,1)) barplot(score = "bw", height = cs, mean = TRUE, error.bar = "ci", las = 1, names.arg = c("Color","Animal breed","Brand","Animal welfare","Tenderness", "Quality appearance","Price","Origin","Cleanliness", "Taste","Halal label","Freshness", "Food safety certification"), col=c("darkred", "firebrick", "red", "salmon", "lightsalmon", "lightpink", "mistyrose", "lightblue", "skyblue", "cornflowerblue", "royalblue", "blue", "darkblue")) ##### FURTHER ANALYSIS OF MEAT ATTRIBUTES PREFERENCES ## Summary table for the counting approach # The summary() method calculates (1) aggregated B, W, BW, and standardized BW scores, #(2) item ranks based on the BW score, #(3) means of B, W, BW, and standardized BW scores, and #(4) the square root of the ratio of aggregated B to aggregated W scores and the corresponding standardized scores summary(cs) sum(cs) ############ Modeling approach for the Analysis of BWS questions ######### ########### Measuring Preferences based on conditional logit model using mlogit function # an arbitrary item, in this case, the least important one, which is color, has been removed ## Fitting the conditional logit model on the basis of maxdiff model #systematic component of the utility function mf <- RES ~ Origin + Price + Breed + Brand + Taste + Tenderness + Halal_Label + Freshness + Cleanliness + Animal_Welfare + Food_Safety_Certification + Quality_Appearance - 1 ## The dataset created by the bws.dataset() function has to be converted into that of the S3 class "dfidx" # with the dfidx() function in the dfidx package library(dfidx) # Data set for the maxdiff model / prepare dataset for conditional logit model meat.data.dfidx <- dfidx(data = meat.data, idx = list(c("STR", "id"), "PAIR"), choice = "RES") ## After preparing the dataset for the Conditional Logit Model (CLM) analysis, # fit the maxdiff model using mlogit() with the model formula mf and the dataset meat.data.dfidx: meat.out <- mlogit(formula = mf, data = meat.data.dfidx) summary(meat.out) ## To see the extent of relative importance among the thirteen characteristics, ## the share of preferences for them are calculated using bws.sp(): sp.meat <- bws.sp(meat.out, base = "color") sp.meat ## Calculate the goodness of fit measure # The value of McFadden’s R-squared for our model and its adjusted version are calculated as bellow: LL0 <- - 227 * 13 * log(12) # the value of log-likelihood at zero LLb <- as.numeric(meat.out$logLik) # the value of log-likelihood at convergence 1 - (LLb/LL0) # McFadden's R-squared #result: 0.1146 1 - ((LLb - 6)/LL0) # adjusted McFadden's R-squared #result" 0.1138 ######### Random Parameters Logit Model ########### ## the model does not show a good fit, because the heterogeneity is not considered in the model. ## We apply random parameters logit model to the data # Assuming that the coefficients of the item variables are normally distributed, Halton sequences are used in the simulation process of the random parameters logit model estimation. rpl1 <- mlogit( formula = RES ~ Origin + Price + Breed + Brand + Taste + Tenderness + Halal_Label + Freshness + Cleanliness + Animal_Welfare + Food_Safety_Certification + Quality_Appearance - 1| 0, data = meat.data.dfidx, rpar = c(Origin = "n", Price = "n", Breed = "n", Brand = "n", Taste = "n", Tenderness = "n", Halal_Label = "n", Freshness = "n", Cleanliness = "n", Animal_Welfare = "n", Food_Safety_Certification = "n", Quality_Appearance = "n"), R = 100, halton = NA, panel = TRUE) summary(rpl1) ##Table 4, Random parameters logit model estimates ## The value of McFadden's R-squared and its adjusted version for the fitted model are given as follows: LLb.rpl <- as.numeric(rpl1$logLik) 1 - (LLb.rpl/LL0) # [1] 0.154 1 - ((LLb.rpl - 12)/LL0) # [1] 0.153 ## the goodness-of-fit measures for the random parameters logit model are improved compared with those for the conditional logit model. ## Calculate the shares of preferences for the thirteen items using the estimated mean coefficients, from origin to qualityappearance, are as follows: bws.sp(object = rpl1, base = "color", coef = c("Origin", "Price", "Breed", "Brand", "Taste", "Tenderness", "Halal_Label", "Freshness", "Cleanliness", "Animal_Welfare", "Food_Safety_Certification", "Quality_Appearance")) ## Table 4, Share of Preferences ## When fitting the random parameters logit model, the output from mlogit() includes individual-specific parameters in the object indpar: head(rpl1$indpar) ## We did not report this in the paper ### NOTE: the above RPL model produced some standard deviations with negative signs. # However, the negative signs are not a problem — they are estimation artifacts due to unconstrained optimization in mlogit. # The true standard deviations are their absolute values. ##########Uncertainty of Preference Share Estimates - To calculate Standard Error of Share of Preferences ######## ##18.4.2025 ## Uncertainty of Preference Share estimates library(MASS) # 1. Extract coefficients and variance-covariance matrix coef_estimates <- coef(rpl1) vcov_matrix <- vcov(rpl1) # 2. Define the order of coefficients (same order as in bws.sp) coef_names <- c("Origin", "Price", "Breed", "Brand", "Taste", "Tenderness", "Halal_Label", "Freshness", "Cleanliness", "Animal_Welfare", "Food_Safety_Certification", "Quality_Appearance") # 3. Simulate random draws of coefficients set.seed(123) # for reproducibility simulated_betas <- mvrnorm( n = 1000, mu = coef_estimates[coef_names], Sigma = vcov_matrix[coef_names, coef_names] ) # 4. Add a column for Color (fixed utility = 0) simulated_betas_with_color <- cbind(Color = 0, simulated_betas) # 5. Calculate preference shares exp_betas <- exp(simulated_betas_with_color) simulated_shares <- exp_betas / rowSums(exp_betas) # 6. Compute Mean Share and Standard Error mean_shares <- colMeans(simulated_shares) se_shares <- apply(simulated_shares, 2, sd) # 7. Prepare a clean results table preference_shares_results <- data.frame( Attribute = c("Color", coef_names), MeanShare = mean_shares, SE = se_shares ) # 8. View the result print(preference_shares_results) ## Table 4, Share of Preferences with Standard Errors in Brackets #########Random parameters logit model with gmnl() package ######### library(gmnl) meat.data.ml <- mlogit.data(data =meat.data, choice = "RES", shape = "long", alt.var = "PAIR", chid.var = "STR", id.var = "id") summary(meat.data.ml) ##### Random Parameter model with subscales of MAQ ##### ##8.2.25## # Mixed Logit Model incorporating multiple individual characteristics #Model1 meatFull.rpl <- gmnl(RES ~ origin + price + breed + brand + taste + tenderness + label + freshness + cleanliness + animalwelfare + foodsafety + appearance | 0 | 0 | female + age2 + educ3 + educ4 + educ5 + Inc_mid + Inc_high + undisclosed + Occasional + Never_eat + Risk_Pref - 1, data = meat.data.ml, ranp = c(origin = "n", price = "n", breed = "n", brand = "n", taste = "n", tenderness = "n", label = "n", freshness = "n", cleanliness = "n", animalwelfare = "n", foodsafety = "n", appearance = "n"), mvar = list(origin = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), price = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), breed = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), brand = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), taste = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), tenderness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), label = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), freshness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), cleanliness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), animalwelfare = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), foodsafety = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), appearance = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref")), model = "mixl", R = 100, halton = NA, panel = TRUE) ## with 100 halton draws # Display results summary(meatFull.rpl) ##Model2 meatFull2.rpl <- gmnl(RES ~ origin + price + breed + brand + taste + tenderness + label + freshness + cleanliness + animalwelfare + foodsafety + appearance | 0 | 0 | female + age2 + educ3 + educ4 + educ5 + Inc_mid + Inc_high + undisclosed + Occasional + Never_eat + Risk_Pref - 1, data = meat.data.ml, ranp = c(origin = "n", price = "n", breed = "n", brand = "n", taste = "n", tenderness = "n", label = "n", freshness = "n", cleanliness = "n", animalwelfare = "n", foodsafety = "n", appearance = "n"), mvar = list(origin = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), price = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), breed = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), brand = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), taste = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), tenderness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), label = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), freshness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), cleanliness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), animalwelfare = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), foodsafety = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref"), appearance = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref")), model = "mixl", R = 500, halton = NA, panel = TRUE) ## with 500 halton draws # Display results summary(meatFull2.rpl) ##Interacting MAQ into RPL model meatFull3.rpl <- gmnl(RES ~ origin + price + breed + brand + taste + tenderness + label + freshness + cleanliness + animalwelfare + foodsafety + appearance | 0 | 0 | hedonism + affinity + dependence + entitlement - 1, data = meat.data.ml, ranp = c(origin = "n", price = "n", breed = "n", brand = "n", taste = "n", tenderness = "n", label = "n", freshness = "n", cleanliness = "n", animalwelfare = "n", foodsafety = "n", appearance = "n"), mvar = list(origin = c("hedonism", "affinity", "dependence", "entitlement"), price = c("hedonism", "affinity", "dependence", "entitlement"), breed = c("hedonism", "affinity", "dependence", "entitlement"), brand = c("hedonism", "affinity", "dependence", "entitlement"), taste = c("hedonism", "affinity", "dependence", "entitlement"), tenderness = c("hedonism", "affinity", "dependence", "entitlement"), label = c("hedonism", "affinity", "dependence", "entitlement"), freshness = c("hedonism", "affinity", "dependence", "entitlement"), cleanliness = c("hedonism", "affinity", "dependence", "entitlement"), animalwelfare = c("hedonism", "affinity", "dependence", "entitlement"), foodsafety = c("hedonism", "affinity", "dependence", "entitlement"), appearance = c("hedonism", "affinity", "dependence", "entitlement")), model = "mixl", R = 500, halton = NA, panel = TRUE) # Display results summary(meatFull3.rpl) ################# ##New RPL Model with all covariates in a single Table ##Model2 meatFull2.rpl <- gmnl(RES ~ Origin + Price + Breed + Brand + Taste + Tenderness + HalalLabel + Freshness + Cleanliness + AnimalWelfare + FoodSafetyCertification + QualityAppearance | 0 | 0 | female + age2 + educ3 + educ4 + educ5 + Inc_mid + Inc_high + undisclosed + Occasional + Never_eat + Risk_Pref + hedonism + affinity + dependence + entitlement - 1, data = meat.data.ml, ranp = c(Origin = "n", Price = "n", Breed = "n", Brand = "n", Taste = "n", Tenderness = "n", HalalLabel = "n", Freshness = "n", Cleanliness = "n", AnimalWelfare = "n", FoodSafetyCertification = "n", QualityAppearance = "n"), mvar = list(Origin = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Price = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Breed = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Brand = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Taste = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Tenderness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), HalalLabel = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Freshness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), Cleanliness = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), AnimalWelfare = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), FoodSafetyCertification = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement"), QualityAppearance = c("female", "age2", "educ3", "educ4", "educ5", "Inc_high", "Inc_mid", "undisclosed", "Occasional", "Never_eat", "Risk_Pref", "hedonism", "affinity", "dependence", "entitlement")), model = "mixl", R = 500, halton = NA, panel = TRUE) ## with 500 halton draws # Display results summary(meatFull2.rpl) warnings()