Circadian rhythms of parasites and their hosts can influence processes such as transmission, pathology and life cycle evolution. For trematode parasites that depend on free-living infectious stages (i.e. cercariae) to move among host species, the timing of parasite release is hypothesized to increase the likelihood of contacting a host. Yet, a persistent challenge in studying such biorhythms involves selection of appropriate analytical techniques. Here, we extend a generalized linear mixed modelling (GLMM) framework to cosinor analyses, thereby allowing flexibility in the statistical distribution of the response variable, incorporation of multiple covariates and inclusion of hierarchical grouping effects. By applying this approach to 93 snails infected with trematode parasites from freshwater pond ecosystems, we detected non-random rhythms in six of eight species, with variation in both the timing of peak cercariae release (between 5:10 and 21:46 h) and its magnitude (between 13 and 386). The use of GLMM yielded more accurate and precise estimates of the cosinor parameters compared with classical least-squares (LS) based on a simulation-based sensitivity analysis. The sensitivity analysis revealed that the amplitude and rhythm-adjusted mean values from the LS models diverged from the true values at some limits. We highlight the importance of novel analytical approaches for evaluating parasite circadian rhythms and investigating their underlying mechanisms.