[**NOTE**: I’m not finished writing this thing yet — there’s a lot more that could be done, and sadly many other projects currently demand both my CPU and attention. I’m cutting a few corners in the interests of time and convenience (e.g. it’s a complete-cases analysis that assumes data is MCAR, doesn’t take into account measurement uncertainty, doesn’t incorporate multivariate outcomes, does the bare minimum w.r.t. diagnosing MCMC health and performance, exploring results beyond superficially, and examining and averaging over a fairly limited number of plausible models). And as a standard caveat, I *might* have also made a nontrivial mistake in my code or reasoning somewhere; though I did check things over fairly carefully, to err is human. Regardless, let this post serve as a placeholder for now. Also, WordPress likes making pictures tiny unless your browser’s really zoomed in, so if you want to see a figure in full-res, open it in its own window and change “*w=*” to some large value]

[**NOTE2: **I think I have an idea for implementing the models I’d like to proceed, both in the simple case and the complex case (some of which I’ve done and need to write up). Still, since this was a side project that nobody but me really cares about (I’ve gotten some… rather explicit statements of disinterest from relevant parties lol — I think people are maybe just opposed to Bayesian inference?), it’s hard to motivate myself to do the boring stuff when I could be doing *other* boring stuff. I’ll get back to it someday!

So a little while ago the organization “Mercy for Animals” (MFA) ran a “survey experiment” to see how effective online veg*n outreach was at getting people to change their dietary habits and/or their minds about animal welfare. Since I know a wee bit of Bayesian Regression (not much, but enough to be dangerous!) and donate a decent amount of money to Mercy for Animals (and other ACE-recommended charities, as well as the HSA), and *also* since I recently had some free time (spring break woooooo!), I decided to play around with the data produced by MFA’s experiment. My step-by-step analysis and discussion can be seen below, and you can see MFA’s write-up and description, as well as download all the original raw data, by clicking on the following link, with a more thorough description of the experiment here.

To summarize what I did in **R** and **Stan**: after recoding/cleaning some of the data (shortening sentences to phrases/words, standardizing predictors, reordering levels, basic housekeeping stuff), I did some quick Exploratory Data Analysis (EDA), fit several multilevel generalized linear mixed models (GLMMs) that I thought reasonable (bordering on dredging, though, since I find lotsa stuff reasonable) with outcomes including the ordered categorical response to the Likert-scale attitude questions (using a cumulative link function – modeling outcomes with log cumulative odds and an ordered logit likelihood), and the integer “Servings of X in the past two days” questions (with a zero-inflated Poisson mixed model, since outcomes are non-negative integers, rates/mean are low, and people can forego different animal products both because they have a veg*n diet and because they just didn’t happen to have any over two days). Plausible, weakly informative, regularizing priors were used, sampling from the joint posterior was done using Hamiltonian (Markov Chain) Monte Carlo (HMC), and results were interpreted using posterior predictive and counterfactual simulation while averaging across models with WAIC (the “Widely Applicable Information Criterion”) weights. I’d wanted to impute missing data and accommodate measurement uncertainty, but Stan doesn’t seem to permit discrete variables as stochastic nodes, and almost all the data is discrete by design (Likert attitudes, # of servings, “more”-“less”-“same amount”, etc.), and the most continuous variable is not one that’s likely to be too mis-measured (age). As a result, I unfortunately had to throw a lot (a quarter-ish) of the data away (and assume that it’s MCAR); for anyone hoping to draw hard and fast conclusions from any of this, you’d really need to impute it properly, as excluding measurements could easily introduces bias.

I’d also initially wanted to fit models with varying slopes, but they were taking too dang long to run, so I restricted myself to models with varying intercepts (aka adaptive regularizing priors, aka hyperpriors). The handful I did fit for pork consumption did pretty poorly in model comparison (receiving 0 WAIC weight), but that doesn’t mean they’d be overfit for any of the other outcomes. A more complete analysis should definitely aim to include them.

Some benefits of this approach (imo) include allowing me to accommodate uncertainty in model selection, parameter estimation, measurement (by modeling data as distributions – though I didn’t do this here, I’d like to when I revisit this material later), and prediction (I get to use the whole posterior for predictive simulation, instead of point or interval estimates. Simulation also allows for much more intuitive interpretation of results compared to trying to make sense of huge tables of partially represented estimates), avoids overfitting by using using (adaptive) regularizing priors and WAIC to estimate out-of-sample deviance and average across model uncertainty, statistically “account” for stuff like masking, model the data-generating process a lot more realistically (e.g. by preserving ordinality in stuff like “Intended Future Meat Consumption” and ordering in age, instead of having to lose a ton of resolution by binning), and so on. Doing multiple regression lets me more capably account for and reduce “noise” in simulated outcomes by allowing me to condition on different combinations of predictors (e.g. relating to gender, age, country, etc.) and more clearly see the effect of the variable I’m actually interested in (i.e. exposure to the experimental treatment).

All in all, this was a fun exercise and let me refresh some rusty data analysis skills! My code is included below and commented/annotated from start to finish, so feel free to run it yourself (maybe only fit those models you’ve special interest in – collectively, they all still took many nights to fit on my moderately overclocked i5-4690K, even distributing the workload across multiple cores). Non-centered parameterization could help with that, I think. **Any comments, criticisms, and suggestions are welcome!**

So to start, lets boot up R and load up some of the packages we’ll be using:

#load packages library(rethinking) library(mcmcplots) library(ggplot2) library(tidyr) library(dplyr) library(gridExtra)

Make sure you have Stan installed and configured too. All of the model fitting below is done through Richard’s R package “**rethinking**” (McElreath, R. (2016). Rethinking: Statistical Rethinking book package, version 1.58. GitHub project: rmcelreath/rethinking). Richard’s a brilliant and fun dude; thanks for putting *rethinking* together!

Now load the data. I cleaned it up a little in Excel, since the raw data were a bit messy and scattered across multiple files. You can find the csv file with the combined data here. Save it somewhere, load it into R, and examine it.

d <- read.csv("C:\\Users\\Nikolai\\Dropbox\\MFA_Survey\\CombinedData.csv") str(d)

What do all these variables represent?

**Variable Name: **Description of Variable

**Treat**: Whether the individual was exposed to the 12-minute pro-veg*n advertisement (“What Cody Saw Will Change Your Life“) two to four months prior; a 1 indicates that they were subject to this “experimental treatment” and a 0 indicated they were subject to the “control” (an unrelated video on fighting tropical diseases). Note that a 1 does not mean that they actually watched the video, but rather that they clicked on a link that took them to the video, which they might have immediately closed, or closed after a minute, or whatever.

**Platform**: What website the survey was administered on

**Age**: The age of the person responding to the survey; it doesn’t look like ages below 12 were permissible

**X2Days…**: The person’s response to the question: In the __past two days__, how many servings have you had of the following foods? Please give your best guess.

**X2DaysPork**: Pork (ham, bacon, ribs, etc.)

**X2DaysBeef**: Beef (hamburgers, meatballs, in tacos, etc.)

**X2DaysDairy**: Dairy (milk, yogurt, cheese, etc.)

**X2DaysEggs**: Eggs (omelet, in salad, etc.)

**X2DaysBird**: Chicken and Turkey (fried chicken, turkey sandwich, in soup, etc.)

**X2DaysSea**: Fish and Seafood (tuna, crab, baked fish, etc.)

**P4MonIncr**: Response to the question: “At any time in the past four months, have you increased or decreased the total amount of meat, chicken and fish you ate? This may be a choice you made for a few days or for more than a few days.” An answer in this column indicates that the individual selected “I increased the total amount of meat, chicken and fish I ate”.

**P4MonNC**: Response to the question: “At any time in the past four months, have you increased or decreased the total amount of meat, chicken and fish you ate? This may be a choice you made for a few days or for more than a few days.” An answer in this column indicates that the individual selected “I did not change the total amount of meat, chicken and fish I ate”.

**P4MonDecr**: Response to the question: “At any time in the past four months, have you increased or decreased the total amount of meat, chicken and fish you ate? This may be a choice you made for a few days or for more than a few days.” An answer in this column indicates that the individual selected “I decreased the total amount of meat, chicken and fish I ate”.

**NoMeatAlready**: Response to the question: “At any time in the past four months, have you increased or decreased the total amount of meat, chicken and fish you ate? This may be a choice you made for a few days or for more than a few days.” An answer in this column indicates that the individual selected “I did not eat meat, chicken, and fish then and I do not eat meat, chicken, and fish now”.

**NowVs4MonPast**: Response to the question: “Compared to four months ago, are you eating more, less, or about the same amount of meat, chicken and fish (in total)?” Responses were one of “I am eating (*more/about the same amount of/less*) meat, chicken and fish (in total) as I was four months ago”.

**NowVs4MonFut**: Response to the question: “Four months from now, how much meat, chicken and fish (in total) do you think you will be eating?” Responses were one of “I will probably be eating (*more/less/about the same amount of/no*) meat, chicken and fish (in total) as I am now.

**MeatUniqueSmart**: Response were on a typical 5-point Likert scale (*Strongly Disagree, Somewhat Disagree, Neither Agree nor Disagree, Somewhat Agree, Strongly Agree*) to the following statement: “Cows, pigs and chickens are intelligent and emotional with unique personalities.”

**LessMeatGood**: Response were on a typical 5-point Likert scale to the following statement: “Eating less meat, chicken and fish (in total) is the right thing to do.”

**MeatReplace**: Response were on a typical 5-point Likert scale to the following statement: “I know how to replace meat, chicken and fish dishes with appealing non-meat options.”

**Gender**: What is the responder’s gender? Options given were (*male/female/other*).

**Country**: What country the responder lives in.

There’s still come tidying we need to do before getting started, so run the following code. See comments in code or text for justification or explanation:

#let's reorder some levels d$NowVs4MonPast <- factor(d$NowVs4MonPast, levels = c("less", "same", "more")) d$NowVs4MonFut <- factor(d$NowVs4MonFut, levels = c("none", "less", "same", "more"))

Ah, now I think this is a mite goofy. Why have a *none* in there? *None* distinguishes an absolute amount (i.e., 0), but all the other ones are relative amounts! Someone eating none can be eating less or same, as well, so how do they know what to pick? And the NowVs4MonPast variable doesn’t have a “none” option? There *is* a “none” option describing consumption over the past four months for a different question (answers under “NoMeatAlready”), but answers there don’t entirely match up to answers for “NowVs4MonPast”, as we shall see later.

d$MeatUniqueSmart <- factor(d$MeatUniqueSmart, levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")) d$LessMeatGood <- factor(d$LessMeatGood, levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")) d$MeatReplace <- factor(d$MeatReplace, levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree"))

#now let's look at and get a sense of the data d$anyChangeP4M <- as.factor(paste0(d$P4MonIncr, d$P4MonNC, d$P4MonDecr)) summary(d$anyChangeP4M)

Since this asks about any change in the past 4 months, we get some rather flip-floppy answers, where some people increased, decreased, and did not change their behaviors, all at the same time (or over the past 4 months). I dunno that I really want to look at this variable, especially since something like “NowVs4MonPast” has a much more clear-cut interpretation, and the number of servings variables clearer interpretation still. Here, we don’t know if it’s a lasting or super ephemeral change, or if they did it purely by chance or for unrelated reasons (like going backpacking or out for fieldwork or something — I know my diet can change rather dramatically on the road).

summary(d$Gender)

……….Female Male Other Other_NS

372 1480 160 13 1

There’s an “Other_NS” option; I dunno how it differs from “Other”, but I’ll just recode it as “Other” for now.

d$Gender[which(d$Gender == "Other_NS")] <- "Other" d$Gender <- factor(d$Gender)

summary(d$Age)

#Why is there a "12 and under" category? Age is a perfectly good numeric variable! d$c.Age <- as.character(d$Age) d$c.Age[which(d$c.Age == "12 or younger" | d$c.Age == "12")] <- "12" d$c.Age <- as.numeric(d$c.Age) summary(d$c.Age)

Min. 1st Qu. Median Mean 3rd Qu. Max.

12.00 15.00 17.00 17.57 20.00 25.00

Can you not ask people younger than 12 how old they are? Don’t you need to be at least 13 to be on facebook anyway? Whatever, there are only 85 of them, for now I’ll recode the “12 and under”s as 12, which seems like a reasonable assumption, as there aren’t going to be 5 year olds taking surveys on the interner. I want age to be numeric and not categorical. I can later try coding it as missing or accommodating fuzziness in age assignment statistically (by modeling) or dumping all the data from 12 and unders, since it might have resulted from a different person taking the survey than the one who watched the video from facebook (in the case of multiple individuals using the same computer). If it doesn’t take too long, I can later rerun the analysis under each approach to assess sensitivity to this decision. Let’s also standardize this variable to ease later interpretation and improve MCMC performance:

d$Age.s <- (d$c.Age - mean(d$c.Age)) / sd(d$c.Age) #also called "z-scores"

So now that we’ve done all of that, we can continue to visualizing different aspects of the data.

#not looking at NAs for now consumption <- d[,c("Treat", "X2DaysPork", "X2DaysBeef", "X2DaysDairy", "X2DaysEggs", "X2DaysBird", "X2DaysSea")] %>% gather(Type, Amount, -Treat) ggplot(data = consumption, aes(x=Type, y=Amount, color=Treat)) + geom_jitter() + theme_bw() #do note -- this removes over half the data, which we can try to impute later

Hmm, interesting. I guess when the survey asked “In the past two days, how many servings have you had of Port, Beef, Dairy, Eggs, Chicken, Fish?” it didn’t allow responses above 9? I don’t think it coded them as NAs, since all the variables have 228 NAs. This might be failing to capture some of the effect (I’ve known a lot of people to eat several dozen servings of dairy, beef, eggs, chicken, etc. over 2 days), though maybe young women on facebook are less likely to do so than 250 lb male powerlifters. One quick informal test could be to look for an excess of 9s relative to 8s, 7s, 6s, etc., since people who ate more than 9 servings would just put 9 absent a 10, 11, 12, etc. option. It would be easier to visualize this as a series of histograms, rather than as jittered scatterplots. We can do this for the entire dataset using:

par(mfrow = c(2,3)) hist(d$X2DaysPork); hist(d$X2DaysBeef); hist(d$X2DaysDairy); hist(d$X2DaysEggs); hist(d$X2DaysBird); hist(d$X2DaysSea)

In these plots, I see more 9s than I’d expect if the data were actually ~ZIP distributed (see below), since outcomes of a Poisson process should dwindle in their frequency as you get to higher and higher values (on average/for higher samples, ignoring simulation variance). An alternative explanation could be that a Poisson process does not adequately approximate eating behavior, and that people just naturally gravitate towards eating 9 servings every 2 days, such that our distributions of observed outcomes are and should indeed be mildly bimodal. That seems somewhat implausible, though. And of course, observed outcomes needn’t *actually* be “Poisson distributed” for them to have been generated by a Poisson process, but rather the residual outcomes once different predictors have been taken into account. The observed outcomes could be a combination of multiple Poisson processes; imagine that you have two groups, one of which generates Poisson-distributed events with rate of 20 and a the other with a rate of 5. Plotting hist(c(rpois(1e6, 5), rpois(1e6, 20))) clearly gives us a bimodal distribution in observed frequencies, but if you condition on group the ordinary Poisson distribution falls out. Still, I don’t think any combination of Poisson processes could give us that little upward blip at 9 without a tail stretching out into 10, 11, and 12…

But in any case, it’s probably not too terrible for our purposes here — except for the dairy outcome — since there are so few people eating 9 servings overall (and dairy’s what I’d expect it to be most problematic for, since milk, cheese, butter, etc. are in everything).

Also, eyeballing the amounts consumed, I don’t see any effect of treatment. It *is* clear that people eat a lot more dairy than they do fish, though.

On a side note, I’m not sure *servings* were the best unit of measurement to use when asking people about their eating habits. Speaking personally, I think of food I consume in units kcal, oz, grams, macronutrient amounts, etc., and given how most meat I encounter is labeled by weight (e.g. a burger with two *quarter-pound* patties, an *8 oz. * steak, chicken breasts that are *$4.99/lb*, and so on), I’d reckon most others think similarly. I’m not too familiar with the literature on this, but I also vaguely recall seeing plenty of articles describing how notoriously bad Americans are at thinking in servings, with illustrations involving steaks and packs of playing cards (e.g. most will say a 16 oz. steak is one serving, where it’s actually ~6ish, and an actual serving is the size of a standard deck). Additionally, weren’t “standard serving sizes” invented by the USDA for (outdated) Food-Pyramid-related reasons? They might then not have been appropriate for non-American respondents, and units of mass could have worked better (or, better still, allow survey takers to answer in whatever quantified units they were most comfortable with, and convert between units afterwards). You might also include a question to gauge certainty: something like “how confident are you in these numbers?” or “is your reported consumption here typical?” Then, since outcomes of mass or weight would be continuous and not discrete, you could accommodate measurement uncertainty according to these answers in Stan and use a Gaussian likelihood (though I dunno that you could account for zero-inflation with off-the-shelf stuff there). Finally, meats can obviously come in intermediate amounts — people don’t eat them one serving at a time, and can just as well eat 1.75 servings or whatever, so allowing responses to more faithfully reflect reality seems preferable (then again, we don’t want to fetishize precision). Oh well, it may not be too important, since the question clarifies that “one serving of meat or eggs is 3 ounces, about the size of a deck of cards”, but it might be something to consider when designing future surveys.

Another question could maybe also ask about how much money is spent on each meat on average, since “economics” are the mechanism by which future meat production gets incentived (so buying *filet mignon* could be a lot “worse” from the perspective of causing future suffering than buying the same amount of ground beef). Though this might be a little tricky in case exposure to the ad results in people buying more expensive “happy meat” out of increased concern for animal welfare.

Moving on, let’s look at consumption with respect to gender with the code:

consumption <- d[,c("Gender", "X2DaysPork", "X2DaysBeef", "X2DaysDairy", "X2DaysEggs", "X2DaysBird", "X2DaysSea")] %>% gather(Type, Amount, -Gender) consumption <- consumption[consumption[,1] != "",] ggplot(data = consumption, aes(x=Type, y=Amount, color=Gender)) + geom_jitter() + theme_bw()

And tell pretty much nothing, since the sample is mostly female.

What about a relationship with age? With the code:

consumption <- d[,c("c.Age", "X2DaysPork")] %>% gather(Type, Amount, -c.Age) porkAge <- ggplot(data = consumption, aes(x=c.Age, y=Amount)) + geom_jitter(alpha = .2) + theme_bw() + ggtitle("Pork") consumption <- d[,c("c.Age", "X2DaysBeef")] %>% gather(Type, Amount, -c.Age) beefAge <- ggplot(data = consumption, aes(x=c.Age, y=Amount)) + geom_jitter(alpha = .2) + theme_bw() + ggtitle("Beef") consumption <- d[,c("c.Age", "X2DaysDairy")] %>% gather(Type, Amount, -c.Age) dairyAge <- ggplot(data = consumption, aes(x=c.Age, y=Amount)) + geom_jitter(alpha = .2) + theme_bw() + ggtitle("Dairy") consumption <- d[,c("c.Age", "X2DaysEggs")] %>% gather(Type, Amount, -c.Age) eggsAge <- ggplot(data = consumption, aes(x=c.Age, y=Amount)) + geom_jitter(alpha = .2) + theme_bw() + ggtitle("Eggs") consumption <- d[,c("c.Age", "X2DaysBird")] %>% gather(Type, Amount, -c.Age) birdAge <- ggplot(data = consumption, aes(x=c.Age, y=Amount)) + geom_jitter(alpha = .2) + theme_bw() + ggtitle("Bird") consumption <- d[,c("c.Age", "X2DaysSea")] %>% gather(Type, Amount, -c.Age) seaAge <- ggplot(data = consumption, aes(x=c.Age, y=Amount)) + geom_jitter(alpha = .2) + theme_bw() + ggtitle("Sea") grid.arrange(porkAge, beefAge, dairyAge, eggsAge, birdAge, seaAge, ncol = 3)

We see no immediately discernible relationship with age (maybe that people in their early twenties are underrepresented). Oh well. There also looks to be an excess of 25-year-olds relative to 24, 23, etc. year-olds, too; hopefully options for age weren’t presented as a drop down menu with a maximum of 25!

Let’s look at the attitude questions, now, at least with respect to treatment:

attitude <- d[,c("Treat", "MeatUniqueSmart", "LessMeatGood", "MeatReplace")] %>% gather(Question, LikertAnswer, -Treat) attitude$LikertAnswer <- factor(attitude$LikertAnswer, levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")) ggplot(data = na.omit(attitude), aes(x=Question, y=LikertAnswer, color=Treat)) + geom_jitter() + theme_bw() #again, we omit NAs

Once more, it doesn’t look like there’s a very strong relationship with treatment, if any.

Do we care about people who already didn’t eat meat and still don’t eat meat? I guess we might, in case watching the ad persuaded them to consider trying it in the future or changed their opinions on the moral value of animals. Unfortunately, the only data we have on pre-existing vegetarianism is the question about whether people who don’t eat meat now didn’t eat it back then. It’s a shame there wasn’t a question that simply asked “Were you a vegetarian 4 months ago?”, since the data we do have necessarily underestimates the amount of pre-existing vegetarians (since some will have stopped eating meat, which is what we’re interested in). I guess we can assume that nobody stopped eating meat if they didn’t before, but that seems iffy (if it’s a false assumption, we partially open ourselves up to concomitant variable bias). Ah well, let’s briefly look at the people who didn’t eat meat before and don’t eat it now.

veg <- which(d$NoMeatAlready == "DidNotEatAlready") d$AlrVeg <- ifelse(d$NoMeatAlready == "DidNotEatAlready",1,0) #let's create another column of this too d$AlrVeg[which(d$anyChangeP4M == "" & d$AlrVeg == 0)] <- NA #preserve missing data attitudeVeg <- d[veg ,c("Treat", "MeatUniqueSmart", "LessMeatGood", "MeatReplace")] %>% gather(Question, LikertAnswer, -Treat) attitudeVeg$LikertAnswer <- factor(attitudeVeg$LikertAnswer, levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")) ggplot(data = na.omit(attitudeVeg), aes(x=Question, y=LikertAnswer)) + geom_jitter() + theme_bw()

Interesting, and about what we’d expect. Let’s also check:

summary(d$NowVs4MonPast[veg])

less same more NA’s

75 62 9 5

summary(d$NowVs4MonFut[veg])

none less same more NA’s

104 10 24 9 4

So despite not eating meat then and now, some people are eating more meat and many less meat. Weird! Also, it looks like some of our non-meat eaters will be eating no meat 4 months from now, but most will be eating none and some less. That’s good, at least!

I wonder how much meat the folks who don’t “eat meat, chicken, and fish” have actually eaten in the last few days?

par(mfrow = c(2,3)) hist(d$X2DaysPork[veg]); hist(d$X2DaysBeef[veg]); hist(d$X2DaysDairy[veg]); hist(d$X2DaysEggs[veg]); hist(d$X2DaysBird[veg]); hist(d$X2DaysSea[veg])

#or, alternatively consumptionVeg <- d[veg,c("Treat", "X2DaysPork", "X2DaysBeef", "X2DaysDairy", "X2DaysEggs", "X2DaysBird", "X2DaysSea")] %>% gather(Type, Amount, -Treat) ggplot(data = consumptionVeg, aes(x=Type, y=Amount)) + geom_jitter() + theme_bw()

So most but not all of them consume no servings, though many are consuming dairy or eggs. I’m not super sure how to deal with these, especially the ones who then report actually eating meat, but there are only 151 of them — enough for each to get a first gen, Kanto Pokemon.

The massaged data on pre-existing veg*nism are rather suspect, and I’m reluctant to include it in our models as a result, at least in fear of post-treatment bias — if the video affected rates of retention of veg*n dieting, then including it as a predictor could obscure the causal importance of treatment.

Anyway, let’s create a new data frame with just the variables we’re interested in, and then clean it up a bit more.

d2 <- d[, c("Treat", "X2DaysPork", "X2DaysBeef", "X2DaysDairy", "X2DaysEggs", "X2DaysBird", "X2DaysSea", "NowVs4MonPast", "NowVs4MonFut", "MeatUniqueSmart", "LessMeatGood", "MeatReplace", "Gender", "Country", "Age.s", "AlrVeg")] d2$Country[d2$Country == ""] <- NA d2$country_id <- coerce_index(d2$Country) d2$Gender[d2$Gender == ""] <- NA d2$gender_id <- coerce_index(d2$Gender) d2$Treat <- as.integer(d2$Treat) #change from factor to integer d2$Treat[d2$Treat == 1] <- 0 d2$Treat[d2$Treat == 2] <- 1 d2$Treat <- as.integer(d2$Treat) #change from factor to integer d2$AlrVeg <- as.integer(d2$AlrVeg)

We’re unfortunately doing a complete-cases analysis, as mentioned earlier, since Stan has trouble with discrete variable imputation. Let’s hope the data’s MCAR! I plan to revisit this issue later and incorporate multivariate outcomes, measurement uncertainty, missing data imputation, etc. I know how to do the latter two in Stan for continuous variables, but not for discrete, and I know in principle how to do it for discrete variables, so it shouldn’t be too hard (my code initially assumed you could, but R/Stan kept throwing errors, so I consulted my friendly neighborhood statistician and he informed me of Stan’s limitations).

dcc <- d2[ complete.cases(d2) , ]

The most interesting outcome variable to examine is, I think, the amount of meat consumed in the past 2 days, which I’ll say is zero-inflated Poisson distributed. This seems reasonable, as our outcomes are non-negative integers, and the event rate (serving consumption) and mean are low. The mean is supposed to be a bit less than the variance, too, but that’s in reference to the distribution of residuals after conditioning on the predictors, and not the raw outcomes; we can probably check into that later. Why a zero-inflated Poisson (ZIP)? Some people can have eaten 0 servings both by chance, and because they don’t eat particular sorts of meat due to veg*n (or religious, cultural, etc.) dietary preferences. A ZIP mixed model has a likelihood that’s involves a binomial process determining whether an individual is a candidate for accumulating events according to a Poisson process. A picture might serve to better illustrate:

Damn good special effects, eh? I’m not quite sure yet how to do regression with multivariate outcomes (to accommodate correlations between parameter estimates for each outcome; if outcomes are tightly correlated but noisy, I think you run an overfitting risk that traditional model comparison like IC are not robust too, since you can just keep fitting models with different outcomes until you get a result you want). Furthermore, I’m not even sure how I’d write a multivariate model — I think there’d be positive association between different values for lambda (on average), since individuals who eat a ton of meat of one sort would be likely to eat a ton of meat of another sort. On the other hand, realized outcomes on shorter timescales would probably be inversely correlated, since, if an individual chooses a particular meat to consume, that means they won’t be eating a different type of meat. Actual consumption of each type of meat is not independent, so we can’t just model multivariate outcomes as multiple independent Poisson processes with positively correlated lambdas (e.g. consider the extreme case, where somebody cooks a turkey for Thanksgiving and then eats turkey sandwiches every day for the next week. People will eat until a certain point, and if they’ve filled up on one thing won’t separately decide whether to eat another).

Similarly, the “probability of abstention” parameter *p* for each sort of meat would have a component shared between meats (e.g. if individuals have a veg*n diet) and a component independent for each meat (e.g. pig’s not ḥalāl or kosher, but beef can be, I think. Or allergies, or something). I reckon you’d want to infer two multivariate distributions of *p*s and *lambda*s; for the covariance matrix of the *p*s, off-diagonal elements would represent veg*n dieting, whereas for the lambdas they’d represent the “trade-off” of choosing to eat one sort of meat over another on shorter time-scales. But it would average out on longer time scales, so you’d need to keep redrawing day-by-day lambdas. Diagonal elements would then represent overall rates of abstention and consumption for each sort of meat, and maybe you could pull a scalar multiplier out for the lambdas as a general term for how much a person eats. Alternatively, maybe you could do a ZIP for total meat consumption, and then draw from a multinomial or multivariate hypergeometric distribution (depending on whether we think people get tired of a meat or not?) to determine which meat-types people ate. Maybe consumption amounts and probabilities of abstention for each type of meat could then be inferred?

IDK that *rethinking* and Stan currently have any of this functionality, though, so let’s just look at outcomes on a univariate, meat-by-meat basis.

But for now which outcomes should we then focus on? From the google doc, it seems the vegetarian advertisement used was the one at “meatvideo.com” titled “What Cody Saw Will Change Your Life“. Watching this video, pigs are the first animals to be discussed, followed by chickens. It seems reasonable that viewers will have been more likely to watch the first part of the video than any other distinct part if they clicked away without finishing it. Furthermore, pigs seem pretty smart — smarter than chickens, cattle, and fish, at least to my uneducated eye — and so perhaps are more deserving of our (or rather *my*) moral concern. Chickens, meanwhile, are tiny, at least compared to pigs and cows (as I discussed briefly here), and so perhaps each reduction in serving of chicken is proportionally more important.

Anyway, let’s just look at pigs and chickens/birds (it’ll also cut computation time down a fair bit). A few attempts to examine the data have summed across all food categories to get composite consumption scores, but this seems silly in that it 1) throws away a lot of resolution in the data, and 2) assumes equal weighting of servings, which strikes me as wrong for reasons mentioned in the previous paragraph. Different mechanisms might also apply to different animals, which lumping them together might obscure (e.g. shifts in diet whereby you replace one type of meat with another, or a differential influence of age, country, etc. on what foods you eat).

There are a lot of possible predictor variables (age, country, gender, treatment, maybe-biased-pre-existing vegetarianism), and so a TON of possible models we can consider, especially since variables can interact in a bunch of crazy ways; e.g. maybe treatment effect is only substantial in older men in the UK. That’s getting a bit silly, but it certainly seems plausible that, say, younger people might be more easily persuaded than older people (having had less time to internalize “carnism” or whatever), or that, I dunno, non-gender-binary people will be more able to sympathize with the plights of oppressed animals and so more likely to change their dietary habits. We don’t want to just fit every possible model, though because we might eventually find one that fits really well purely by chance in a way that fools information criteria and other model comparison tools. So instead, I’ll fit a bunch of models that seem intuitively plausible to me, and then average them together after the fact. Hopefully this approach can capture the complexity of the relationships between different variables, while also avoiding peculiar sorts of overfitting.

That said, I fit the below multilevel/hierarchical mixed models with Hamiltonian Monte Carlo (since something like quadratic approximation wouldn’t work — in exchange, having varying instead of fixed intercepts should help guard against overfitting!). As there were a fair few models to consider, I didn’t go super in-depth with respect to checking MCMC health and instead just visually inspected trace plots to check that they fit the “fuzzy caterpillar” criterion, ensured that the Gelman-Rubin Convergence Diagnostic was one for almost all my parameters (once or twice for relatively unimportant parameters it was 1.01, which is suspicious but not catastrophic), and looked at how many effective samples I drew (rarely in the low hundreds despite running my chain for 10-20k iterations, but almost always I got at least several thousand samples) (though do note: THIS IS BAD! IF THIS WERE SOMETHING I WERE ACTUALLY INTERESTED IN PUBLISHING, I TOTALLY SHOULD HAVE DONE MORE! But since this is just a fun side project, I was content with limited assurance and the slight possibility of pathological chains). Sometimes, things looked a wee grim, so I reran the chains with substantially more samples and a greater warm-up. I also used only a single chain, which is ALSO BAD from the perspective of diagnosing MCMC pathology since multiple chains are a good way to check that you’re converging on the true joint posterior (I will remark, however, that my asshole computer crashed while I was inspecting stuff before I could save my workspace, and inferred parameters and WAIC weights were very similar to the ones I inferred when rerunning the code again, and that similar structures in different models had very similar distributions). I and many others eagerly await Mike’s release of **BONSAI** (Bayesian Output Needs Semi-Automated Inspection), as it oughtta really streamline the process of diagnosing chain performance and health.

I used adaptively regularizing/multilevel priors, but when I had to put money down I chose conservatively. They were, in some sense, empirical, in that they represented values I found plausible (e.g. my prior for the intercept on *p* was “normal” before the logit link, and after it had an actual median of about ~12% and a +/- 47.5% interval ranging from ~5% to 27%. My prior for the intercept on lambda after the log link (for individuals who did not abstain from particular sorts of meat) had a median of *e *with 95% intervals from ~.4 to 20). Coefficients on predictors were always centered on 0 with decently small variance, since I don’t expect gender, country, treatment, age (in sds), or any interactions to have too great an effect, and having tight regularizing priors avoids overfitting (do note — I used adaptively regularizing priors/varying intercepts for discrete variables with more than 2 states, i.e. gender and country. In other words, I set hyperpriors on those intercept and interaction parameters and learned the prior from the data itself, or alternatively inferred the population-level distribution from which individual country and gender main and interaction effects were drawn. The prior on sd was a fairly tight half-Cauchy, and I think it worked OK overall, but since the Cauchy has infinite variance and thick tails it might be a bit overfit with the limited number of categories available. An exponential prior may have been preferable, and should be explored in sensitivity analysis). You can see my specification of the actual distributions below, along with some quick visualizations.

*Usually* the prior’s overwhelmed by the likelihood, in any case, if you have reasonably informative data, and I’ll do a full sensitivity analysis later to check that I get roughly the same results. I did, however, do a very *basic* sort of sensitivity analysis wherein I ran my best-supported pork model with very flat (“uninformative”, though not really — weighting values equally doesn’t make them “unweighted”, imo) priors, and compared parameter estimates from there to those obtained from my “empirical”, more informative priors. They looked ok. You can see for yourself (this is a bit jumping the gun, because we’ve not “actually” run any models yet, but bear with me):

**# Flat, Awful Priors (m10cc.p.sens)**

*pI ~ dnorm(0, 2),*

* lI ~ dnorm(0,5),*

* c(lT, lA, lAT) ~ dnorm(0,5),*

* lG[gender_id] ~ dnorm(0,sigGl),*

* sigGl ~ dcauchy(0,1.5),*

* lGT[gender_id] ~ dnorm(0,sigGTl),*

* sigGTl ~ dcauchy(0,1.5),*

* lC[country_id] ~ dnorm(0,sigCl),*

* sigCl ~ dcauchy(0,1.5)*

**#Informative, Decent Priors (m10cc.p)**

* pI ~ dnorm(-2, 0.5),*

* lI ~ dnorm(1,1),*

* c(lT, lA, lAT) ~ dnorm(0,1),*

* lG[gender_id] ~ dnorm(0,sigGl),*

* sigGl ~ dcauchy(0,0.5),*

* lGT[gender_id] ~ dnorm(0,sigGTl),*

* sigGTl ~ dcauchy(0,0.5),*

* lC[country_id] ~ dnorm(0,sigCl),*

* sigCl ~ dcauchy(0,0.5)*

plot(coeftab(m10cc.p, m10cc.p.sens))

So not *too* shabby (note: the open circles are MAP estimates and the lines represent 89% percentile intervals). The estimates are all highly similar, at least qualitatively, suggesting that there’s more than enough information in the data to overwhelm a crappy prior (presumably, we don’t actually assign a substantial amount of probability to people eating more than the US’ entire pig population every day). Some estimates are more confident that others; this is probably in part a reflection to how much information in the data is capable of informing those estimates (e.g. genders don’t vary too much in the sample, since it’s mostly female, so it’s hard to get a sense how much gender affects consumption). Later on, I’d like to rerun the entire analysis with different priors and examine how much the actual results change.

It was requested that I provide graphs of the priors used in the pork and bird models, in case their visualization is tricky or not immediate. You can play around with the following code to see how the priors shift visually as we define them in different ways. Priors on effects are graphed in units of relative change (with a line through 1, signifying no change), because absolute change depends on the total sum of the components of the linear model (since we have those pesky logit and log link functions).

#Plotting priors #pI #Prior for the intercept on p dens(logistic(rnorm(100000, -2, 0.5)) * 100, xlab = ("Percent Abstention")) title("Prior for the intercept on p")

#pT, pA, pAT, etc. #Prior for relative effect of Treatment, Age, etc. dens(exp(rnorm(100000, 0, 0.5)), xlab = "Proportional Odds") title("Prior for effect sizes of predictor coefficients on p") abline(v = 1, lty = 2)

#Prior for pC, pG, pGT, etc. sigs <- rcauchy(n = 3000, location = 0, scale = 0.5) #stan looks at half-cauchy defined on + reals by setting <lower = 0> sigs <- sigs[sigs > 0] dens(exp(rnorm(1000, 0, sigs[1])), xlab = "Proportional Odds", xlim = c(0,10), ylim = c(0,4), col=col.alpha(1,0.05)) for (i in 2:length(sigs)){ dens(exp(rnorm(1000, 0, sigs[i])), add = T, col=col.alpha(1,0.05)) } abline(v = 1, lty = 2) title(main = "Gender/Country on p Prior (defined by hyperprior)")

#lI dens(exp(rnorm(100000, 1, 1)), xlab = "Number of Servings Rate (per 2 days)", xlim = c(0,20)) title("Prior for intercept lI on lambda")

#lT, lA, lAT dens(exp(rnorm(100000, 0, 1)), xlab = "Proportional Change in Servings Rate (per 2 days)", xlim = c(0,20)) abline(v = 1, lty = 2) title("Prior for effects (lT, lA, lAT) on lambda (proportional change)")

And the priors on lG, lC, lGT, etc. are the same as for pC and pG, only now it’s proportional change rather than proportional odds. Ultimately, though, the priors don’t really matter *too* much, and if you disagree with my choice in them, perhaps you’ll be assured by the similarity in estimates obtained when I made the priors much, much flatter above.

So hopefully that’s all clear! For the actual analysis, I first fit the following models (including models without predictors for treatment, to check for overfitting). Sorry for the weird numbering!:

m0cc.p <- map2stan( #just intercepts, most simple model alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI, log(lambda) <- lI, c(pI) ~ dnorm(-2,0.5), c(lI) ~ dnorm(1,1) ) , data= list(X2DaysPork = dcc$X2DaysPork), iter = 25000, warmup = 5000 ) precis(m0cc.p)

m1cc.p <- map2stan( #intercepts and treatment alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat, log(lambda) <- lI + lT*Treat, pI ~ dnorm(-2, 0.5), pT ~ dnorm(0,0.5), lI ~ dnorm(1,1), lT ~ dnorm(0,1) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork), iter = 25000, warmup = 5000 ) precis(m1cc.p)

m2cc.p <- map2stan( #no interactions, just main effects for treatment, gender, age, and country alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat + pG[gender_id] + pA*Age_s + pC[country_id], log(lambda) <- lI + lT*Treat + lG[gender_id] + lA*Age_s + lC[country_id], pI ~ dnorm(-2, 0.5), c(pT, pA) ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), c(lT, lA) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 40000, warmup = 8000 ) precis(m2cc.p, depth = 2)

m3cc.p <- map2stan( #no interactions just Treatment and Gender effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat + pG[gender_id], log(lambda) <- lI + lT*Treat + lG[gender_id], pI ~ dnorm(-2, 0.5), pT ~ dnorm(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), lT ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id), iter = 25000, warmup = 5000 ) precis(m3cc.p, depth = 2)

m4cc.p <- map2stan( #no interactions just Treatment and Country effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat + pC[country_id], log(lambda) <- lI + lT*Treat + lC[country_id], pI ~ dnorm(-2, 0.5), pT ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), lT ~ dnorm(0,1), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m4cc.p, depth = 2)

m5cc.p <- map2stan( #no interactions, just Treatment and Age effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat + pA*Age_s, log(lambda) <- lI + lT*Treat + lA*Age_s, pI ~ dnorm(-2,0.5), c(pT, pA) ~ dnorm(0,0.5), lI ~ dnorm(1,1), c(lT, lA) ~ dnorm(0,1) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, Age_s = dcc$Age.s), iter = 25000, warmup = 5000 ) precis(m5cc.p, depth = 2)

m6cc.p <- map2stan( #A-T interaction and main effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat + pAT*Treat*Age_s + pG[gender_id] + pA*Age_s + pC[country_id], log(lambda) <- lI + lT*Treat + lAT*Age_s*Treat + lG[gender_id] + lA*Age_s + lC[country_id], pI ~ dnorm(-2, 0.5), c(pT, pA, pAT) ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m6cc.p, depth = 2)

m7cc.p <- map2stan( #A-T and A-G interaction and main effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pATG*Treat + pG[gender_id] + pA*Age_s + pC[country_id], pATG <- pT + pAT * Age_s + pGT[gender_id], log(lambda) <- lI + lATG*Treat + lG[gender_id] + lA*Age_s + lC[country_id], lATG <- lT + lAT * Age_s + lGT[gender_id], pI ~ dnorm(-2, 0.5), c(pT, pA, pAT) ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), pGT[gender_id] ~ dnorm(0,sigGTp), sigGTp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lGT[gender_id] ~ dnorm(0,sigGTl), sigGTl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m7cc.p, depth = 2)

m8cc.p <- map2stan( #all interactions with T and main effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pATGC*Treat + pG[gender_id] + pA*Age_s + pC[country_id], pATGC <- pT + pAT * Age_s + pGT[gender_id] + pCT[country_id], log(lambda) <- lI + lATGC*Treat + lG[gender_id] + lA*Age_s + lC[country_id], lATGC <- lT + lAT * Age_s + lGT[gender_id] + lCT[country_id], pI ~ dnorm(-2, 0.5), c(pT, pA, pAT) ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), pGT[gender_id] ~ dnorm(0,sigGTp), sigGTp ~ dcauchy(0,0.5), pCT[country_id] ~ dnorm(0,sigCTp), sigCTp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lGT[gender_id] ~ dnorm(0,sigGTl), sigGTl ~ dcauchy(0,0.5), lCT[country_id] ~ dnorm(0,sigCTl), sigCTl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m8cc.p, depth = 2)

m9cc.p <- map2stan( #A-T and A-G interaction on lambda and main effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pT*Treat, log(lambda) <- lI + lATG*Treat + lG[gender_id] + lA*Age_s + lC[country_id], lATG <- lT + lAT * Age_s + lGT[gender_id], pI ~ dnorm(-2, 0.5), pT ~ dnorm(0,0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lGT[gender_id] ~ dnorm(0,sigGTl), sigGTl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m9cc.p, depth = 2)

m10cc.p <- map2stan( #A-T and A-G interaction on lambda and main effects; no effect on p of treatment alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI, log(lambda) <- lI + lATG*Treat + lG[gender_id] + lA*Age_s + lC[country_id], lATG <- lT + lAT * Age_s + lGT[gender_id], pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lGT[gender_id] ~ dnorm(0,sigGTl), sigGTl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m10cc.p, depth = 2)

m12cc.p <- map2stan( #A-T, G-T, C-T interactiosn on lambda and main effects; no effect on p of treatment alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI, log(lambda) <- lI + lATGC*Treat + lG[gender_id] + lA*Age_s + lC[country_id], lATGC <- lT + lAT * Age_s + lGT[gender_id] + lCT[country_id], pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lGT[gender_id] ~ dnorm(0,sigGTl), sigGTl ~ dcauchy(0,0.5), lCT[country_id] ~ dnorm(0,sigCTl), sigCTl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m12cc.p, depth = 2)

m13cc.p <- map2stan( #A-T interaction on lambda and main effects; no effect on p of treatment alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI, log(lambda) <- lI + lT*Treat + lAT * Age_s * Treat + lG[gender_id] + lA*Age_s + lC[country_id], pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), c(lT, lA, lAT) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m13cc.p, depth = 2)

m14cc.p <- map2stan( #A-G interaction on lambda and main effects; no effect on p of treatment alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI, log(lambda) <- lI + lTG*Treat + lG[gender_id] + lA*Age_s + lC[country_id], lTG <- lT + lGT[gender_id], pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), c(lT, lA) ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lGT[gender_id] ~ dnorm(0,sigGTl), sigGTl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(Treat = dcc$Treat, X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m14cc.p, depth = 2)

# ----- #NO TREATMENT EFFECTS ----- # m2cc.p.nt <- map2stan( #no interactions, just main effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pG[gender_id] + pA*Age_s + pC[country_id], log(lambda) <- lI + lG[gender_id] + lA*Age_s + lC[country_id], pI ~ dnorm(-2, 0.5), pA ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), lA ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m2cc.p.nt, depth = 2)

m3cc.p.nt <- map2stan( #no interactions, just Gender effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pG[gender_id], log(lambda) <- lI + lG[gender_id], pI ~ dnorm(-2, 0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5) ) , data= list(X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id), iter = 25000, warmup = 5000 ) precis(m3cc.p.nt, depth = 2)

m4cc.p.nt <- map2stan( #no interactions, just Country effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pC[country_id], log(lambda) <- lI + lC[country_id], pI ~ dnorm(-2, 0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(X2DaysPork = dcc$X2DaysPork, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m4cc.p.nt, depth = 2)

m5cc.p.nt <- map2stan( #no interactions, just Age effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pA*Age_s, log(lambda) <- lI + lA*Age_s, pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), pA ~ dnorm(0,0.5), lA ~ dnorm(0,1) ) , data= list(X2DaysPork = dcc$X2DaysPork, Age_s = dcc$Age.s), iter = 25000, warmup = 5000 ) precis(m5cc.p.nt, depth = 2)

m7cc.p.nt <- map2stan( #gender, age, and country main effects alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI + pG[gender_id] + pA*Age_s + pC[country_id], log(lambda) <- lI + lG[gender_id] + lA*Age_s + lC[country_id], pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), pA ~ dnorm(0,0.5), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,0.5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,0.5), lI ~ dnorm(1,1), lA ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 25000, warmup = 5000 ) precis(m7cc.p.nt, depth = 2)

m11cc.p.nt <- map2stan( #no treatment effects on p or l, no interactions, just main effects on lambda alist( X2DaysPork ~ dzipois( p , lambda ), logit(p) <- pI, log(lambda) <- lI + lG[gender_id] + lA*Age_s + lC[country_id], pI ~ dnorm(-2, 0.5), lI ~ dnorm(1,1), lA ~ dnorm(0,1), lG[gender_id] ~ dnorm(0,sigGl), sigGl ~ dcauchy(0,0.5), lC[country_id] ~ dnorm(0,sigCl), sigCl ~ dcauchy(0,0.5) ) , data= list(X2DaysPork = dcc$X2DaysPork, gender_id = dcc$gender_id, Age_s = dcc$Age.s, country_id = dcc$country_id), iter = 30000, warmup = 6000 ) precis(m11cc.p.nt, depth = 2)

save.image("C:/Users/Nikolai/Dropbox/R Files/Pork.RData")

Phew, that took a while. Some of the models/chains didn’t converge/mix well, so I upped the warmup and total number of iterations and reran them, after which they did. I looked through the trace plots of my different parameters and the chains all looked pretty healthy (occasionally with a spike here or there); I won’t reproduce them all here, since that would involve sharing hundreds of graphs, but here are the traces for the two best supported models (according to WAIC weight; see below):

That little blip around iteration ~17k is a teensy bit worrisome, but ultimately harmless, imo (there were similar blips in some of the trace plots for the other models). We can also look at the marginal posterior distributions for and correlations between different parameters; for m10cc.p, we get:

Not ideal (which would be very little correlation between parameters when making draws from the joint posterior), but not terrible, specifically because I don’t actually care what the estimates for each parameter are. You might think, wait, aren’t we interested in what the coefficient for the main effect of treatment on lambda is? Doesn’t the multicollinearity in parameters estimates preclude us from drawing firm inferences from the marginal distributions of each parameter? To which I’d say, sure, but who does that, anyway? Multicollinearity and non-identifiability are issues (especially when they slow down your MCMC, which they probably did here, tbh), but in this case we’re not going to try to look at marginal posteriors and try to divine what’s going on — we’re just going to simulate from it, which, if (and since) our chains were healthy, should work just fine.

But first let’s compare the models with WAIC (feel free to use whatever model comparison device you like — again, in a actual analysis, I would use a bunch of them — BIC, DIC, Bayes Factors, rjMCMC, etc. — and see where they agree and disagree, but I won’t here since that’d be too much work. WAIC seems to me one of the better information criteria, since it allows for informative priors and a non-multivariate normal joint posterior, and also makes use of the entire posterior, and not just its tip. Personally, I think rjMCMC is a more elegant approach to model-averaging, but I’m not currently sure how to do that within a linear regression framework).

compare(m1cc.p, m2cc.p, m3cc.p, m4cc.p, m5cc.p, m6cc.p, m7cc.p, m8cc.p, m9cc.p, m10cc.p, m11cc.p.nt, m12cc.p, m13cc.p, m14cc.p, m2cc.p.nt, m3cc.p.nt, m4cc.p.nt, m5cc.p.nt, m7cc.p.nt)

So it looks like we’re pretty uncertain about which model is most consistent with the data, but models without treatment predictor variables have all the WAIC weight, so they’re overall better supported. What sorts of parameters have been estimated, and with what MAP (*maximum a* *posteriori — *i.e. the marginal posterior mode) estimates?

coeftab(m13cc.p, m12cc.p, m9cc.p, m10cc.p, m7cc.p, m6cc.p, m8cc.p)

Pretty tricky, right! Look, for example, at lT — the main effect of treatment on lambda, i.e. the rate parameter of our Poisson process. It’s positive in all models, suggesting that exposure to the treatment increases the rate at which people consume pigs. But it varies somewhat from model to model, because it represents different things in each — in one, it might represent the average effect of treatment across all ages (if no age predictor were included in the equation for lambda); in others, it might be the effect when Age_s = 0, that is to say, when an individual is of average age (*roughly*, at least. I standardized before removing incomplete cases, so *0* might not actually represent the average age in the data that I’m actually using. It doesn’t really matter for the results, since I can arbitrarily change the scale of my predictors however I want, but it does make a interpretation on the altered scale a little iffier and the model will adjust accordingly so long as I ensure that my priors also roughly match). Additionally, there are tons of interaction terms that make drawing meaningful conclusions from reading off tables of estimates impossible (and also consider the issue of multicollinearity mentioned before). Plus, since these are MAPs, we have no sense of the uncertainty about each estimate. Let’s instead simulate *p*s and *lambda*s from these models and use them to simulate individuals who eat pork who were or were not exposed to the treatment.

But first, who do we want to simulate? Let’s focus on American women of average age (in the dataset), since they’re the largest contingent to target here (with a slight bit more work, we could generate observations from a multivariate distribution of genders, ages, and countries that reflect that seen in the data, but let’s just go with this for now).

d.pred <- data.frame( Treat = c(0,1), gender_id = 1, country_id = 6, Age_s = 0, )

Normally, we can get our *lambda*s and *p*s using the ensemble() function, but unfortunately because of the way I constructed the models (by pulling interaction terms out of the main equation), the outputs don’t have the same number of dimensions, so I have to go in and do it manually. I’ll use the top 5 models here, collectively accounting for 96% of the WAIC weight:

#simulate values for p and lambda a <- list() a[[1]] <- as.data.frame(link(data = d.pred, m10cc.p, n = 1e4)) #.41 a[[2]] <- as.data.frame(link(data = d.pred, m13cc.p, n = 1e4)) #.31 a[[3]] <- as.data.frame(link(data = d.pred, m9cc.p, n = 1e4)) # .12 a[[4]] <- as.data.frame(link(data = d.pred, m12cc.p, n = 1e4)) # .06 a[[5]] <- as.data.frame(link(data = d.pred, m7cc.p, n = 1e4)) #.06 #sample proportional to WAIC weights b <- list() b[[1]] <- sample_n(as.data.frame(a[[1]]), size = 4100) b[[2]] <- sample_n(as.data.frame(a[[2]]), size = 3100) b[[3]] <- sample_n(as.data.frame(a[[3]]), size = 1200) b[[4]] <- sample_n(as.data.frame(a[[4]]), size = 600) b[[5]] <- sample_n(as.data.frame(a[[5]]), size = 600) #put it all in one data frame values <- b[[1]][,c("p.1", "p.2", "lambda.1", "lambda.2")] for (i in 2:5) { values <- rbind(values, b[[i]][,c("p.1", "p.2", "lambda.1", "lambda.2")])} #plot the distributions par(mfrow = c(2,1)) #ps dens(values$p.1) dens(values$p.2, add = T, col = 2) title("Model-Averaged Estimates for p for Pork Consumption") legend(.403,25, c("Control","Treatment"), lty=c(1,1), lwd=c(2.5,2.5),col=c(1,2)) #lambdas dens(values$lambda.1, xlim = c(1.4, 2.3)) dens(values$lambda.2, add = T, col = 2) title("Model-Averaged Estimates for Lambda for Pork Consumption") legend(2.2,5, c("Control","Treatment"), lty=c(1,1), lwd=c(2.5,2.5),col=c(1,2))

So it looks like exposure to the treatment didn’t really affect rates of abstention from pork (by persuading people to go veg*n or converting them religiously or something), but did motivate individuals to consume higher amounts of pork. How much higher? We can’t really guesstimate it from looking at the bottom plot, since the visual overlap doesn’t actually represent the difference between control and treatment. Instead, we must compute the distribution of the difference explicitly.

#difference in lambda between control and treatment par(mfrow = c(1,1)) dens(values$lambda.2 - values$lambda.1) title("Estimated Difference in Lambda Between Control and Treatment")

So it looks like the models are decently confident that watching the vegetarian advertisements causes average-aged American women to increase their consumption of pork to the tune of about a fifth of a serving per 2 days. What’s the probability that it increases pork consumption?

diff <- values$lambda.2 - values$lambda.1 sum( diff > 0 ) / length( diff )

About 91% (that’s how much probability is piled to the right of 0 in the figure above).

Lots of the models WAIC preferred had terms for an interaction between age and treatment, so let’s see how the effects of treatment on *p* and *lambda* varied across different ages. To do this, we can repurpose the code from before but have it make predictions across a variety of ages, and then plot those predictions alongside intervals to represent our uncertainty. For ease of understanding, I also convert standardized ages back into their original scale (in years) when I do the plotting:

## Estimate lambda and p for a variety of ages ## age.seq <- seq(from = -1.5, to = 2, by = .1) d.pred.NT <- data.frame( Treat = 0, gender_id = 1, country_id = 6, Age_s = age.seq ) d.pred.T <- data.frame( Treat = 1, gender_id = 1, country_id = 6, Age_s = age.seq ) #simulate values for p and lambda for NT aNT <- list() aNT[[1]] <- as.data.frame(link(data = d.pred.NT, m10cc.p, n = 1e4)) #.41 aNT[[2]] <- as.data.frame(link(data = d.pred.NT, m13cc.p, n = 1e4)) #.31 aNT[[3]] <- as.data.frame(link(data = d.pred.NT, m9cc.p, n = 1e4)) # .12 aNT[[4]] <- as.data.frame(link(data = d.pred.NT, m12cc.p, n = 1e4)) # .06 aNT[[5]] <- as.data.frame(link(data = d.pred.NT, m7cc.p, n = 1e4)) #.06 #sample proportional to WAIC weights for NT bNT <- list() bNT[[1]] <- sample_n(as.data.frame(aNT[[1]]), size = 4100) bNT[[2]] <- sample_n(as.data.frame(aNT[[2]]), size = 3100) bNT[[3]] <- sample_n(as.data.frame(aNT[[3]]), size = 1200) bNT[[4]] <- sample_n(as.data.frame(aNT[[4]]), size = 600) bNT[[5]] <- sample_n(as.data.frame(aNT[[5]]), size = 600) #simulate values for p and lambda for T aT <- list() aT[[1]] <- as.data.frame(link(data = d.pred.T, m10cc.p, n = 1e4)) #.41 aT[[2]] <- as.data.frame(link(data = d.pred.T, m13cc.p, n = 1e4)) #.31 aT[[3]] <- as.data.frame(link(data = d.pred.T, m9cc.p, n = 1e4)) # .12 aT[[4]] <- as.data.frame(link(data = d.pred.T, m12cc.p, n = 1e4)) # .06 aT[[5]] <- as.data.frame(link(data = d.pred.T, m7cc.p, n = 1e4)) #.06 #sample proportional to WAIC weights for T bT <- list() bT[[1]] <- sample_n(as.data.frame(aT[[1]]), size = 4100) bT[[2]] <- sample_n(as.data.frame(aT[[2]]), size = 3100) bT[[3]] <- sample_n(as.data.frame(aT[[3]]), size = 1200) bT[[4]] <- sample_n(as.data.frame(aT[[4]]), size = 600) bT[[5]] <- sample_n(as.data.frame(aT[[5]]), size = 600) #put it all in one data frame for NT valuesNT <- bNT[[1]][,c(paste0("p.", seq(1:36)), paste0("lambda.", seq(1:36)))] for (i in 2:5) { valuesNT <- rbind(valuesNT, bNT[[i]][,c(paste0("p.", seq(1:36)), paste0("lambda.", seq(1:36)))])} #put it all in one data frame for T valuesT <- bT[[1]][,c(paste0("p.", seq(1:36)), paste0("lambda.", seq(1:36)))] for (i in 2:5) { valuesT <- rbind(valuesT, bT[[i]][,c(paste0("p.", seq(1:36)), paste0("lambda.", seq(1:36)))])} valuesNT.mean <- apply(valuesNT, 2, mean) valuesT.mean <- apply(valuesT, 2, mean) valuesNT.HPDI <- apply(valuesNT, 2, HPDI) valuesT.HPDI <- apply(valuesT, 2, HPDI) #convert back to original age values trueAgeSeq <- age.seq * sd(d$c.Age) + mean(d$c.Age) #plotting change in p with age plot(trueAgeSeq, valuesNT.mean[1:36], type = "l", ylim = c(0.3, 0.37), lwd = 2, xlab = "Age", ylab = "Predicted p") lines(trueAgeSeq, valuesT.mean[1:36], col = 2) shade(valuesNT.HPDI[,1:36], trueAgeSeq) shade(valuesT.HPDI[,1:36], trueAgeSeq, col = col.alpha(2, alpha = .1)) legend(22.8,.373, c("Control","Treatment"), lty=c(1,1), lwd=c(2.5,2.5),col=c(1,2), bty = "n") title("Estimated Effects of Exposure to Ad on Probability of Pork Abstention") #plotting change in lambda with age plot(trueAgeSeq, valuesNT.mean[37:72], type = "l", ylim = c(1, 3), lwd = 2, xlab = "Age", ylab = "Predicted λ") lines(trueAgeSeq, valuesT.mean[37:72], col = 2, lwd = 2) shade(valuesNT.HPDI[,37:72], trueAgeSeq) shade(valuesT.HPDI[,37:72], trueAgeSeq, col = col.alpha(2, alpha = .1)) legend(22.5,3, c("Control","Treatment"), lty=c(1,1), lwd=c(2.5,2.5),col=c(1,2)) title("Estimated Effects of Exposure to Ad on Rate of Pork Consumption")

In both these graphs, solid lines represent mean posterior estimates for *p* and *lambda*, and shaded regions represent 89% HPDIs (*highest posterior density intervals*). We can see that treatment has very little estimated effect on *p*, with the control group staying about constant with respect to inferred probabilities of abstention and the treatment group having slightly higher rates of abstention in older individuals and slightly lower rates of abstention in younger individuals. Whatever effect that’s present, however, is very small and rather uncertain.

A much clearer relationship between age and treatment effect is seen in our estimates for *lambda*, the Poisson rate parameter. Here, we observe a distinct inverse relationship in the control group between pork consumption and age, with non-abstaining younger individuals eating more pork than older individuals. However, when exposed to the treatment, younger individuals consume less pork and older individuals consume more pork. This is pretty bizarre, but lots of explanatory hypotheses immediate jump to mind — perhaps younger individuals are easier to persuade, and older individuals more likely to take offense (consider “for every pig you don’t eat, I’ll eat three!”)? We can explore these later on, when we look at the attitude data.

Let’s also look at the actual distribution of treatment effects on lambda as individuals age (code adapted from before):

Of course, since we sampled from the posterior with MCMC, the true posterior distributions for lambda probably aren’t quite so bumpy. But you get the gist (also note: negative values here represent a reduction in consumption rate and positive values an increase in consumption rate, with approximate age in the top right corner. The number after the “% Neg:” is the proportion of the distribution that’s negative — i.e. to the left of the dashed line at 0).

Throughout all this, we’ve been acting at the level of parameters, and not individuals and **servings**. Let’s simulate individuals and examine their behavior after being exposed to the ad (averaging across the entire posterior and our uncertainty in model choice). Thankfully, we can do this with ensemble():

sims <- ensemble(m13cc.p, m12cc.p, m9cc.p, m10cc.p, m7cc.p, m6cc.p, m8cc.p, data = d.pred, do_link = F, do_sim = T, n = 5e3)$sim

Now we have 5,000 simulated individuals who were exposed to the control, and then alternatively to the treatment. How many servings of pork did each alternative eat? On average, those exposed to the treatment consumed:

mean(sims[,2]) - mean(sims[,1])

0.057 more servings of pork in the two day period, once you account for the potentially very minor (perhaps correlated in some weird way) difference in rates of abstention. This may seem small, but if it represents consistent behavior the numbers could very quickly add up. It’s also consistent with our earlier examination of lambda.

Or, if you prefer a visual representation:

hist(sims[,1], col=rgb(1,0,0,0.3), main = "Simulated Pork Consumption", xlab = "Number of Servings") hist(sims[,2], col=rgb(0,0,1,0.3), add=T) legend(7,1800, c("Control","Treatment"), lty=c(1,1), lwd=c(2.5,2.5),col=c(rgb(1,0,0,0.3),rgb(0,0,1,0.3)))

So, the pork story is a bit disappointing; it looks like clicking on this particular vegetarian ad causes average-aged American women to increase their pork consumption. I wonder why? Later on, when I explore the effect of treatment on attitude, perhaps we’ll see (if I had to guess and armchair psychoanalyze, I’d say exposure to the ad could cause people to experience cognitive dissonance, which they resolve by deciding that non-human animals are incapable meaningfully suffering or feeling).

Given how much WAIC likes models with *age* in them, let’s also explore variation in consumption at different ages. To do this, I’ll simulate American women across the age range seen in our data (predicting out-of-sample doesn’t work too well, obviously) who both did and did not receive the treatment. Maybe some ages are more susceptible to persuasion?

First, let’s define the range of ages we’ll explore (we did this before, too. Here, I’ll leave ages standardized):

age.seq <- seq(from = -1.5, to = 2, by = .1)

Next, let’s simulate individuals across this age range who either did or did not receive the treatment, and calculate mean servings consumed and an 89% percentile interval of simulated consumption about that mean:

d.pred.NT <- data.frame( Treat = 0, gender_id = 1, country_id = 6, Age_s = age.seq ) d.pred.T <- data.frame( Treat = 1, gender_id = 1, country_id = 6, Age_s = age.seq ) sims.a.T <- ensemble(m13cc.p, m12cc.p, m9cc.p, m10cc.p, m7cc.p, m6cc.p, m8cc.p, data = d.pred.T, do_link = F, do_sim = T, n= 5e3)$sim consMean.T <- apply(sims.a.T, 2, mean) consPI.T <- apply(sims.a.T, 2, PI, prob = .89) sims.a.NT <- ensemble(m13cc.p, m12cc.p, m9cc.p, m10cc.p, m7cc.p, m6cc.p, m8cc.p, data = d.pred.NT, do_link = F, do_sim = T, n= 5e3)$sim consMean.NT <- apply(sims.a.NT, 2, mean) consPI.NT <- apply(sims.a.NT, 2, PI, prob = .89)

Now, let’s plot our predictions (means and PIs), with treatment in red and control in black, overlain on the raw, jittered data:

plot( jitter(X2DaysPork) ~ jitter(Age.s) , dcc , col=col.alpha(1,0.1) ) lines( age.seq , consMean.T , col = 2, lwd = 2) shade( consPI.T , age.seq, col = col.alpha(2, 0.1) ) lines( age.seq , consMean.NT , col = 1, lwd = 2) shade( consPI.NT , age.seq) legend(1.3,9.4, c("Control","Treatment"), lty=c(1,1), lwd=c(2.5,2.5),col=c(1, 2)) title("Predicted Effects of Age on Consumption")

Because outcomes are discrete, our PIs are all blocky and gross. But in any case, from this we can once more see that older individuals are the ones for whom exposure to a pro-veg*n advertisement *increases* pork consumption, and younger individuals are the ones for whom the treatment *decreases* pork consumption. This pattern is reflected both in average consumption, and in the percentile interval about that average. Looking at this graph, you might also think, wait, these lines aren’t identical to those seen in the graph for lambda — isn’t the mean outcome of a Poisson process just its rate? But recall that our model’s a zero-inflated Poisson, so our outcomes take into account the fact that a decent chunk of the simulated population (about 33%) doesn’t eat pork, whatever their reasons, thus bringing average pork consumption down.

But perhaps **chicken/bird consumption** fares “better”? I fit the same exact models as above with # of bird servings as my outcome variable (making some basic checks to ensure MCMC health by visually inspecting trace plots, ensuring the GR diagnostic was 1, or at worst 1.01 in a couple cases, and that my effective number of parameters was always at least several hundred), and got the following WAIC weights (I won’t reproduce my code here, since it’s essentially identical to what’s above and this document’s already really long without the redundancy):

So a different set of models is supported in the chicken data — indeed, most of the weight is held by models without any treatment variables (denoted .nt in the model names), suggesting that treatment has little effect on this outcome. Nevertheless, some of the supported models do have treatment effects (including the 2nd best supported model), though generally on *lambda* and not *p*.

The MAP coefficient estimates for those models with non-zero weight were:

So we see a small positive main effect of the treatment on *p, *and a mixed, but mostly negative main effect of the treatment on lambda. Importantly, however, the interaction of treatment and woman (lGT[1]) is positive, and larger than the main effect of treatment. Since interpreting results from a table of MAP estimates is a total crapshoot, let’s do as we did before and simulate *p*s and *lambda*s as well as outcomes for average-aged American women:

As expected, it looks like there’s almost no effect of treatment on *p, *since very few of the models we averaged had coefficients representing that relationship. Like with pork, it looks like the treatment increases bird consumption — by how much?

Probably by nothing, since most of our best-supported models had no effect of the treatment on lambda, since they had no effect of treatment on anything. Most of the probability’s centered around 0, indicating that the model’s fairly confident that the treatment had very little effect on non-abstainers bird eating habits.

Looking at how *p* and *lambda *change with respect to bird consumption for the treatment and control groups across different ages, we get (again, with solid lines representing means and shaded regions 89% HPDIs):

So it doesn’t look like there’s any real relationship with age, unlike in the pork models.

Simulating individuals from our inferred models, we get the following histogram of consumption:

With an average difference in consumption between treatment and control of 0.047 servings per 2 days (that is to say, exposure to treatment on average seems to increase bird consumption by 0.047 servings per 2 days, but there’s a fairly good chance it doesn’t do anything at all. So maybe the bird story is not as bleak as the pig story, in that it could well be that exposure to the vegetarian ads has no effect o young American women’s dietary habits).

We can look at the interaction between age and treatment on simulated outcomes, too, with the same sort of plot as before:

Though there’s a clear (though uncertain) downward trend with age and individuals consume fewer birds the older they are (as we saw for predicted *lambdas*), it still doesn’t look like treatment has much, if any effect (since the red line is above the black one, treatment looks like it may increase bird consumption by a small amount for essentially all ages, similar to what we saw when inspecting simulation output earlier).

Ok, so that’s enough of that! Let’s look at some of the attitude questions. For this, I’d like to primarily look at the **LessMeatGood** outcome, which represents individuals’ attitudes (on a typical 5-point Likert scale) to the statement: “Eating less meat, chicken and fish (in total) is the right thing to do.”

We’ll fit a variety of ordered logit models to the data to explore the effects of treatment on attitude. To do this, we’ll be working on the scale of log cumulative odds, so that we can fit linear models on a cumulative probability scale between 0 and 1 but, with a link function, have predictor variables that can “contribute” an effect above 1 and below 0. Fundamentally, though, our likelihood is the ordered distribution, which is parameterized with a vector of numbers between 0 and 1 describing the probability of observing a response less than or equal to some corresponding maximum response value. We’ll incorporate a linear model in addition to intercepts (/cutpoints) for each categorical outcome, which we’ll assign to *phi *(which will be related to the vector of probabilities with a link function. It might be helpful for intuition-building purposes to play around with the rordlogit(100, phi = 0, a = c(1,2,3,4,5)) function in R, perhaps plotting simplehist(rordlogit(100, phi = 0, a = c(1,2,3,4,5))) with different values for *phi* and *a*. A picture might serve well here, too (note, I’ve recoded the responses as numbers, such that “Strongly Disagree” becomes a 1, “Somewhat Disagree” a 2, and so on until “Strongly Agree” becomes a 5).:

#recode responses as numbers dcc$LMG_i <- as.numeric(dcc$LessMeatGood) #doing some plotting par(mfrow = c(1,3)) simplehist(dcc$LMG_i, xlim = c(1,5), xlab = "response") title("Histogram of Responses") # discrete proportion of each response value pr_k <- table( dcc$LMG_i ) / nrow(dcc) # cumsum converts to cumulative proportions cum_pr_k <- cumsum( pr_k ) # plot plot( 1:5 , cum_pr_k , type="b" , xlab="response" , ylab="cumulative proportion" , ylim=c(0,1) ) title("Cumulative Proportion of Responses") logit <- function(x) log(x/(1-x)) # convenience function 11.4 lco <- logit( cum_pr_k ) plot( 1:5 , lco , type="b" , xlab="response" , ylab="log cumulative odds") title("Log Cumulative Odds of Responses")

Let’s fit some models! Starting with the following, at least for now (I’m getting impatient waiting for models to run…). Priors chosen were pretty flat, since there’s a lot of data here from which the model can learn:

# Fitting Models # #simplest model, just intercepts m1 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- 0, #note: link function already embedded in dordlogit likelihood cutpoints ~ dnorm(0,5) ) , data= list(LMG_i = dcc$LMG_i) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) #start the chain off somewhere decent precis(m1, depth = 2)

#effect of treatment on phi m2 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT, cutpoints ~ dnorm(0,5), pT ~ dnorm(0,5) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m2, depth = 2)

m3 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + Age_s * pA, cutpoints ~ dnorm(0,5), c(pT, pA) ~ dnorm(0,5) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m3, depth = 2)

m4 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + Age_s * pA + pG[gender_id], cutpoints ~ dnorm(0,5), c(pT, pA) ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s, gender_id = dcc$gender_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m4, depth = 2)

m5 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + Age_s * pA + pG[gender_id] + pC[country_id], cutpoints ~ dnorm(0,5), c(pT, pA) ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m5, depth = 2)

m6 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + pG[gender_id] + pC[country_id], cutpoints ~ dnorm(0,5), pT ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m6, depth = 2)

m7 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + Treat * Age_s * pTA + Age_s * pA + pG[gender_id] + pC[country_id], cutpoints ~ dnorm(0,5), c(pT, pA, pTA) ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m7, depth = 2)

m8 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + Treat * Age_s * pTA + Age_s * pA, cutpoints ~ dnorm(0,5), c(pT, pA, pTA) ~ dnorm(0,5) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m8, depth = 2)

## No effect of treatment ## m9n <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Age_s * pA, cutpoints ~ dnorm(0,5), pA ~ dnorm(0,5) ) , data= list(LMG_i = dcc$LMG_i, Age_s = dcc$Age_s) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m9n, depth = 2)

m10n <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Age_s * pA + pG[gender_id], cutpoints ~ dnorm(0,5), pA ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Age_s = dcc$Age_s, gender_id = dcc$gender_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m10n, depth = 2)

m11n <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Age_s * pA + pG[gender_id] + pC[country_id], cutpoints ~ dnorm(0,5), pA ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Age_s = dcc$Age_s, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m11n, depth = 2)

m12n <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- pG[gender_id] + pC[country_id], cutpoints ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m12n, depth = 2)

m13 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + Treat * Age_s * pTA + Age_s * pA + pG[gender_id] + pC[country_id] + Treat*pGT[gender_id], cutpoints ~ dnorm(0,5), c(pT, pA, pTA) ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pGT[gender_id] ~ dnorm(0,sigGTp), sigGTp ~ dcauchy(0,1), pC[country_id] ~ dnorm(0,sigCp), sigCp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m13, depth = 2)

m14 <- map2stan( alist( LMG_i ~ dordlogit( phi , cutpoints ), phi <- Treat * pT + pG[gender_id] + Treat*pGT[gender_id], cutpoints ~ dnorm(0,5), pT ~ dnorm(0,5), pG[gender_id] ~ dnorm(0,sigGp), sigGp ~ dcauchy(0,1), pGT[gender_id] ~ dnorm(0,sigGTp), sigGTp ~ dcauchy(0,1) ) , data= list(LMG_i = dcc$LMG_i, Treat = dcc$Treat, Age_s = dcc$Age_s, gender_id = dcc$gender_id, country_id = dcc$country_id) , iter = 25000, warmup = 5000, start=list(cutpoints = c(-1.88, -.76, .62, 1.75)) ) precis(m14, depth = 2)

And then we can compare them (after checking MCMC health) with WAIC, like before:

compare(m1, m2, m3, m4, m5, m6, m7, m8, m9n, m10n, m11n, m12n, m13, m14)

So it looks like models with a treatment predictor aren’t doing as hot, since the top two models lack it (denoted with an *n* in the model name). This suggests that including a treatment predictor results in overfitting, and that treatment does not have much influence on the data (including it will, of course, improve within-sample deviance, but we’re interested in estimated out-of-sample deviance). Still, a non-insignificant portion of WAIC weight is held by treatment-containing models, so let’s simulate from them. But first, we can quickly examine estimates from the models with non-zero WAIC weight:

coeftab(m12n, m11n, m6, m5, m10n, m7, m13, m14, m4)

Estimated intercepts are all pretty similar, and all ordered in the right way (increasing from one to four). The main effect of treatment is alternatively both positive and negative, but we can’t actually tell much from this because we can’t see the uncertainty about each point or immediately grasp the effects of tricky interaction terms. And I dunno about you, but I can’t intuitively grasp values on a log cumulative odds scale. On to simulation!

…I’m having some trouble getting link(), sim(), and ensemble() to work with these models. Dunno why. So let’s do it manually, instead, starting with the treatment model that received the most WAIC weight, m6 (a piddling 15%, sure, so the model isn’t too well supported, but whatever! We’re interested in the effects of treatment, and that’s the best we’ve got! For what it’s worth, it gets 45% of the WAIC weight when we exclude the no-treatment-effects models from consideration).

Run the following code to extract 200 samples from the posterior and loop over predictions from m6 across the variation in treatment:

#extract samples from the posterior nsamp <- 200 #number of samples to extract from the posterior and subsequently plot post <- extract.samples( m6 , n = nsamp) plot( 1 , 1 , type="n" , xlab="Treatment" , ylab="probability" , xlim=c(0,1) , ylim=c(0,1) , xaxp=c(0,1,1) , yaxp=c(0,1,2) ) Treat <- 0:1 #range of treatment values to iterate over gender_id <- 1 #value for gender (i.e. woman) country_id <- 6 #value for country (i.e. USA) for ( s in 1:nsamp ) { p <- lapply(1:length(post), function (x) if(length(dim(post[[x]])) == 2) post[[x]][s,] else post[[x]][s]) names(p) <- names(post) ak <- as.numeric(p$cutpoints[1:4]) phi <- p$pT*Treat + p$pG[gender_id] + p$pC[country_id] pk <- pordlogit( 1:4 , a=ak , phi=phi ) for ( i in 1:4 ) lines( Treat , pk[,i] , col=col.alpha(rangi2,0.1) ) } title("Effect of Treatment on Cumulative Probability of Response") abline(h = 1, lty = 2) abline(h = 0, lty = 2) text(x = .8, y= .05, labels = "Strongly Disagree") text(x = .8, y= .2, labels = "Somewhat Disagree") text(x = .8, y= .5, labels = "Neither Agree nor Disagree") text(x = .8, y= .75, labels = "Somewhat Agree") text(x = .8, y= .95, labels = "Strongly Agree")

It looks like the model’s pretty confident that treatment has very little effect on disposition with respect to the statement “Eating less meat, chicken and fish (in total) is the right thing to do.” So even models with a treatment effect suggest that treatment has very little effect, at least with respect to an American woman’s beliefs regarding the moral correctness of consuming animals. The lines are mostly flat, and you can examine their distribution (across a greater proportion of the posterior) with the following code (where the number outside of the parentheses indicates the coded response and the number inside the parentheses is the proportion of negative slopes):

nsamp <- 5000 #number of samples to extract from the posterior and subsequently plot predictions of post <- extract.samples( m6 , n = nsamp) Treat <- 0:1 #range of treatment values to iterate over gender_id <- 1 #value for gender country_id <- 6 #value for country (i.e. USA) diffs <- matrix(nrow = 4, ncol = nsamp) for ( s in 1:nsamp ) { p <- lapply(1:length(post), function (x) if(length(dim(post[[x]])) == 2) post[[x]][s,] else post[[x]][s]) names(p) <- names(post) ak <- as.numeric(p$cutpoints[1:4]) phi <- p$pT*Treat + p$pG[gender_id] + p$pC[country_id] pk <- pordlogit( 1:4 , a=ak , phi=phi ) for ( i in 1:4 ) diffs[i,s] <- pk[2,i] - pk[1,i] } dens(diffs[1,], lwd = 2) text(x = mean(diffs[1,]), y = 2, labels = paste0("1 (", sum(diffs[1,] < 0) / length(diffs[1,]), ")"), lwd = 2) for(i in 2:4){ dens(diffs[i,], add = T, col = i, lwd = 2) text(x = mean(diffs[i,]), y = i*2, labels = paste0(i, " (", sum(diffs[i,] < 0) / length(diffs[i,]), ")"), col = i, lwd = 2) }

So our eyes did not deceive us — the slopes were pretty evenly distributed around zero, but with slightly more slopes being negative than positive, and there being more uncertainty with respect to the slopes of 2 and 3 than 1 and 4, as visible in the previous graph.

I also just noticed something — *meat* is usually defined broadly as animal flesh or narrowly as muscle tissue, yet this statement specifies *meat, chicken, **and fish*, suggesting to me that chicken and fish are not made of meat. I wonder if this confused any of the survey takers?

Anyway, Model 6 got the most WAIC weight when looking at only models with treatment effects, but Model 5 was hot on its tail, receiving 26% of the weight compared to m6.

compare(m1, m2, m3, m4, m5, m6, m7, m8, m13, m14)

Model 5 also had an effect for age, where Model 6 did not. So by looking at Model 5, we can test the hypothesis from earlier, that older individuals might have resolved in favor of eating meat after being exposed to treatment, and that this was driving the unfortunate trend seen there. With some small modifications to our earlier code, let’s plot predictions from the posterior to test it out!

#testing the effects of age given treatment nsamp <- 200 #number of samples to extract from the posterior and subsequently plot post <- extract.samples( m5 , n = nsamp) age.seq <- seq(from = -1.5, to = 2, by = .2) par(mfrow = c(6,3)) for(ages in 1:length(age.seq)){ plot( 1 , 1 , type="n" , xlab="Treatment" , ylab="probability" , xlim=c(0,1) , ylim=c(0,1) , xaxp=c(0,1,1) , yaxp=c(0,1,2) ) Treat <- 0:1 #range of treatment values to iterate over gender_id <- 1 #value for gender (i.e. woman) country_id <- 6 #value for country (i.e. USA) Age_s <- age.seq[ages] for ( s in 1:nsamp ) { p <- lapply(1:length(post), function (x) if(length(dim(post[[x]])) == 2) post[[x]][s,] else post[[x]][s]) names(p) <- names(post) ak <- as.numeric(p$cutpoints[1:4]) phi <- p$pT*Treat + p$pG[gender_id] + p$pC[country_id] + p$pA*Age_s pk <- pordlogit( 1:4 , a=ak , phi=phi ) for ( i in 1:4 ) lines( Treat , pk[,i] , col=col.alpha(rangi2,0.1) ) } title(paste0("Age_s = ", Age_s)) abline(h = 1, lty = 2) abline(h = 0, lty = 2) text(x = .8, y= .05, labels = "StD") text(x = .8, y= .2, labels = "SoD") text(x = .8, y= .5, labels = "NAND") text(x = .8, y= .75, labels = "SoA") text(x = .8, y= .95, labels = "StA") }

Not much to be gleaned from here, besides perhaps the insight that predictions for older ages is progressively more and more uncertain (probably just due to variation in sampling intensity). If I cross my eyes I can maaaaaaybe see a downward trend in some of them, overall (indicating that people shift their answers from the *disagree* end of the spectrum to the *agree* end), but (keeping in mind that models with treatment effects are not well supported) it’d be better to examine the distribution of differences explicitly, with similar code as before:

nsamp <- 5000 #number of samples to extract from the posterior and subsequently plot post <- extract.samples( m5 , n = nsamp) age.seq <- seq(from = -1.5, to = 2, by = .2) par(mfrow = c(6,3)) for(ages in 1:length(age.seq)){ Treat <- 0:1 #range of treatment values to iterate over gender_id <- 1 #value for gender (i.e. woman) country_id <- 6 #value for country (i.e. USA) Age_s <- age.seq[ages] diffs <- matrix(nrow = 4, ncol = nsamp) for ( s in 1:nsamp ) { p <- lapply(1:length(post), function (x) if(length(dim(post[[x]])) == 2) post[[x]][s,] else post[[x]][s]) names(p) <- names(post) ak <- as.numeric(p$cutpoints[1:4]) phi <- p$pT*Treat + p$pG[gender_id] + p$pC[country_id] + p$pA*Age_s pk <- pordlogit( 1:4 , a=ak , phi=phi ) for ( i in 1:4 ) diffs[i,s] <- pk[2,i] - pk[1,i] } dens(diffs[1,], lwd = 2) title(paste0("Age_s = ", Age_s)) text(x = mean(diffs[1,]), y = 2, labels = paste0("1 (", sum(diffs[1,] < 0) / length(diffs[1,]), ")"), lwd = 2) for(i in 2:4){ dens(diffs[i,], add = T, col = i, lwd = 2) text(x = mean(diffs[i,]), y = i*2, labels = paste0("1 (", sum(diffs[i,] < 0) / length(diffs[i,]), ")"), col = i, lwd = 2) } }

So slopes appear to still be distributed roughly around zero — indeed, the effect of treatment is almost identical to the one we saw before. Looking at how the distribution shifts with age, it seems 4 and 3 gets more uncertain at higher ages, and 1 and 2 get more certain. Not by much, though.

It might be easier to examine predictions not as a series of plots but as an animation. Changing our code from above slightly (and drawing a few more samples from the posterior, setting a smaller interval for the differences between ages), we get:

Of course, model 5 did *not* have an explicit interaction between treatment and age, just a main effect for both; recall the linear equation for phi:

phi <- Treat * pT + Age_s * pA + pG[gender_id] + pC[country_id] #m5

This isn’t to say that Age and Treatment couldn’t interact — draws in parameter values for each could have been correlated, and I think floor and ceiling effects make all parameters interact with each other implicitly. Nevertheless, we might want to examine a model with an explicit interaction term between treatment and age.

Which model 7, our third best WAIC-approved model (that included treatment terms), did have. We’re sorta skimming the bottom of the barrel here — m7 received a whopping 4% of the WAIC weight when all the models were considered. But whatever — let’s probe what this model has to say about variation in the effects of treatment across different ages (examining the marginal posterior of pTA for model 7 gives a mean of -0.06, but with a high sd, 0.09, representing a fairly large amount of uncertainty about that estimate):

AHA! There *is* a ton of uncertainty about these estimates, and we *did *have to use a rather poorly (thought non-trivially) supported model, but we can totally see the trend I predicted earlier — exposure to the survey convinces *younger* people that “eating less meat, chicken and fish (in total) is the right thing to do”, by shifting their opinions from the disagree to the agree end of the spectrum, and conversely increases the probabilities with which older people disagree with the statement (this manifests as a shift from downward to upward slopes as our simulated individuals get older). From the animated version, it looks like most of the effect is coming from the “control” opinions shifting as individuals age. When we examine the distributions of slopes explicitly, we see the same story:

Note how the entire distributions shift rightward as you age towards more positive values (and the numbers in parentheses — indicating the proportion of negative slopes — get smaller). It’s not a *huge* trend — the advertisement isn’t turning younger people into animal activists or convincing older people to slap “People Eating Tasty Animals” stickers on their cars, but it does seem to be fairly clear. It’s a shame we had to examine a model with such a low WAIC weight to see it: an explicit age and treatment interaction is probably overfit, so if we actually averaged predictions across model uncertainty the effect would be even more minuscule. But if you have high prior confidence in the model, perhaps you’d be willing to cast caution to the wind and believe what m7 has to say about how the effects of exposure to pro-veg*n advertisement vary as American women age.

From the animation, we also see that the model predicts a greater range of simulated slopes at extreme ages, probably a reflection of the fewer data points available at those ages. I’m a little concerned about sensitivity to my choice of priors, here (the one I chose for the Treatment/Age interaction was really flat — norm(0,5) — and I *think* “uninformative” and symmetric, but I’m not as familiar with the ordered distribution as I’d like to be…), though I think, in this case, they’re completely overwhelmed by the likelihood and information present in the data. If we examine both the prior and posterior distributions of pTA, we can see how much and how confidently the model’s learned, and be at least a little assured until a full sensitivity analysis can be performed later:

[TO BE CONTINUED]

(Header Image Source: Wikimedia Commons)

## One thought on “MFA Ads Analysis”