I am conducting a meta-analysis of diagnostic test accuracy studies focusing on myocardial perfusion imaging.
I have used first Meta-Disc, but only for descriptive purposes, as it is clear that univariate approaches such as those provided by this package are biased (eg Takwoingi et al). Then, I have found the following bivariate methods, and used several of them:


Bayesian bivariate model using bamdit in R; 
Bayesian bivariate model using meta4diag in R;
Bayesian HSROC using HSROC in R; 
frequentist bivariate model using metamisc in R; 
frequentist bivariate model using metandi in Stata;
frequentist copula mixed model using CopulaREMADA in R;
frequentist hierarchical summary receiver operating characteristic (HSROC) model using metandi in Stata;
frequentist proportional hazard model using mada in R;
frequentist Reitsma model using mada in R;
frequentist Reitsma model using Metatron in R.


Results are similar across many of these methods, albeit obviously not identical. Yet, I would favor the Reitsma model as available in mada as it gives me more comprehensive analytical and graphical results. 

My questions stem from this actual project but are quite more general. 

Is there a method which is best for meta-analysis of diagnostic test accuracy studies? Or are they more or less similar? Is there any other method not listed above which is better still?

Your clue to figuring this out should be that the parameter estimates from the scikit-learn estimation are uniformly smaller in magnitude than the statsmodels counterpart. This might lead you to believe that scikit-learn applies some kind of parameter regularization. You can confirm this by reading the scikit-learn documentation.
There is no way to switch off regularization in scikit-learn, but you can make it ineffective by setting the tuning parameter C to a large number. Here is how that works in your case:
# module imports
from patsy import dmatrices
import pandas as pd
from sklearn.linear_model import LogisticRegression
import statsmodels.discrete.discrete_model as sm

# read in the data & create matrices
df = pd.read_csv("http://www.ats.ucla.edu/stat/data/binary.csv")
y, X = dmatrices('admit ~ gre + gpa + C(rank)', df, return_type = 'dataframe')

# sklearn output
model = LogisticRegression(fit_intercept = False, C = 1e9)
mdl = model.fit(X, y)
model.coef_

# sm
logit = sm.Logit(y, X)
logit.fit().params

UPDATE: As correctly pointed out in the comments below, now you can switch off the relularization in scikit-learn by setting penalty='none' (see the docs).
