Motivation:
A permutation of a set with no fixed points is called a derangement. The number of derangements of $n$ elements is notated as $!n$ or "$n$ subfactorial".
The relation $$!n=(n-1)(!(n-1)+!(n-2))$$ implies that $!n\equiv0\pmod{n-1}$. Likewise, $$!n=n(!(n-1))+(-1)^n$$ implies that $!n\equiv(-1)^n\pmod{n}$. I wanted to see if there was any similar relation involving other factors $n+k$. I wasn't able to find one; instead, I stumbled upon something weirder.
A weird pattern:
When looking at the table of $!n\pmod{n+1}$ (where by $\text{mod}$ we mean the least non-negative residue), there's no apparent pattern. However, if you look close enough, you'll notice that a disproportionate amount of elements are multiples of $2$. After calculating the first $5000$ values, this number seems to converge to around $76\%$, greater than the $50\%$ you'd expect if this sequence was random.
Likewise, a disproportionate amount of elements in the table of $!n\pmod{n+2}$ are multiples of $3$: $56\%$, contrasted to the expected $33\%$.
This trend continues: there's no single $k\leq1000$ for which there are less than the expected $\frac1{k+1}$ from the first $\require{cancel}\cancel{5000}$ $50000$ values of $!n\pmod{n+k}$ that are multiples of $k+1$. Looking at the data, I'd even argue that the trend gets more pronounced the farther you go. Here's a Pastebin with my data [edit: and with the predicted data from Haran's answer].
My question:
Why does this pattern hold? Why is $!n\pmod{n+k}$ a multiple of $k+1$ so much more often than it would be by chance? And where do these percentages come from?
Here is a table of the first few derangement values: $$ \begin{array}{c}\\ n & !n \\ \hline 1 & 0 \\ 2 & 1\\ 3 & 2 \\ 4 & 9 \\ 5 & 44 \\ 6 & 265 \\ 7 & 1854 \\ 8 & 14833 \\ 9 & 133496 \\ 10 & 1334961 \\ \end{array} $$
First, we can observe that all odd inputs of $n$ are giving even outputs while all even inputs of $n$ are giving odd outputs. This can be justified using induction. It clearly works for the first few cases from the above table. Now, by induction hypothesis, $!n+!(n-1) \equiv 1 \pmod{2}$ which shows that $!(n+1) \equiv n \pmod{2}$ which shows that odd inputs give even outputs while even inputs give odd outputs.
For all odd values of $n$, we can see that $!n$ and $n+1$ are even which shows that the residue $!n \pmod{n+1}$ is also even. Thus, we have already acquired $50\%$ of the multiples of $2$. On the other hand, for even values of $n$, it seems like there is a $50\%$ chance of the least residue being odd or even, as expected. This gives us an expected ratio of $50\%+\frac{1}{2}(50\%)=75\%$ of even least residues which agrees with your data.
For modulo $3$, you can observe that $!n \pmod{3}$ is $0,1,2,0,2,1 \pmod{3}$ based on $n \pmod{6}$. As our formula is based on recursion, it is easy to see that similar to the above case, we can evaluate $!n \pmod{3}$ using $n \pmod{6}$. Again, like the above case, whenever $n \equiv 1,4 \pmod{6}$ i.e. whenever $n \equiv 1 \pmod{3}$, we have $!n \equiv 0 \pmod{3}$. Thus, the least residue of $!n \pmod{n+2}$ will be divisible by $3$. This is $1$ out of every $3$ cases. That gives an expected ratio of $\frac{1}{3}(100\%)+\frac{2}{3}(33.33\%)=55.55\%$, agreeing with your data.
We would like to show that in general, $!n \equiv 0 \pmod{m}$ for $n\equiv 1 \pmod{m}$. This holds for base case $n=1$ as $!1=0$. We can use the summation formula for derangements to prove this easily- $$!(qm+1)=(qm+1)!\sum_{i=0}^{qm+1} \frac{(-1)^i}{i!}$$ We only need to consider the last two terms as the rest of the terms are divisible by $qm$. So- $$!(qm+1) \equiv (qm+1)(-1)^{qm}+(-1)^{qm+1} \equiv (-1)^{qm}+(-1)^{qm+1} \equiv 0 \pmod{m}$$
We have proved our hypothesis that if $n \equiv 1 \pmod{m}$ then $!n \equiv 0 \pmod{m}$. When we take $!n \pmod{n+k}$ and check its value $\pmod{k+1}$, if $\gcd(k+1,n+k)=1$, we can always expect equal chance for all values $\pmod{k+1}$. So, if $k+1$ is prime, we would expect the expected ratio to be- $$\frac{1}{k+1}(100\%)+\frac{k}{k+1}\bigg(\frac{1}{k+1} \cdot 100\%\bigg) = \frac{2k+1}{(k+1)^2}(100\%)=\bigg(1-\frac{k^2}{(k+1)^2}\bigg)(100\%)>\frac{100\%}{k+1}$$
which agrees with values such as $k=1,2,4,6$ where $k+1$ is prime.
When $k+1=\prod_{i=1}^tp_i^{a_i}$ is composite, we cannot always say that $\gcd(k+1,n+k)=1$. Thus, we must check for all possible $\gcd$ values. If we have $\gcd(k+1,n+k)=\prod_{i=1}^tp_i^{b_i}$, we have $n \equiv 1 \pmod{\prod_{i=1}^tp_i^{b_i}}$ and $n+k \equiv 0 \pmod{\prod_{i=1}^tp_i^{b_i}}$. This $\gcd$ is seen for $\phi(\prod_{i=1}^tp_i^{a_i-b_i})$ cases. Moreover, we must also have divisibility by $\prod_{i=1}^tp_i^{a_i}$ given divisibility by $\prod_{i=1}^tp_i^{b_i}$, which happens with a chance of $1$ in $\prod_{i=1}^tp_i^{a_i-b_i}$. Thus our expected value is: $$E=\frac{1}{\prod_{i=1}^tp_i^{a_i}}\sum \frac{\phi(\prod_{i=1}^tp_i^{a_i-b_i})}{\prod_{i=1}^tp_i^{a_i-b_i}}=\frac{1}{k+1} \cdot \sum_{d \mid (k+1)} \frac{\phi(d)}{d}=\frac{1}{k+1} \cdot \prod_{i=1}^t \bigg(1+\frac{p_i-1}{p_i} \cdot a_i\bigg)$$
Thus, we conclude that for $k+1=\prod_{i=1}^tp_i^{a_i}$, the expected value (in fraction) is: $$E=\frac{1}{k+1} \cdot \prod_{i=1}^t \bigg(1+\frac{p_i-1}{p_i} \cdot a_i\bigg)$$
This seems to agree with all values of $k$ settling the problem.