Schistosomiasis is a helminth (worm) infection with approximately 200 million people infected worldwide. There is still controversy on whether differing worm burdens between individuals is a result of differing contact rates or of acquired immunity. In this paper, we present a stochastic modelling framework to address these issues. By using appropriate approximations for the higher moments of the joint distributions, differential equations for the means, variances and co-variances of infection and immunity can be obtained. Numerical solutions of these equations to obtain age profiles of the above properties were compared with Monte Carlo simulations of the stochastic process. Simulations showed that the results depended on whether between host heterogeneity was generated by differing contact rates or differing immune responses. Heterogeneity in contact rates produced a highly aggregated distribution of parasites with a large variance/mean ratio. Heterogeneity in the immune response had very little effect on the overall dynamics. This agrees with the predominant field evidence which would suggest that infection is mainly determined by ecology with a smaller contribution of immunity.