The little I did understand from googling on R is that the model here is that we have a binary outcome $Y_{ij}$ some sort of bounded variable $\theta_{ij}$. Where $j$ is the group within which random intercept is kept constant. So the model is: $$ \begin{align} c_j&\sim N\left(\mu,\,\sigma^2\right) \\ \eta_{ij}&=c_j+w_1\sin\theta_{ij}+w_2\cos\theta_{ij}\\ \pi_{ij}&=logit\left(\eta_{ij}\right)\\ Y_{ij}&\sim Bernoulli(\pi_{ij}) \end{align} $$ Where $\mu$, $\sigma$, $w_{1,2}$ are to be fitted With some work you could write the likelihood for your data under this model and then use something like [Fisher Information](https://en.wikipedia.org/wiki/Fisher_information) to extract the standard errors. It also, so happens that for Bernoulli, the algebraic expressions can be relatively simple Instead, it may be easier to supplement this model with $\theta_{ij}\sim Uniform\left(0,2\pi\right)$ generate some $Y_{ij}$, and then use your current fitting procedure to get estimates of $w_{1,2}$. Their distribution, and distribution of their sum of squares will tell you the distributions of the statistics under null hypothesis.