Let's start by unpacking Eq. (1), which we rewrite here:
$$\begin{equation}
(H_{SCF} - \varepsilon_n)\lvert \Delta \psi_n \rangle = -(\Delta V_{SCF} + \Delta \varepsilon_n) \lvert \psi_n \rangle.
\tag{1}\end{equation}$$
Well, we know (Eq. (28) of Baroni et al.) that the variation in wavefunction $n$ can be written as
$$\begin{equation}
\Delta \psi_n(r) = \sum_{m \neq n} \frac{\langle \psi_m | \Delta V_{SCF} | \psi_n \rangle}{\varepsilon_n - \varepsilon_m} \psi_m(r),
\tag{28}\end{equation}$$
where $n, m$ range over both occupied states, $\lvert \psi \rangle \in occ$, and virtual states, denoted $\lvert \phi \rangle \not\in occ$. (Note that our expressions will lack a factor of 2 compared to Baroni et al. because we are counting up and down spin orbitals separately.)
Now, the response to the electron density is only summed over occupied states. So, taking $z^*$ as the complex conjugate of $z$, Eqs. (23) and (29) of Baroni et al. read
$$\begin{align}
\Delta n(r) &= \sum_{n \in occ} \left[ \Delta \psi_n^*(r) \psi_n(r) + \psi_n^*(r) \Delta \psi_n(r) \right] \tag{23} \\\\
&= \sum_{n \in occ} \sum_{m \neq n} \left[ \psi_m^*(r) \frac{\langle \psi_n | \Delta V_{SCF} | \psi_m \rangle}{\varepsilon_m - \varepsilon_n} \psi_n(r) + \psi_n^*(r) \frac{\langle \psi_m | \Delta V_{SCF} | \psi_n \rangle}{\varepsilon_n - \varepsilon_m} \psi_m(r) \right]. \tag{29}
\end{align}$$
If we let $m$ range over occupied states, we notice that the contributions to $\Delta n$ coming from orbitals $m,n$ and from $n,m$ cancel; so we can take $m$ to range over the virtual states only. This is succinctly described a little later in Baroni et al.:
the response of the system to an external perturbation depends only on the component of the perturbation that couples the occupied-state manifold with the empty-state one.
Let's return to Eq. (1) above. Implicit in passing to Eq. (2) is the idea that $\lvert \psi_n \rangle$ is an occupied state: that means that $\lvert \Delta \psi_n \rangle$ lives in the virtual manifold. That's why we multiply the right side of Eq. (2) onto $P_c$, the unoccupied projector. What we really meant was
$$\begin{equation}
(H_{SCF} - \varepsilon_n)\lvert \Delta \psi_n \rangle = - P_v (\Delta V_{SCF} + \Delta \varepsilon_n) \lvert \psi_n \rangle - P_c (\Delta V_{SCF} + \Delta \varepsilon_n) \lvert \psi_n \rangle,
\end{equation}$$
but the components of $\lvert \Delta \psi_n \rangle \in \text{span}\ {P_v}$ are all zero, so we can just neglect to write it. This gives us the original Eq. (2):
$$\begin{equation}
(H_{SCF} + \alpha P_v - \varepsilon_n) \lvert \Delta \psi_n \rangle = -P_c \Delta V_{SCF} \lvert \psi_n \rangle, \qquad n \in occ.
\tag{2}\end{equation}$$
But the occupied-virtual coupling goes both ways! If we let $n$ in Eq. (29) index the (infinitely many, so we just handwave away the functional analysis issues) virtual states, then $m$ has to index occupied states. So if we assume $\lvert \psi_n \rangle$ is an unoccupied state in Eq. (1), then $\lvert \Delta \psi_n \rangle$ must belong to the occupied manifold, and we can write
$$\begin{equation}
(H_{SCF} - \varepsilon_n) \lvert \Delta \psi_n \rangle = -P_v \Delta V_{SCF} \lvert \psi_n \rangle \qquad n \not\in occ;
\tag{3}\end{equation}$$
(note that $\Delta \epsilon_n = \langle \psi_n | \Delta V_{SCF} | \psi_n \rangle$ goes away upon the projection just like in the question's equation $1 \to 2$), and we can fix singularities in $H_{SCF} - \varepsilon_n$ by adding a small amount of an operator we know acts trivially on $\lvert \Delta \psi_n \rangle$. So at long last
$$\begin{equation}
(H_{SCF} + \alpha P_c - \epsilon_n) \lvert \Delta \psi_n \rangle = -P_v \Delta V_{SCF} \lvert \psi_n \rangle, \qquad n \not\in occ. \tag{4}
\end{equation}$$
