Let $d$ be a positive number.
There is a two dimensional recurrence relation as follow:
$$R(n,m) = R(n-1,m-1) + R(n,m-d)$$
where
$R(0,m) = 1$ and $R(n,0) = R(n,1) = \cdots = R(n, d-1) = 1$ for all $n,m>0$.


How to analyze the asymptotics of $R(n, kn)$ for fixed $k$?
It is easy to see that
$$R(n, kn) = O\left( c\_{k,d}^{n} \cdot (n+k+d)^{O(1)} \right)$$
Is there a way (or an algorithm) to find $c\_{k,d}$ given $k$ and $d$?


PS: I have calcuated the bivariate generating function of $R(\cdot, \cdot)$:
\begin{align}
f(x,y) 
&= \frac{1 - xy - y^{d} + xy^{d}}{(1 - x)(1 - y)(1 - xy - y^{d})} \\
&= \frac{1}{(1 - x)(1 - y)} + \frac{xy^{d}}{(1 - x)(1 - y)(1 - xy - y^{d})} \\
\end{align}



How to find the asymptotics of a linear two-dimensional recurrence relation
Your generating function $f(x,y)$ is convergent on a polydisk $|x|<\epsilon$, $|y|< \delta$, for some $\epsilon, \delta < 1$. We can reduce $\epsilon$ so that $\epsilon \ll \delta^k$ and the domain of $f(x,y)$ to the product of the polydisk $|x|<\epsilon$ and the annulus $\delta' < |y| < \delta$, with $\delta' \sim \delta$. Then $g(z,y) = f(z/y^k, y)$ is well-defined and analytic on the product domain of $|z|\lesssim \epsilon/\delta^k$ and $\delta' < |y| < \delta$. The generating function $g(z,y)$ has a unique expansion that is a power series in $z$ and a Laurent series in $y$. Since $x^n y^m = (z/y^k)^n y^m = z^n y^{m-kn}$, the coefficients $R(n,kn)$ are the coefficients of $z^n y^0$ in this double expansion. Hence, using the Cauchy integral formula, the generating function for $R(n,kn)$ is now given by the contour integral
$$
 \sum\_{n=0}^\infty R(n,kn) z^n = \frac{1}{2\pi i} \oint\_{|y|\sim \delta} g(z,y) \frac{dy}{y} .
$$
Of course, the integration contour is now allowed to deform, as long as it doesn't hit any singularities of the integrand.


The above integrand has a convenient partial fraction expansion expansion with respect to $z$,
$$
 \frac{g(z,y)}{y} = \frac{y^k}{(1-y-y^d)(y^k-z)}
 - \frac{y^{d+k-1}(1-y^d)}{(1-y)(1-y-y^d) (y^{k-1} (1-y^d)-z)} .
$$
Expanding in $z$ under the integral gives the formula
$$
 R(n,kn) = \frac{1}{2\pi i} \oint \frac{dy}{y^{kn} (1-y-y^d)}
 - \frac{1}{2\pi i} \oint \frac{y^d dy}{(1-y)(1-y-y^d) y^{n(k-1)} (1-y^d)^n} ,
$$
whose leading asymptotics can be estimated by residues or steepest descent.


Using residues on the first term, all the contributions come from the roots of $1-y-y^d=0$. The root with the smallest magnitude $y\_\*$ will give the leading contribution. Some experimentation shows that $y\_\*$ is the unique positive real root. The leading asymptotic term then looks like
$$
 C\_\* y\_\*^{-kn} .
$$
In the special case $d=1$, $y\_\* = 1/2$. In the special case $k=1$, the second term has no poles inside its contour and hence evaluates to zero. So in the sequel we can assume that $k>1$.


Applying [steepest descent](https://en.wikipedia.org/wiki/Method_of_steepest_descent) to the second integral, we find the stationary phase points at the roots $y\_\star$ of $y^d = (k-1)/(k+d-1) + O(1/n)$. Let $y\_\star$ be that root (actually we just need its $n\to \infty$ limit) which minimizes the magnitude of $y\_\star^{k-1} (1-y\_\star^d)$. It is easy too see that $y\_\star = \sqrt[d]{(k-1)/(d+k-1)}$. Then the leading asymptotic term is
$$
 C\_\star n^{-1/2} y\_\star^{-n(k-1)} (1-y\_\star^d)^{-n}
 = C\_\star n^{-1/2} w\_\star^{-kn} ,
$$
where $w\_\star = \left(\frac{k-1}{d+k-1}\right)^{\frac{1-1/k}{d}} \left(\frac{d}{d+k-1}\right)^{1/k} = \frac{s^{s/k}}{(s+1)^{(s+1)/k}}$ with $s=(k-1)/d$, as pointed out in the comments.


The conclusion is that
$$
 R(n,kn) = O\left((\min(y\_\*,w\_\star)^{-k})^n\right) .
$$


Experimentation suggests that $y\_\* < w\_\star$ in all cases, except for $(k,d)=(2,1)$, when $y\_\* = w\_\star = 1/2$. Otherwise the $y\_\*$ contribution always seems to dominate over the $w\_\star$ contribution. At least that's how it seemed for not too large values of $d$.


Doing some rudimentary asymptotic calculations for large $d$, it seems that $y\_\* \sim 1- \frac{\log d}{d} + o(\log(d)/d)$. A better approximation is
$$
 y\_\* \sim 1 - \frac{W(d)}{d} + \frac{W(d)^3}{2(W(d)+1) d^2} + O(W(d)^3/d^3) ,
$$
where $W(d)$ is the [Lambert W function](https://en.wikipedia.org/wiki/Lambert_W_function). On the other hand, $w\_\star$ has a minimum as a function of $k$ (treating it as a continuous variable) around $k\_\star \sim W(d) + 1 + O(W(d)^2/d)$, the minimum reaches roughly
$$
 w\_\star \sim 1 - \frac{W(d)}{d-1} \lesssim y\_\* .
$$
One should push the asymptotics of $k\_\star$ one more order to get a better estimate of the difference, but experiments do show that the minimum of $w\_\star$ does dip below $y\_\*$, for instance this happens near $(k,d) = (5, 205)$ or $(6,700)$. So it looks like the $w\_\star$ asymptotic contribution will dominate in small ranges of $k \sim k\_\star(d)$ if it happens to fall near an integer. To get the size of that window, one could go to a quadratic Taylor approximation of $w\_\star$.


With a bit of extra work, one could also extract the coefficients $C\_\*$ and $C\_\star$ from the integral formula.


