Estimating expected values for correlated data using random effects models Statement of the problem: in a study, continuous and dichotomous variables were measured for both eyes for 60 individuals. The researchers needs estimates of expected values (means and proportions) for those measurements for all 60 subjects across bot eyes. In order to do this, the 120 eyes from the 60 subjects must be used to provide a pooled estimate.
The proposed random effects models to achieve this are as follows:
$E(y_{ij})=\mu+\alpha_j+\epsilon_{ij}$
and
$Logit(p_{ij})=\gamma+\omega_j$
Where $\mu$ is the overall mean for a continuous variable $y_{ij}$, $\gamma$ is the overall log odds of the probability for the dichotomous variables, $\alpha_j, \omega_j, \epsilon _{ij}$ are uncorrelated random effects with normal distributions ($\alpha_j \sim N(0,\sigma_{\gamma}), \;\omega_j \sim N(0,\sigma_{\omega}), \; \epsilon_{ij} \sim N(0,\sigma_{\epsilon}), Cov(\alpha_j,\epsilon_{ij})=0$). Index $j$ stands for subject and index $i$ stands for eye nested in the subject.
A more complex nested random effects model could be appropriate, however, for the sake of simplicity will be ignored.
I have made a github project with the data and code in R to do this (https://github.com/nmolanog/glmer_question).
Now I present the main issue of this post: for  dichotomous variables I am observing huge differences in estimates ignoring correlation of eyes nested in subjects vs estimates provided by the random effect models. Those differences are so important, that researchers are questioning and mistrusting the approach and its results. For continuous variables the differences in estimates are almost none existent and (as expected) the main differences are found in confidence intervals, where random effects models provides wider CI's (see figure).

See for example variables M and N, the differences between approaches are huge. In the github repo I explored a nested random effects models for variable K obtaining very similar results to those provided by the simpler random effects model.
How those differences could be explained? Is there any problem with the approach?
Update-Sample code:
###estimate proportion for variable K using glm
mk_glm<-glm(K~1,data = ldf, family = binomial(link = "logit"))
mk_glm_ci<-inv.logit(confint(mk_glm))

##arrange result from glm model
(res_df<-data.frame(method="glm",estimate=inv.logit(mk_glm$coefficients),LCI=mk_glm_ci[1],UCI=mk_glm_ci[2]))

#compare to  raw estimate:
ldf$K%>%table()%>%{.[2]/sum(.)}

###estimate proportion for variable K using glmer model 1
mk_glmer<-glmer(K~1+(1|Id),data = ldf, family = binomial(link = "logit"),control=glmerControl(optimizer = "bobyqa"),nAGQ = 20)
mk_glmer_ci<-confint(mk_glmer)
#add result to res_df
(res_df<-rbind(res_df,data.frame(method="glmer",estimate=inv.logit(fixef(mk_glmer)),LCI=inv.logit(mk_glmer_ci[2,1]),UCI=inv.logit(mk_glmer_ci[2,2]))))

###estimate proportion for variable K using glmer model 2, nested random effects
mk_glmer_2<-glmer(K~1+(1|Id/eye),data = ldf, family = binomial(link = "logit"),control=glmerControl(optimizer = "bobyqa"))
mk_glmer_2_ci<-confint(mk_glmer_2)
(res_df<-rbind(res_df,data.frame(method="glmer2",estimate=inv.logit(fixef(mk_glmer_2)),LCI=inv.logit(mk_glmer_2_ci[3,1]),UCI=inv.logit(mk_glmer_2_ci[3,2]))))

Outuput
             method  estimate       LCI       UCI
(Intercept)     glm 0.7083333 0.6231951 0.7846716
(Intercept)1  glmer 0.9230166 0.7399146 0.9990011
(Intercept)2 glmer2 0.9999539 0.9991883 0.9999995

The dataset and code can be found in https://github.com/nmolanog/glmer_question
 A: The largest variation in the width of your confidence intervals occurs in the estimates for the dichotomous outcome variables, so I will focus mostly on that part of the model.  I will speak to the models for the continuous outcome variables at the end.  The phenomenon you are observing is fairly easy to explain in the present case; it arises from the "externalising" effect that adding a random effect has in a GLM.

Models for dichotomous outcome variables: You fit one model that is a standard GLM and another that is a random effects model that includes a random effect on the subject index:$^\dagger$
$$\begin{matrix}
\text{GLM}   & & & \text{Logit}(p_{ij}) = \gamma_* \quad \quad \\[6pt]
\text{GLMER} & & & \text{Logit}(p_{ij}) = \gamma + \omega_j \\[6pt]
\end{matrix}$$
This leads you to the following estimates for the intercept terms $\gamma_*$ (in red) and $\gamma$ (in blue).

When you fit the initial GLM, the parameter $\gamma_*$ is an estimate of the location of the true probability $p_{ij}$ for the dichotomous outcome, taking into account both the variation over eyes and also the variation across subjects.  Since this is using a lot of information, it gives quite a tight estimate for the parameter, as shown by the relatively narrow confidence interval.  Contrarily, when you add a random effects term across subjects in the latter model, the variation of the outcome across subjects is "externalised" into the random effects term, so now the new parameter $\gamma$ is an estimate of the location of the true probability $p_{ij}$ taking into account only the variation over eyes.  Since this is very little information, it gives a very poor estimate for the parameter, as shown by the extremely wide confidence intervals.
This result is really quite unsurprising --- if you add a random effects term across subject then you are "externalising" the variation across subjects so it no longer affects the intercept parameter.  The particular reason that you get very wide confidence intervals in this case is that, presumably, the eye variable is only weakly associated with the dichotomous outcome variables.  If there is low association between these variables then the former gives little information on the latter, and so the range of estimates of the relevant coefficient parameter is large.  (It is also useful to note that the relationship is mediated through the logit function, so it is not the linear association that is germane here.)  If you look "under the hood" at the likelihood functions for each model, you will see that the intercept parameter in the second model is relatively insensitive to changes across subjects (in terms of derivatives, etc.) and this manifests in major differences in the estimated standard error of the intercept parameters in the two models.
As you can see from the above, the issue here is that you are using two very different models to estimate "the same" underlying parameter.  One model incorporates variation across subjects into the estimator, and therefore estimates relatively accurately.  The other model intentionally excludes this information (by externalising it into random effects terms) and therefore gives an estimate using much less information.  It is unsurprising that the results of the two exercises are very different.  Although they are estimating "the same" parameter, they are effectively using two very different sets of information.

Models for continuous outcome variables: In these cases you can see that the same phenomenon is occurring to some extent ---i.e., the confidence intervals under the random effects model are wider than under the corresponding models without those random effects.  The size of the effect is substantially smaller in this case, and as you can see, the difference in the width of the confidence intervals is much smaller.  Presumably this occurs because the eye variable is giving more information about the continuous outcome variables than the dichotomous outcome variables, and so the "remaining information" is larger in the continuous case.  It is also worth bearing in mind that this model posits a linear association between the variables, so the coefficient is more sensitive to outcome in the extremes of the range, and this may lead to the eye variable being more "informative" in the continuous case.

$^\dagger$ Note that I have used $\gamma_*$ instead of $\gamma$ for the GLM, in order to differentiate the parameters of different models.
A: In the model for the continuous outcome $y$,
$$E(y_{ij})=\mu+\alpha_j+\epsilon_{ij}$$
$\alpha_j$ is measured in units of whatever your outcome variable is.
In the model for the binary outcome $p$,
$$Logit(p_{ij})=\gamma+\alpha_j$$
$\alpha_j$ is measured in units of the log odds.
This is clearly a problem!
I think this could be addressed by adding a scaling parameter to the first model,
$$E(y_{ij})=\mu+\beta\alpha_j+\epsilon_{ij}$$
where $\beta$ captures the mapping between the random effects in the binary model, measured in log-odds, and those in the continuous model, measured in units of $y$.
A: You are not supposed to compare models that have different assumptions. The classical GLM assumes independent data which you stated that this assumption is violated! So, you can't trust the output of such model.
The other point about GLMM (glmer) models, you have to come out with the best fit for the models first, for example compare the two models that have different structures of random effects using the
-2 * logLik(fit1) + 2 * logLik(fit2)
then decide which fit is better.
You can also use model diagnostics such as in "DHARMa" package to be more sure about the fit and the assumptions.
Note:
The number of random effects units should be at least 5-6, but you have only two ~(eyes) and this could pose a problem in the CIs, check out : http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#inference-and-confidence-intervals
"Clark and Linzer (2015)... One point of particular relevance to ‘modern’ mixed model estimation (rather than ‘classical’ method-of-moments estimation) is that, for practical purposes, there must be a reasonable number of random-effects levels (e.g. blocks) – more than 5 or 6 at a minimum"
