Glassy liquid crystalline systems are expected to show significant history-dependent effects. Two model glassy systems are the RAN and SSS lattice models. The RAN model is a Lebwohl-Lasher lattice model with locally coupled nematic spins, together with uncorrelated random anisotropy fields at each site, while the SSS model has a finite concentration of impurity spins frozen in random directions. Here Brownian simulation is used to study the effect of different sample histories in the low temperature regime in a three dimensional (d = 3) model intermediate between SSS and RAN, in which a finite concentration p < pc (pc the percolation threshold) of frozen spins interacts with neighboring nematic spins with coupling W . Simulations were performed at temperature T ∼ TNI /2 (TNI the bulk nematic-isotropic transition temperature) for temperature-quenched and field-quenched histories (TQH and FQH respectively), as well as for temperature-annealed histories (AH). The first two of these limits represent extreme histories encountered in typical experimental studies. Using long-time averages for equilibrated systems, we calculate orientational order parameters and two-point correlation functions. Finite size-scaling was used to determine the range of the orientational ordering, as a function of coupling strength W, p and sample history. Sample history plays a significant role; for given concentration p, as disorder strength W is increased, TQH systems sustain quasi-long-range (QLRO) and short-range-order (SRO). The data are also consistent with a long-range order (LRO) phase at very low disorder strength. By contrast, for FQH, only LRO and QLRO occur within the range of parameters investigated. The crossover between regimes depends on history, but in general, the FQH phase is more ordered than than the AH phase, which is more ordered than the TQH phase. We detect also in the QLRO phase a domain-type structural pattern, consistent with ideas introduced by Giamarchi and Doussal [Phys. Rev. B 52, 1242Rev. B 52, (1995] on superconducting flux lattices. In phases in which short-range exponential order occurs, the orientational correlation length in the weak-disorder limit obeys Larkin-Imry-Ma scaling ξ ∼ D −2/(4−d) .