Your residual structure is totally expected with this model specification and an indication of an ill-specified model. What you basically are trying to do is to fit a linear line through points that can only take values of 0 and 1 on the $y$-axis.
Let's look at a simple example with arbitrarily generated variables:
#-----------------------------------------------------------------------------
# Generate random data for logistic regression
#-----------------------------------------------------------------------------
set.seed(123)
x <- rnorm(1000)
z <- 1 + 2*x
pr <- 1/(1+exp(-z))
y <- rbinom(1000,1, pr)
#-----------------------------------------------------------------------------
# Plot the data
#-----------------------------------------------------------------------------
par(bg="white", cex=1.2)
plot(y~x, las=1, ylim=c(-0.1, 1.3))
#-----------------------------------------------------------------------------
# Fit a linear regression (nonsensical) and plot the fit
#-----------------------------------------------------------------------------
linear.mod <- lm(y~x)
segments(-2.32146, 0, 1.24196, 1, col="steelblue", lwd=2)
segments(1.24196, 1, 100, 28.71447, col="red", lwd=2)
segments(-100, -27.41153, -2.32146, 0, col="red", lwd=2)
As you can see, a linear line is fitted through the data. One problem of this is that the line predicts outcomes that are outside the interval $[0,1]$ (illustrated by the red lines outside that interval). Let's have a look at the residuals:
#-----------------------------------------------------------------------------
# Add the residual lines
#-----------------------------------------------------------------------------
x.y0 <- sample(which(y==0), 50, replace=F)
x.y1 <- sample(which(y==1), 50, replace=F)
pre <- predict(linear.mod)
segments(x[x.y0], y[x.y0], x[x.y0], pre[x.y0], col="red", lwd=2)
points(x[x.y0], y[x.y0], pch=16, col="red", las=1)
segments(x[x.y1], y[x.y1], x[x.y1], pre[x.y1], col="blue", lwd=2)
points(x[x.y1], y[x.y1], pch=16, col="blue", las=1)
I randomly picked some values to show the pattern. The red and blue lines are depicting the residuals, which is the difference between the predicted value of the line and the actual observed value (red and blue dots). The blue lines correspond to the residuals where $y=1$ whereas the red residuals correspond to the situation where $y=0$. Because the outcome can only be either 0 or 1, the residuals are simply the distances between the regression line and either 0 or 1. The residuals take exactly the form that you see in your data:
#-----------------------------------------------------------------------------
# Plot the residuals
#-----------------------------------------------------------------------------
res.linear <- residuals(linear.mod, type="response")
par(bg="white", cex=1.2)
plot(predict(linear.mod)[y==0], res.linear[y==0], las=1,
xlab="Fitted values", ylab = "Residuals",
ylim = max(abs(res.linear))*c(-1,1), xlim=c(-0.4, 1.6), col="red")
points(predict(linear.mod)[y==1], res.linear[y==1], col="blue")
abline(h = 0, lty = 2)
The colors correspond to the residuals shown above: the blue dots are the residuals where $y=1$ and the red dots are the residuals where $y=0$. In normal linear regression, the residuals are assumed to be approximately normally distributed. But in this case, the residuals can hardly be normal. They are binomial.
We need a transformation that transformes the probability, which is bound within $[0,1]$ into a variable that ranges over $(-\infty, \infty)$. One such transformation is the logit (this is not the only possibility: we could also use probit or the complementary log-log function). Let's fit a logistic regression with a logit-link and again plot the binned residuals (explained on page 97 by Gelman and Hill (2007)). Plotting the raw residuals vs. fitted values are generally not useful after logistic regression:
#-----------------------------------------------------------------------------
# Fit a logistic regression
#-----------------------------------------------------------------------------
glm.fit <- glm(y~x, family=binomial(link="logit"))
#-----------------------------------------------------------------------------
# Plot the binned residuals as recommended by Gelman and Hill (2007)
#-----------------------------------------------------------------------------
library(arm)
par(bg="white", cex=1.2, las=1)
binnedplot(predict(glm.fit), resid(glm.fit), cex.pts=1, col.int="black")
The residuals in logistic regression can be define -$~$as in linear regression$~$- as observed minus expected values:
$$
\text{residual}_{i}=y_{i}-\mathrm{E}(y_{i}|X_{i})=y_{i}-\text{logit}^{-1}(X_{i}\beta)
$$
Because the data $y_{i}$ are discrete, so are the residuals. In the plot above, the residuals are binned by dividing the data into categories based on their fitted values, and are then plotted against the average residual versus the average fitted value for each category (bin). The lines indicate $\pm2$ standard-error bounds, within which one we would expect about 95% of the binned residuals to fall, under the assumption that the model is true.
So the remedy for your immediate problem is to fit a mixed effects logistic regression by typing:
model <- glmer(error~is_frisian*condition*person+(1|subject_id),
data=output, family="binomial")
For a good introduction to mixed effects logistic regression in R, see here. For a good overview of diagnostics in linear and generalized linear models, see here.
Best Answer
Maybe I am stating some obvious but at the following link http://glmm.wikidot.com/faq under the section How do I compute a coefficient of determination (R2), or an analogue, for (G)LMMs? provides a quite good treatment of the issue. Just to quote some of the final lines from there to avoid possible link rot and provide the reference they present:
A rather helpful reference of the subject seems to be Xu's: Measuring explained variation in linear mixed effects models (2003) where the author computes the residual variance of the full model against the residual variance of a (fixed) intercept-only null model (in R syntax):
1-var(residuals(m))/(var(model.response(model.frame(m)))