Occupancy models involve both the probability a site is occupied and the probability occupancy is detected. The homogeneous occupancy model, where the occupancy and detection probabilities are the same at each site, admits an orthogonal parameter transformation that yields a two-stage process to calculate the maximum likelihood estimates so that it is not necessary to simultaneously estimate the occupancy and detection probabilities. The two-stage approach is examined here for the heterogeneous occupancy model where the occupancy and detection probabilities now depend on covariates that may vary between sites and over time. There is no longer an orthogonal transformation but this approach effectively reduces the parameter space and allows fuller use of the R functionality. This permits use of existing vector generalised linear models methods to fit models for detection and allows the development of an iterative weighted least squares approach to fit models for occupancy. Efficiency is examined in a simulation study and the full maximum likelihood and two-stage approaches are compared on several data sets.