I know this is a pretty old question, but I've been looking for a derivation of the expressions for group delay and phase delay on the internet. Not many such derivations exist on the net so I thought I'd share what I found. Also, note that this answer is more of a mathematical description than an intuitive one. For intuitive descriptions, please refer to the above answers. So, here goes:
Let's consider a signal
$$x(t) = a(t) \cos(\omega_0 t)$$
and pass this through an L.T.I. system with frequency response
$$H(j\omega) = e^{j\phi(\omega)}$$
We have considered the gain of the system to be unity because we are interested in analyzing how the system alters the phase of the input signal, rather than the gain. Now, given that multiplication in time domain corresponds to convolution in frequency domain, the Fourier Transform of the input signal is given by
$$X(j\omega) = {1 \over 2\pi}A(j\omega) \cdot \big( \pi\delta(\omega - \omega_0) + \pi\delta(\omega + \omega_0) \big)$$
which amounts to
$$X(j\omega) = {A(j(\omega-\omega_0))+A(j(\omega+\omega_0)) \over 2}$$
Therefore, the output of the system has a frequency spectrum given by
$$B(j\omega) = {e^{j\phi(\omega)} \over 2} \big( A(j(\omega-\omega_0))+A(j(\omega+\omega_0)) \big)$$
Now, to find the inverse Fourier Transform of the above expression, we need to know the exact analytical form for $\phi(\omega)$. So, to simplify matters, we assume that the frequency content of $a(t)$ includes only those frequencies which are significantly lower than the carrier frequency $\omega_0$. In this scenario, the signal $x(t)$ can be viewed as an amplitude modulated signal, where $a(t)$ represents the envelope of the high frequency cosine signal. In the frequency domain, $B(j\omega)$ now contains two narrow bands of frequencies centered at $\omega_0$ and $-\omega_0$ (refer to the above equation). This means that we can use a first order Taylor series expansion for $\phi(\omega)$.
$$\begin{align}
\phi(\omega) &= \phi(\omega_0) + \frac{d\phi}{d\omega}(\omega_0)(\omega - \omega_0) \\
&= \alpha + \beta\omega \\
\end{align}$$
where
$$\begin{align}
\alpha &= \phi(\omega_0) - \omega_0\frac{d\phi}{d\omega}(\omega_0) \\
\beta &= \frac{d\phi}{d\omega}(\omega_0) \\
\end{align}$$
Plugging this in, we can calculate the Fourier transform of the first half of $B(j\omega)$ as
$$\frac{1}{2\pi}\int_{-\infty}^{\infty} \frac{1}{2} A(j(\omega - \omega_0)) e^{j(\omega t + \alpha + \beta\omega)} d\omega$$
Substituting $\omega - \omega_0$ for $\omega'$, this becomes
$$\frac{1}{2\pi}\int_{-\infty}^{\infty} \frac{1}{2} X(j(\omega')) e^{j((\omega' + \omega_0) (t + \beta) + \alpha)} d\omega'$$
which simplifies to
$$a(t + \beta)\frac{e^{j(\omega_0 t + \omega_0\beta + \alpha)}}{2}$$
Plugging in the expressions for $\alpha$ and $\beta$, this becomes
$$a(t + \beta)\frac{e^{j(\omega_0 t + \phi(\omega_0))}}{2}$$
Similarly the other half of the inverse Fourier Transform of $B(j\omega)$ can be obtained by replacing $\omega_0$ by $-\omega_0$. Noting that for real signals, $\phi(\omega)$ is an odd function, this becomes
$$a(t + \beta)\frac{e^{-j(\omega_0 t + \phi(\omega_0))}}{2}$$
Thus, adding the two together, we get
$$b(t) = x(t + \frac{d\phi}{d\omega}(\omega_0)) \cos(\omega_0(t + \tfrac{\phi(\omega_0)}{\omega_0}))$$
Notice the delays in the envelope $a(t)$ and the carrier cosine signal. Group delay $(\tau_g)$ corresponds to the delay in the envelope while phase delay $(\tau_p)$ corresponds to the delay in the carrier. Thus,
$$\tau_g = -\frac{d\phi}{d\omega}(\omega_0)$$
$$\tau_p = -\frac{\phi(\omega_0)}{\omega_0}$$
