[R-sig-ME] question regarding multiple subjects in lmer

[Francisco Juretig]

I am trying to understand how to handle random effects that come from multiple subjects in R for nonlinear models. I am able to formulate a simple linear model, integrate it via Gaussian Quadrature and compare it against lmer using just one subject. This works perfectly, since I can replicate the log-likelihood that lmer reports perfectly. But when I have a model that contains two subjects (non-nested) I observe a discrepancy. In particular, for the following example I use 54 points and get a log-likelihood of -420.6831 while lmer reports -419.8. I suspect I am not formulating correctly the likelihood. Of course I know that lmer does not maximize the log-likelihood like this, but understanding why I have this discrepancy helps me understand what I am doing.
This is the code

########reading datapoliteness       = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")politeness=politeness[-39,]politeness = politeness[,-1]politeness = politeness[,-3]politeness = politeness[order(politeness$gender,politeness$scenario),]library(lme4)
summary(lmer(frequency ~ 1+ (1|gender) + (1|scenario) ,REML=F, data=politeness))
#######quadrature points and weightsquadrature_weights=c(9.616569769829950237E-41,1.2653200760169131064E-35,1.2587808691128592114E-31,2.9626275871495619963E-28,2.59233715574559135377E-25,1.07064905864221394747E-22,2.41469199997268789E-20,3.27825196978982733433E-18,2.8709892632049633876E-16,1.70734067334684107569E-14,7.1713449978551800189E-13,2.19439846768108755285E-11,5.0147398704734680389E-10,8.7334448358337263426E-9,1.1786188465819040326E-7,1.24981002841129022132E-6,1.05354444418257434272E-5,7.1294347272484865949E-5,3.90516172663158571989E-4,0.00174352542196549101376,0.0063820594170464546055,0.0192464413705395102212,0.04801030616129549637501,0.099382787160472611675,0.17114652799601409275,0.2456424705042445872716,0.294199084527291925209,0.2941990845272919252086,0.245642470504244587272,0.17114652799601409275,0.09938278716047261167547,0.048010306161295496375,0.0192464413705395102212,0.0063820594170464546055,0.00174352542196549101376,3.9051617266315857199E-4,7.12943472724848659488E-5,1.05354444418257434272E-5,1.24981002841129022132E-6,1.1786188465819040326E-7,8.7334448358337263426E-9,5.0147398704734680389E-10,2.19439846768108755285E-11,7.1713449978551800189E-13,1.7073406733468410757E-14,2.87098926320496338764E-16,3.27825196978982733433E-18,2.41469199997268789E-20,1.07064905864221394747E-22,2.5923371557455913538E-25,2.96262758714956199629E-28,1.2587808691128592114E-31,1.2653200760169131064E-35,9.616569769829950237E-41)
quadrature_points=c(-9.58416554752346979315,-8.93368833209619334348,-8.39452562332720216575,-7.912775297812494749828,-7.467875283483621332888,-7.04923709867837791084,-6.65050917856649298516,-6.2675031293697518237,-5.897270957641121499901,-5.5376356924033013317,-5.18692897959288955054,-4.843833634581987470744,-4.507283811819901067553,-4.17639883047848521271,-3.850437668366106474864,-3.528766681491085761193,-3.210836082528293830418,-2.89616239045714748026,-2.584315051873627521821,-2.274906037541697572368,-1.967581597392067095654,-1.662015602658996821714,-1.357904066244222484384,-1.054960541873898933872,-0.7529121774748870925153,-0.4514962498062427320625,-0.150457042920862972965,0.150457042920862972965,0.451496249806242732062,0.752912177474887092515,1.054960541873898933872,1.357904066244222484384,1.662015602658996821714,1.967581597392067095654,2.274906037541697572368,2.584315051873627521821,2.896162390457147480258,3.210836082528293830418,3.52876668149108576119,3.85043766836610647486,4.176398830478485212713,4.507283811819901067553,4.84383363458198747074,5.186928979592889550535,5.537635692403301331704,5.897270957641121499901,6.2675031293697518237,6.650509178566492985164,7.04923709867837791084,7.467875283483621332888,7.912775297812494749828,8.39452562332720216575,8.93368833209619334348,9.58416554752346979315)
##Gaussian Quadratureqz <- unique(politeness$gender)mu            = 192.53sdx1          = 54.43sdx2          = 13.4sdx3          = 35.06loop          = 0;f_likelihood  = 0;obs_used      = 0;loop          = 0;qpoints       = 54;for (obs in 1:2){ level_1 = qz[obs] xl      = subset(politeness,politeness$gender == level_1); lev2unique = unique(politeness$scenario)
 total_integral= 0;  for (simulation2 in 1:qpoints){  xx1=quadrature_points[simulation2] likelihood=numeric(length(lev2unique)) for (simulation1 in 1:qpoints){  for(j in 1:length(lev2unique)){ x   = subset(xl,xl$scenario == lev2unique[j]); integral  = 1;  xx2=quadrature_points[simulation1] for (i in 1:nrow(x)){ integral = integral * dnorm(x$frequency[i]-mu-(sqrt(2)*sdx1*xx1)-(sqrt(2)*sdx2*xx2),0,sdx3); } likelihood[j] = likelihood[j] + quadrature_weights[simulation1]*integral/sqrt(pi); } } multiplic=1; for (m in 1:length(likelihood)){ multiplic = multiplic * (likelihood[m]); } total_integral  = total_integral + quadrature_weights[simulation2]*multiplic /sqrt(pi); } f_likelihood = f_likelihood + log(total_integral); }print(f_likelihood) 
	[[alternative HTML version deleted]]