Your way is already good.
Another way to do it is using the trapezoid rule:
data <- cbind(q$x, q$y)
by(data, q$group, FUN = function(x) trapz(x[, 1], x[, 2]))
The results are nearly the same:
INDICES: 1
[1] 0.9903457
INDICES: 2
[1] 1.000978
INDICES: 3
[1] 0.9811152
This is because at the bandwidth needed to make the graph of the densities look reasonable (interval
in your code), you are very close to what you would get if you could do the actual integral.