<p>My impression is that your problems come from mixing various approaches for various aspects (repeated measurements/correlation vs. heteroscedasticity) that cannot be mixed so easily. Instead of using random effects you might also consider fixed effects, or instead of only adjusting the inference for heteroscedasticity you might consider a Gaussian model and model both mean and variance, etc. For me, it's hard to say what is the best route forward here. Hence, I only comment on some aspects regarding the <code>sandwich</code> package:</p>

<p>The <code>sandwich</code> package is <em>not</em> limited to <code>lm</code>/<code>glm</code> only but it is in principle object-oriented, see <code>vignette("sandwich-OOP", package = "sandwich")</code> (also published as <a href="https://dx.doi.org/10.18637/jss.v016.i09" rel="nofollow noreferrer">doi:10.18637/jss.v016.i09</a>.</p>

<p>There are suitable methods for a wide variety of packages/models but not
for <code>nlme</code> or <code>lme4</code>. The reason is that it's not so obvious for which mixed-effects models the usual sandwich trick actually works. (Disclaimer: But I'm no expert in mixed-effects modeling.)</p>

<p>However, for <code>lme4</code> there is a relatively new package
called <code>merDeriv</code> (<a href="https://CRAN.R-project.org/package=merDeriv" rel="nofollow noreferrer">https://CRAN.R-project.org/package=merDeriv</a>) that
supplies <code>estfun</code> and <code>bread</code> methods so that <code>sandwich</code> covariances can be
computed for <code>lmer</code> output etc. There is also a working paper associated
with that package: <a href="https://arxiv.org/abs/1612.04911" rel="nofollow noreferrer">https://arxiv.org/abs/1612.04911</a></p>
R: mixed model with heteroscedastic data -> only lm function works?
<p><a href="https://stackoverflow.com/questions/26224877/two-fixed-factors-nested-and-crossed-factors-in-r">This question</a> asks the same question, but hasn't been answered. My question relates to how to specify the model with the lm() function and is therefore a programming (not statistical) question.</p>

<p>I have a mixed design (2 repeated and 1 independent predictors). Participants were first primed into group A or B (this is the independent predictor) and then they rated how much they liked 4 different statements (these are the two repeated predictors). 
There are many great online resources how to model this data. However, my data is heterscedastic. So I like to use heteroscedastic-consistent covariance matrices. <a href="https://www.ncbi.nlm.nih.gov/pubmed/24015776" rel="nofollow noreferrer">This paper</a> explains it well. The <a href="https://cran.r-project.org/web/packages/sandwich/sandwich.pdf" rel="nofollow noreferrer">sandwich</a> and <a href="https://cran.r-project.org/web/packages/lmtest/lmtest.pdf" rel="nofollow noreferrer">lmtest</a> packages are great. <a href="https://stats.stackexchange.com/questions/117052/replicating-statas-robust-option-in-r">Here</a> is a good explanation how to do it for a indpendent design in R with lm(y ~ x). </p>

<p>It seems that I have use lm, else it wont work? </p>

<p>Here is the code for a regression model assuming that all variances are equal (which they are not as Levene's test comes back significant). </p>

<pre><code>fit3 &lt;- nlme:::lme(DV ~ repeatedIV1*repeatedIV2*independentIV1, random =  ~1|participants, df) ##works fine
</code></pre>

<p>Here is the code for an indepedent model correcting for heteroscedasticity, which works.</p>

<pre><code>fit3 &lt;- lm(DV ~ independentIV1)
library(sandwich)
vcovHC(fit3, type = 'HC4', sandwich = F)
library(lmtest)
coef(fit3, vcov. = vcovHC, type = 'HC4')
</code></pre>

<p>So my question really is, how to specify my model with lm? 
Alternative approaches in R how to fit my model accounting for heteroscedasticity are welcome too!</p>

<p>Thanks a lot!!! </p>