Estimates of hemodynamic response functions (HRF) are often integral parts of event-related fMRI analyses. Although HRFs vary across individuals and brain regions, few studies have investigated how variations affect the results of statistical analyses using the general linear model (GLM). In this study, we empirically estimated HRFs from primary motor and visual cortices and frontal and supplementary eye fields (SEF) in 20 subjects. We observed more variability across subjects than regions and correlated variation of time-to-peak values across several pairs of regions. Simulations examined the effects of observed variability on statistical results and ways different experimental designs and statistical models can limit these effects. Widely spaced and rapid event-related experimental designs with two sampling rates were tested. Statistical models compared an empirically derived HRF to a canonical HRF and included the first derivative of the HRF in the GLM. Small differences between the estimated and true HRFs did not cause false negatives, but larger differences within an observed range of variation, such as a 2.5-s time-to-onset misestimate, led to false negatives. Although small errors minimally affected detection of activity, time-to-onset misestimates as small as 1 s influenced model parameter estimation and therefore random effects analyses across subjects. Experiment and analysis design methods such as decreasing the sampling rate or including the HRF's temporal derivative in the GLM improved results, but did not eliminate errors caused by HRF misestimates. These results highlight the benefits of determining the best possible HRF estimate and potential negative consequences of assuming HRF consistency across subjects or brain regions.