Can I calculate adjusted means from ANCOVA if data is unbalanced? I have a 10k lines dataframes on which I want to perform ANCOVA so I can get adjusted means.
Please note that I've never done this before so I jump from a tutorial to another, but I still want to make it the right way.
So my model is like Y ~ X * sex, with 


*

*Y the dependant variable (continuous)

*X the continuous independant variable 

*sex the discrete independant variable (here the sex)


Reading this tutorial, I could calculate the Y mean adjusted on X for each sex :
model = aov(Y~sex*X, data=x)
data.predict = data.frame(sex=c("Male", "Female"), X=mean(x$X, na.rm=T))
data.frame(data.predict, Y=predict(model, data.predict))

This gives realistic results, but I realized that anova(aov(Y~sex*X, data=x)) and anova(aov(Y~X*sex, data=x)) give very different results. The calculated means are the same with both models though.
Reading the EdM answer in the question https://stats.stackexchange.com/a/213358/81974, I tried with the car package and Anova(model, type="III"), and this time both give the same results.
I don't really understand how it could matter, but it seems that my data are unbalanced (the aov help "Note" says that it could be misleading).
Knowing this, are the previously calculated adjusted means still usable ?
 A: (Note:  this answer will be mostly about using R, but hopefully the discussion of statistical concepts will keep it on-topic for this site.)
1) [EDIT: Comments on unbalanced designs deleted.  See comments below by @rvl and @SalMangiafico]
2) By default, R uses Type I Sums of Squares.  The Anova function allows you to use Type II or Type III.  For Type III you will need to change the contrasts R uses to code dummy variables.
3) For adjusted (marginal or least square) means, you can use the emmeans package. 
The following is a short example in R.  The marginal means of Y for Male and Female are similar, but their arithmetic means are quite different due to the influence of the variable X.
[EDIT: I changed the code below to produce the e.m. means for the X:Sex interaction rather than for the main effect of Sex, as per comment by @rvl. For this model this change doesn't affect the result.] 
if(!require(car)){install.packages("car")}
if(!require(emmeans)){install.packages("emmeans")}

set.seed(5)

X   = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
Y   = X + rnorm(20, 0, 1)
Sex = c(rep("Male", 10), rep("Female", 10))

Data = data.frame(Y, X, Sex)

mean(Y[Sex=="Male"])

   ### 
   ### [1] 5.421148
   ###

mean(Y[Sex=="Female"])

   ###
   ### [1] 15.01774
   ###

model = lm(Y ~ X + Sex + X:Sex, data = Data)

library(car)

Anova(model)

   ### Anova Table (Type II tests)
   ###
   ### Response: Y
   ###            Sum Sq Df  F value   Pr(>F)    
   ### X         149.664  1 158.3660 1.03e-09 ***
   ### Sex         0.007  1   0.0069   0.9347    
   ### X:Sex       0.104  1   0.1096   0.7449    
   ### Residuals  15.121 16  

library(emmeans)

marginal = emmeans(model, ~ X:Sex)

marginal

   ###  Sex      emmean        SE df lower.CL upper.CL
   ###  Female 10.38104 0.6171581 16 9.072726 11.68936
   ###  Male   10.30839 0.6171581 16 9.000076 11.61671 

