Within-and-between subject design I am planning an experiment where I test whether having a picture affects whether users prefer Coke to Pepsi. (The research question is actually different, but I'm using Coke and Pepsi to simplify things)
Both Coke and Pepsi come in cans and bottles, and maybe different can/bottle designs can interact with the presence of a photo. So I intend to recruit 100 subjects. They will be randomly divided into two groups. Each of whom will undergo two trials:
Group A:
Trial 1: Text question: “Would you rather have a can of Pepsi or a can of Coke?”
Trial 2: Text question: “Would you rather have a bottle of Pepsi or a bottle of Coke?”
Group B:
Trial 1: Same as Group A trial 1, but with pictures of Coke and Pepsi cans.
Trial 2: Same as Group A trial 2, but with pictures of Coke and Pepsi bottles.
I then intend to run the following regression:
$ChoosePepsi_{it}=β_0+β_1 Photo_{it}+β_2 Bottle_{it}+β_3 Photo_{it}Bottle_{it}+ϵ_{it}$
where $i$ is a user and $t$ is a trial. $Photo$ and $Bottle$ are dummy variables, equal to 1 if a photo or bottle was displayed in a trial. $ChoosePepsi$ is equal to 1 if the user chose Pepsi.
$\beta_1$ and $\beta_3$ will tell me whether a photo has an effect on user's preferences.
Questions:


*

*Is this the correct way of running the regression?

*Should I use clustered standard errors? If so, should I cluster by subject?
 A: You can handle this repeated-measures design with multilevel modeling. You have two levels to your data:


*

*Observation. This is where your DV (ChoosePepsi) is measured. This is also where the IV for bottle or can (Bottle) is manipulated.

*Person. This is where the group IV comes in of either text or photo (Photo).
The model is:
Level 1: $y_{op} = \beta_{0p} + \beta_{1p}Bottle_{op} + \epsilon_{op}$
Level 2: $\beta_{0p} = \gamma_{00} + \gamma_{01}Photo_{p} + u_{0p}$ AND $\beta_{1p} = \gamma_{10} + \gamma_{11}Photo_{p} + u_{1p}$
Where $o$ refers to observation and $p$ refers to person. You can see that the dependent variable $y$ varies at both the observation and person (as it has a subscript for both). The two $\beta$s on Level 1 are defined at Level 2, but $\epsilon_{op}$ represents how much someone differs from their expected value at that specific trial. As with any other regression, you will estimate the variance of this.
$\beta_{0p}$ is an intercept that varies at Level 2. It is made up of the average overall intercept $\gamma_{00}$ and the influence of $Photo$ condition ($\gamma_{01}$), which becomes the main effect of $Photo$. The $u_{0p}$ represents how much someone deviates from their predicted intercept. The variance of this will give you the variance of the random intercept.
$\beta_{1p}$ is the effect of $Bottle$ and is predicted by the average effect of $Bottle$, $\gamma_{10}$, which becomes the main effect of $Bottle$. The slope of $Bottle$ is also predicted by $\gamma_{11}$, which will become the cross-level interaction between $Bottle$ and $Photo$.
(All you have to do is plug-in the Level 2 equations into the Level 1 equation to see this). 
Finally, the $u_{1p}$ term is how much someone differs in their slope of $Bottle$ from the expected value of the slope of $Bottle$. The variance of this will tell you the variance of the random slope.
You want to get your data into long format (in dplyr language, you want to gather it). And then you can use the lme4 and lmerTest packages to test your model:
glmer(ChoosePepsi ~ Photo + Bottle + (1+Bottle|ID), family=binomial, data=yourdata)
Where ID is a factor that represents the participant ID.
You want a binomial model, since it seems like you are looking at a logistic regression (i.e., your outcome is dichotomous). If it is continuous, you don't need that argument.
