Different output when running lmer and summarise

Question

I am doing some easy linear mixed effects model as below:

lmer( as.numeric(o.m_20_r_coded)~  site + (1 | village), o.m_20_r_coded_dat) %>% summary()

Random effects:
 Groups   Name        Variance Std.Dev.
 village  (Intercept) 0.3285   0.5732  
 Residual             1.2504   1.1182  
Number of obs: 3580, groups:  village, 59

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   3.6879     0.1194 116.6521  30.887   <2e-16 ***
siteWASEP    -0.1377     0.1106 465.1790  -1.245    0.214  

But when I run descriptives, I get very different results:

o.m_20_r_coded_dat %>% group_by(site) %>% summarise(as.numeric(mean(o.m_20_r_coded,na.rm=T)))

control 3.134864            
WASEP   3.592092    

Why is there such a huge discrepancy? Based on model00, the intercept for WASEP is smaller than the intercept for control, whereas the descriptives show the opposite! What is going on?

There are lots of NAs in the data.