Let's consider a causal IIR filter with the following transfer function:
$$H(z)=\frac{B(z)}{A(z)}=\frac{\displaystyle\sum_{n=0}^{M}b[n]z^{-n}}{1+\displaystyle\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$
Note that $M$ need not be equal to $N$. Let $K$ be the number of samples that must be removed from the impulse response corresponding to the system $(1)$, so the new desired impulse response is
$$\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots\tag{2}$$
Let $H_K(z)$ denote the FIR transfer function corresponding to the first $K$ samples of $h[n]$:
$$H_K(z)=h[0]+h[1]z^{-1}+\ldots +h[K-1]z^{-(K-1)}\tag{3}$$
The transfer function of the new filter is then given by
$$\tilde{H}(z)=z^K\left[\frac{B(z)}{A(z)}-H_K(z)\right]=\frac{z^K\big[B(z)-H_K(z)A(z)\big]}{A(z)}\tag{4}$$
As expected, the denominator polynomial remains unchanged, we just get a new numerator polynomial
$$\tilde{B}(z)=z^K\big[B(z)-H_K(z)A(z)\big]\tag{5}$$
Let's first consider the case $K=1$. In this case, the FIR transfer function $H_1(z)$ defined in $(3)$ is trivially given by
$$H_1(z)=h[0]=b[0]\tag{6}$$
Consequently, the new numerator polynomial becomes
$$\tilde{B}(z)=\big(b[1]-b[0]a[1]\big)+\big(b[2]-b[0]a[2]\big)z^{-1}+\ldots +\\+\big(b[N]-b[0]a[N]\big)z^{-(N-1)}\tag{7}$$
if $M=N$. If $M>N$, there are additional terms $b[k]z^{-(k-1)}$, $k=N+1,\ldots,M$, and if $M<N$, there are additional terms $-b[0]a[k]z^{-(k-1)}$, $k=M+1,\ldots,N$.
The computation of the new numerator coefficients according to $(7)$ is very simple and efficient. For $K>1$ this procedure can be applied iteratively. In Matlab/Octave this is (almost) a one-liner:
for k = 1:K
b = [b(2:M+1);zeros(N-M,1)] - b(1)*[a(2:N+1);zeros(M-N,1)];
M = length(b) - 1;
end
Note that for any $K$, the order of the new numerator polynomial $\tilde{B}(z)$ equals $\max(M-K,N-1)$.
Original Answer (correct, but the presentation and solution given above is probably better; the examples below are still valid and hopefully useful)
A solution can be obtained by recomputing the numerator coefficients as follows:
$$\tilde{b}[n]=(\tilde{h}\star a)[n],\qquad n=0,1,\ldots,L-1\tag{3}$$
where $\star$ denotes convolution, $L=\max(M,N)$, and $\tilde{h}[n]$ is related to the impulse response $h[n]$ of the original IIR filter by
$$\tilde{h}[n]=h[n+1],\qquad n=0,1,\ldots,L-1\tag{4}$$
The method can be summarized as follows:
- Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=1,2,\ldots,L$
- Compute the new numerator coefficients from the impulse response and from the denominator coefficients using Eq. $(4)$
I chose a random filter to show the result of the suggested procedure in the figure below (blue: original impulse repsonse: red: new impulse response shifted to the right by one sample for comparison):

The left column below shows the original numerator coefficients alongside the modified numerator coefficients, aligned in such a way that it becomes obvious that for higher indices the new coefficients equal the original coefficients. (This is the case because in this example $M>N$; c.f. Eq. $(7)$ of the new answer above, including the text below Eq. $(7)$).
-0.67010
0.65673 0.92477
-0.61219 -0.74621
-0.60645 -0.53944
-0.50316 -0.50316
0.82855 0.82855
0.20396 0.20396
1.15680 1.15680
-0.29830 -0.29830
0.10605 0.10605
The denominator coefficients of both filters are
a = [1, 0.4, -0.2, 0.1]
EDIT:
This method can be generalized to remove $K$ initial values of the original impulse response, even if $K>M$:
- Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=K,K+1,\ldots,K+L-1$
- Compute the new numerator coefficients according to Eq. $(4)$ with $\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots,L-1$
The resulting numerator coefficients will generally have several trailing zeros (or values sufficiently close to zero) which can be removed.
The following example shows the result of removing $20$ initial values of the impulse response of the filter with coefficients
[b,a] =
1.000000 1.000000
2.000000 -0.187178
3.000000 -0.075960
-1.000000 0.812250

The numerator coefficients of the new filter are
0.95236
0.90498
0.42252
Of course, these coefficients can no longer be interpreted as (an approximation of) a shifted version of the original numerator coefficients.