This work considers a multifactor linear mixed model under heteroscedasticity in random-effect factors and the skew-normal errors for modeling the correlated datasets. We implement an expectation–maximization (EM) algorithm to achieve the maximum likelihood estimates using conditional distributions of the skew-normal distribution. The EM algorithm is also implemented to extend the local influence approach under three model perturbation schemes in this model. Furthermore, a Monte Carlo simulation is conducted to evaluate the efficiency of the estimators. Finally, a real data set is used to make an illustrative comparison among the following four scenarios: normal/skew-normal errors and heteroscedasticity/homoscedasticity in random-effect factors. The empirical studies show our methodology can improve the estimates when the model errors follow from a skew-normal distribution. In addition, the local influence analysis indicates that our model can decrease the effects of anomalous observations in comparison to normal ones.