assumptions for lmer models subject <- factor(rep(c(1,2,3,4,5,6,7),each=4, times=2))
dep <- c(5,4,9,3,4,4,2,1,10,7,8,7,1,2,1,1,5,10,1,7,3,2,1,4,3,8,7,3,1,1,2,1,15,10,20,11,2,2,1,3,11,12,9,7,2,3,1,2,11,9,8,9,3,4,2,1) 

f1 <- factor(rep(c(rep("Female",times=16),rep("Male",times=12)), times=2))
f2 <- factor(rep(c("day1","day2","day3","day4"),times=14))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2)

m <- lmer(dep ~ f1*f2 + (1|sub), data=data)

I'm trying to understand how I can test the assumptions for mixed models. 
1) In the case of model m, should I look at homogeneity of variances for every combination of f1 and f2 like this plot(resid(m)~fitted(.)|f1:f2) or is it enough to simply do plot(resid(m))?
2) In either case my real model is showing a funnel shape, is this too problematic? What would you do in that case?
3) How can I check the linearity assumption for each categorical variable in R?
4) Apart from homogeneity of variances and normality of residuals, are there any other important assumptions that you think I should be aware of?
 A: The commonly quoted assumptions (or "conditions" as I prefer to call some of them) of linear mixed effects models are:


*

*Linearity of the predictors. This can be checked by plotting the residuals against the response and looking for any systematic shape, and by including non-linear terms (or splines) and comparing the model fit. Very often this will not be an issue, and if it is, then including non-linear terms (such as log, exp or polynomials) in the linear predictor may be sufficient. More importantly, substantive domain knowledge should inform whether the linearity condition is justified. For example, in some domains such as pharmacokinetics, we already know, based on rigorous theory and experimentation, that a linear model is not appropriate in some cases.

*The residuals have constant variance. This can be checked with a plot of residuals against fitted values - there should be no pattern/trend.

*The residuals are independent. This can be checked by plotting residuals against covariates - especially time-varying or spatial covariates. There should not be any systematic pattern

*The residuals are normally distributed. This can be checked in many ways, such as a Q-Q plot and a simple histogram. Statistical tests, such as Anderson-Darling and Kolmogorov–Smirnov are also possible.
Note that "residuals" above refer to both the unit-level residuals (often called "errors") and the random effects. For the random effects, this can be  problematic where only a small number of groups/clusters exist in the sample. In the simulated example given in the OP there are 7 clusters. There are lots of rules of thumb to inform a sufficient number of clusters, and 7 is generally thought to be difficult to draw any conclusions.
It is mentioned in the OP that their actual model exhibits a funneled shape plot of residuals vs fitted values. This indicates heteroskadasticity. One way to proceed is to consider a transformations of variables - ideally this should be informed by expert domain knowledge. With this in mind, Box-Cox transformations may be useful.
A: *

*Existence of variance: Do not need to check, in practice, it is always true.


*Linearity: Do not need to check, because your covariates are categorical.


*Homogeneity: Need to Check by plotting residuals vs predicted values.


*Normality of error term: need to check by histogram, QQplot of residuals, even Kolmogorov-Smirnov test.


*Normality of random effect: Get the estimate of random effect (in your case random intercepts), and check them as check the residual. But it is not efficient because you just have 7 random intercepts.
Another assumption is the independent between subjects. No test, based on your judgement. Subject specific random intercept means the correlation between the response variable from the same subject are the same.
