def METHOD_NAME(hyp, cov_struct=None, mcrep=500):
    # Storage for the simulation results
    scales = [[], []]
    # P-values from the score test
    pv = []
    # Monte Carlo loop
    for k in range(mcrep):
        # Generate random "probability points" u  that are uniformly
        # distributed, and correlated within clusters
        z = np.random.normal(size=n)
        u = np.random.normal(size=n // m)
        u = np.kron(u, np.ones(m))
        z = r * z + np.sqrt(1 - r**2) * u
        u = norm.cdf(z)
        # Generate the observed responses
        y = negbinom(u, mu=mu[hyp], scale=scale)
        # Fit the null model
        m0 = sm.GEE(y,
                    x0,
                    groups=grp,
                    cov_struct=cov_struct,
                    family=sm.families.Poisson())
        r0 = m0.fit(scale='X2')
        scales[0].append(r0.scale)
        # Fit the alternative model
        m1 = sm.GEE(y,
                    x,
                    groups=grp,
                    cov_struct=cov_struct,
                    family=sm.families.Poisson())
        r1 = m1.fit(scale='X2')
        scales[1].append(r1.scale)
        # Carry out the score test
        st = m1.compare_score_test(r0)
        pv.append(st["p-value"])
    pv = np.asarray(pv)
    rslt = [np.mean(pv), np.mean(pv < 0.1)]
    return rslt, scales