Thank you @nirgrahamuk. With that problem solved I can make the original code work with mf. I checked it with and without mf, and I get the same answer for dff1$p. So any thoughts on how to calculate the Mahalanobis distance between males and females?
dff1 = data.frame(score = c(91, 93, 72, 87, 86, 73, 68, 87, 78, 99, 95, 76, 84, 96, 76, 80, 83, 84, 73, 74),
hours = c(16, 6, 3, 1, 2, 3, 2, 5, 2, 5, 2, 3, 4, 3, 3, 3, 4, 3, 4, 4),
prep = c(3, 4, 0, 3, 4, 0, 1, 2, 1, 2, 3, 3, 3, 2, 2, 2, 3, 3, 2, 2),
grade = c(70, 88, 80, 83, 88, 84, 78, 94, 90, 93, 89, 82, 95, 94, 81, 93, 93, 90, 89, 89),
mf = c("m","f","f","m","m","m","f","m","m","f","f","f","m","f","f","m","m","f","f","f"))
dff1$maha1<-mahalanobis(dff1[,1:4], colMeans(dff1[,1:4]), cov(dff1[,1:4]))
dff1$p <- pchisq(dff1$maha1, df=3.0, lower.tail=FALSE) #Mahalanobis distance is distribted chi-square k-1 degrees of freedom where k is the number of variables.
dff1$p