††2010 _Mathematics Subject Classification_ : 39B32, 34M05, 30D35.†† _Key
words and phrases_ : Fermat type equation, delay-differential equations, shift
equation, entire and meromorphic solutions, finite order, Nevanlinna theory.

# On the transcendental solutions of Fermat type delay-differential and
$c$-shift equations with some analogous results

Abhijit Banerjee and Tania Biswas  Department of Mathematics, University of
Kalyani, West Bengal 741235, India<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

In this paper, we mainly investigate on the finite order transcendental entire
solutions of two Fermat types delay-differential and one Fermat type $c$-shift
equations, as these types were not considered earlier. Our results improve
those of [13] in some sense. In addition, we also extend some recent results
obtained in [18]. A handful number of examples have been provided by us to
justify our certain assertion as and when required.

Typeset by AmS-LaTeX

##  1\. Introduction and some basic definitions

At the outset, we assume that the readers are familiar with the basic terms
and notations of Nevanlinna’s value distribution theory of meromorphic
functions in the complex plane $\mathbb{C}$. So for such a meromorphic
function $f$, terms like $T(r,f)$, $N(r,f)$, $m(r,f)$ etc., we refer to [6,
8]. The notation $S(r,f)$, is defined to be any quantity satisfying
$S(r,f)=o(T(r,f))$ as $r\rightarrow\infty$, possibly outside a set $E$ of $r$
of finite logarithmic measure. The order of $f$ is defined by

$\displaystyle\rho(f)=\limsup\limits_{r\longrightarrow\infty}\frac{\log
T(r,f)}{\log r}.$

Moreover, the shift and difference operator of a function are represented by
$f(z+c)$ and $\Delta(f)=f(z+c)-f(z),\;c\in\mathbb{C}\backslash\\{0\\}$,
respectively.

In this paper, by linear $c$-shift operator in $f(z)$, we mean

$\displaystyle L_{c}(z,f)=\displaystyle\sum_{j=0}^{\tau}a_{j}f(z+jc),$ (1.1)

where $\tau\geq 1$, $a_{\tau}\neq 0$, $a_{j}$’s are any constants. On the
other hand, for delay-differential operator in $f$, we mean, a finite sum of
products of $f$, shifts of $f$, derivatives of $f$ and derivative of their
shifts $f(z+jc),(c\in\mathbb{C})$, with constants coefficients.

We organize our paper as follows: In Section 2, we investigate on the entire
solutions of Fermat type delay-differential and $c$-shift equations and extend
some previous results. The non-existence conditions of meromorphic solutions
of certain non-linear $c$-shift equations will be considered in Section 3,
which extend [18].

##  2\. Fermat type Delay-Differential and $c$-Shift Equations

Initially, Fermat type equations were investigated by Gross [3, 4], Montel
[16]. Yang [21] investigated the Fermat type equation and obtained the
following result:

######  Theorem A.

[21] Let $m$, $n$ be positive integers satisfying $\frac{1}{m}+\frac{1}{n}<1$.
Then there are no non-constant entire solutions $f(z)$ and $g(z)$ that satisfy

$\displaystyle a(z)f^{n}(z)+b(z)g^{m}(z)=1,$ (2.1)

where $a(z)$, $b(z)$ are small functions of $f(z)$.

From Theorem A it is clear that either $m\geq 2$, $n>2$ or $m>2$, $n\geq 2$.
So, it is natural that the case $m=n=2$ can be treated when $f(z)$ and $g(z)$
have some special relationship in (2.1). This was the starting point of a new
era about the solution of Fermat type equations. As a result, successively
several papers were published (see [1, 9, 10, 11, 12, 15, 19, 20, 22, 24]).

In 2007, Tang-Liao [19] investigated on the transcendental meromorphic
solutions of the following non-linear differential equations

$\displaystyle f(z)^{2}+P(z)^{2}(f^{(k)}(z))^{2}=Q(z),$ (2.2)

where $P(z)$, $Q(z)$ are non-zero polynomials.

In 2013, Liu-Yang [14] considered the existence of solutions of the analogous
difference equations of (2.2) namely

$\displaystyle f(z)^{2}+P(z)^{2}(f(z+c)-f(z))^{2}=Q(z),$ (2.3)

In this paper, we wish to investigate on the existence of solutions of certain
Fermat type delay-differential equation as follows:

$\displaystyle f^{2}(z)+R^{2}(z)(f^{(k)}(z+c)-f^{(k)}(z))^{2}=Q(z),$ (2.4)

where $R(z)$, $Q(z)$ are and non-zero polynomials.

Liu-Yang [14] proved that (2.3) has no finite order transcendental entire
solution, that’s why in (2.4), it will be natural to investigate the case for
$k\geq 1$.

######  Theorem 2.1.

If the non-linear delay-differential equation (2.4) has a transcendental
entire solution of finite order, then $f(z)$ takes the form

$\displaystyle f(z)=\frac{Q_{1}(z)e^{az+b}+Q_{2}(z)e^{-(az+b)}}{2},$

$a,b\in\mathbb{C}$ such that $Q_{1}(z)Q_{2}(z)=Q(z)$. Moreover, one of the
following conclusions hold:

  * (i)

If $e^{ac}\neq 1$; then $k$ must be odd, $Q(z)$ and $R(z)$ reduce to constants
satisfying the equation $2iRa^{k}+1=0$.

  * (ii)

If $e^{ac}=1$; then $\deg{R(z)}=1$, none of $Q_{1}(z)$, $Q_{2}(z)$ be
constants. Also,

$\displaystyle
R(z)=\frac{Q_{1}(z)}{P(Q_{1})}=\frac{Q_{2}(z)}{(-1)^{k-l+1}P(Q_{2})},$

such that $P(x)=i\sum\limits_{l=0}^{k}{k\choose
l}a^{k-l}[x^{(l)}(z+c)-x^{(l)}(z)].$

The following examples clarify that both the cases of Theorem 2.1 actually
hold.

######  Example 2.1.

The function $f(z)=\displaystyle\frac{2e^{3z+2}+3e^{-(3z+2)}}{2}$ satisfies
the Fermat equation
$f^{2}(z)-\frac{1}{36}(f^{\prime}(z+c)-f^{\prime}(z))^{2}=6$, where $c=\pi i$,
$R=-\frac{1}{6i}$.

######  Example 2.2.

The function $f(z)=\displaystyle\frac{\alpha ze^{az+b}+\beta ze^{-(az+b)}}{2}$
satisfies the Fermat equation
$f^{2}(z)-\frac{z^{2}}{a^{2}c^{2}}(f^{\prime}(z+c)-f^{\prime}(z))^{2}=\alpha\beta
z^{2}$ such that $e^{ac}=1$, $Q_{1}(z)=\alpha z$, $Q_{2}(z)=\beta z$, where
$\alpha,\beta\in\mathbb{C}\backslash\\{0\\}$. Here,
$R=\displaystyle\frac{\alpha z}{ia(\alpha.(z+c)-\alpha.z)}=\frac{\beta
z}{ia(\beta.(z+c)-\beta.z)}$, i.e., $R=\displaystyle\frac{z}{iac}$.

In 2015, Liu-Dong [13] investigated on

$\displaystyle D^{2}f^{2}(z)+(Af(z+c)+Bf(z))^{2}=1$ (2.5)

and proved that if there exist finite order transcendental entire solutions of
(2.5), then $A^{2}=B^{2}+D^{2}$. In that paper, they also discussed about

$\displaystyle f^{2}(z)+(Af^{(m)}(z)+Bf^{(n)}(z))^{2}=1$ (2.6)

that (2.6) admits transcendental entire solution when $m+n$ is an even and
$m$, $n$ are odds.

In this paper, we wish to investigate on the following Fermat type equation as
this type of equation was not dealt earlier:

$\displaystyle f^{2}(z)+R^{2}(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))^{2}=1,$ (2.7)

where $m,n\in\mathbb{N}$, $R(z)$ be a non-zero polynomial, $A$ and $B$ are
non-zero constants and prove the following theorem.

######  Theorem 2.2.

If the non-linear delay-differential equation (2.7) has a transcendental
entire solution of finite order, then $R(z)$ reduces to constant, namely, $R$
and $f(z)$ takes the form

$\displaystyle f(z)=\frac{e^{az+b}+e^{-(az+b)}}{2}$

such that when

  * (I)

$m$, $n$ are even, then $a^{m-n}\neq\pm\frac{B}{A}$,
$R^{2}=\displaystyle\frac{1}{a^{2m}A^{2}-a^{2n}B^{2}}$,
$e^{ac}=\displaystyle\frac{-a^{n}B\pm\sqrt{(a^{n}B)^{2}-(a^{m}A)^{2}}}{a^{m}A}$.
Also, $e^{ac}\not\in\left\\{\pm
1,-\left(\displaystyle\frac{a^{m}A}{a^{n}B}\right)^{\pm 1}\right\\}$;

  * (II)

$m$, $n$ are odd, then $e^{ac}=\pm 1$, $a^{m-n}\neq\mp\frac{B}{A}$ and
$R=\displaystyle\frac{-i}{a^{n}B\pm a^{m}A}$;

  * (III)

$m$ is even, $n$ is odd, then $e^{ac}=\pm i$, $a^{m-n}\neq\pm i\frac{B}{A}$
and $R=\displaystyle\frac{-i}{a^{n}B\pm ia^{m}A}$;

  * (IV)

$n$ is even, $m$ is odd, then $a^{m-n}\neq\pm i\frac{B}{A}$,
$R^{2}=\displaystyle-\frac{1}{a^{2n}B^{2}+a^{2m}A^{2}}$,
$e^{ac}=\displaystyle\frac{-a^{n}B\pm\sqrt{(a^{n}B)^{2}+(a^{m}A)^{2}}}{a^{m}A}$.
Also, $e^{ac}\not\in\left\\{\pm
1,\displaystyle\frac{a^{m}A}{a^{n}B},\displaystyle-\frac{a^{n}B}{a^{m}A}\right\\}$.

######  Remark 2.1.

Adopting the same procedure of Case (I) in Theorem 2.2, for $m=n=0$, if we
choose $f(z)=Dg(z)$ and $R=1/D$ in (2.7), then we have
$R^{2}=\displaystyle\frac{1}{A^{2}-B^{2}}\implies A^{2}=B^{2}+D^{2}$, i.e., we
get [13, Theorem 1.13]. So our theorem is a significant extension of the same.

Following example shows that in Theorem 2.1, each of the cases (I)-(IV)
actually occurs.

######  Example 2.3.

(I) Let $m$, $n$ both be even. Consider the function
$f(z)=\frac{e^{2z+3}+e^{-(2z+3)}}{2}$. It is easy to see that $f$ satisfies
the Fermat equation
$f^{2}(z)-\frac{1}{2^{7}}(f^{\prime\prime}(z+c)+3f^{\prime\prime}(z))^{2}=1$,
such that $e^{2c}=-3+2\sqrt{2}$.
(II) Let $m$, $n$ be both odd. Choose $f(z)=\frac{e^{3z+4}+e^{-(3z+4)}}{2}$,
which satisfies the Fermat equation
$f^{2}(z)-\frac{1}{{12}^{2}}(5f^{\prime}(z+c)+f^{\prime\prime\prime}(z))^{2}=1$,
such that $e^{3c}=-1$.
(III) Let $m$ be even and $n$ be odd. The function
$f(z)=\frac{e^{z+2}+e^{-(z+2)}}{2}$ satisfies the Fermat equation
$f^{2}(z)-\frac{1}{8+6i}(f^{\prime\prime}(z+c)+3f^{\prime\prime\prime}(z))^{2}=1$,
such that $e^{c}=i$.
(IV) Let $m$ be odd, $n$ be even. Consider the function
$f(z)=\frac{e^{2z+1}+e^{-(2z+1)}}{2}$ . The function $f$ satisfies the Fermat
equation
$f^{2}(z)-\frac{1}{100}(3f^{\prime}(z+c)+2f^{\prime\prime}(z))^{2}=1$, such
that $e^{2c}=\frac{1}{3}$.

In [14], Liu-Yang obtained the following result:

######  Theorem B.

There is no finite order transcendental entire solution of (2.3), where P(z),
Q(z) are non-zero polynomials.

In this paper, we partially extend the above result in the following manner.

######  Theorem 2.3.

The non-linear $c$-shift equation

$\displaystyle f^{2}(z)+L_{c}^{2}(z,f)=1$ (2.8)

has finite order transcendental entire solution of the form

$f(z)=\frac{e^{az+b}+e^{-(az+b)}}{2},$

satisfying following two equations:

$\displaystyle\left\\{\begin{array}[]{clcr}a_{0}+a_{1}e^{ac}+a_{2}e^{2ac}+\ldots+a_{\tau}e^{\tau
ac}&=&-i,\\\ a_{0}+a_{1}e^{-ac}+a_{2}e^{-2ac}+\ldots+a_{\tau}e^{-\tau
ac}&=&i.\end{array}\right.$ (2.11)

Here $e^{ac}\neq\pm 1$. Also if $\tau=1$, $a_{0}\neq\pm a_{1}$ is required.

######  Remark 2.2.

In Theorem 2.3, if we take $\tau=1$ and $a_{1}=-a_{0}=1$, from (2.11), it is
clear that when $L_{c}(z,f)\equiv\Delta_{c}f(z)$, there exists no finite order
transcendental entire solution, which includes special case of Theorem B.

######  Example 2.4.

Consider the function $f(z)=\sin(\frac{\pi z}{2c})$. Here $a=\frac{i\pi}{2c}$,
$e^{b}=\frac{1}{i}$. Then $f(z)$ be a solution of the equation
$f^{2}(z)+(L_{c}(z,f))^{2}=1$, provided that $\tau$ is an odd integer say
$\tau=2m+1,\;m\geq 1$ and the coefficients of $L_{c}(z,f)$ satisfy the
following simultaneous equations:

$\displaystyle\left\\{\begin{array}[]{clcr}a_{0}-a_{2}+a_{4}-a_{6}+\ldots+(-1)^{m}a_{2m}&=&0,\\\
a_{1}-a_{3}+a_{5}-a_{7}+\ldots+(-1)^{m}a_{2m+1}&=&-1;\end{array}\right.$

and when $\tau$ is an even integer say $\tau=2m,\;m\geq 1$, the coefficients
of $L_{c}(z,f)$ satisfy the following simultaneous equations:

$\displaystyle\left\\{\begin{array}[]{clcr}a_{0}-a_{2}+a_{4}-a_{6}+\ldots+(-1)^{m}a_{2m}&=&0,\\\
a_{1}-a_{3}+a_{5}-a_{7}+\ldots+(-1)^{m-1}a_{2m-1}&=&-1.\end{array}\right.$

The following lemma plays an important part for the proof in this section:

######  Lemma 2.1.

[23] Suppose $f_{j}(z)$ $(j=1,2,...,n+1)$ and $g_{k}(z)$ $(k=1,2,...,n)$
$(n\geq 1)$ are entire functions satisfying the following conditions:

  * (i)

$\sum_{j=1}^{n}f_{j}(z)e^{g_{j}(z)}\equiv f_{n+1}(z)$,

  * (ii)

The order of $f_{j}(z)$ is less than the order of $e^{g_{k}(z)}$ for $1\leq
j\leq n+1$, $1\leq k\leq n$ and furthermore, the order of $f_{j}(z)$ is less
than the order of $e^{g_{h}(z)-g_{k}(z)}$ for $n\geq 2$ and $1\leq j\leq n+1$,
$1\leq h<k\leq n$.

Then $f_{j}(z)\equiv 0$, $(j=1,2,...,n+1)$.

######  Proof of Theorem 2.1.

Assume that $f(z)$ is a finite order transcendental entire solution of (2.4),
then

$\displaystyle[f(z)+iR(z)(f^{(k)}(z+c)-f^{(k)}(z))][f(z)-iR(z)(f^{(k)}(z+c)-f^{(k)}(z))]=Q(z).$
(2.14)

Thus both of $f(z)+iR(z)(f^{(k)}(z+c)-f^{(k)}(z))$ and
$f(z)-iR(z)(f^{(k)}(z+c)-f^{(k)}(z))$ have finitely many zeros. Combining
(2.14), with Hadamard factorization theorem, we assume that

$f(z)+iR(z)(f^{(k)}(z+c)-f^{(k)}(z))=Q_{1}(z)e^{P(z)}$

and

$f(z)-iR(z)(f^{(k)}(z+c)-f^{(k)}(z))=Q_{2}(z)e^{-P(z)},$

where $P(z)$ is a non-constant polynomial, otherwise $f(z)$ will be a
polynomial and $Q_{1}(z)Q_{2}(z)=Q(z)$, where $Q_{1}(z)$, $Q_{2}(z)$ are non-
zero polynomials. We denote the degrees of $Q_{1}(z)$, $Q_{2}(z)$ and $P(z)$
by $k_{1}$, $k_{2}$ and $k_{3}$, respectively. Thus, we have,

$\displaystyle f(z)=\frac{Q_{1}(z)e^{P(z)}+Q_{2}(z)e^{-P(z)}}{2}$ (2.15)

and

$\displaystyle
f^{(k)}(z+c)-f^{(k)}(z)=\frac{Q_{1}(z)e^{P(z)}-Q_{2}(z)e^{-P(z)}}{2iR(z)}.$
(2.16)

From (2.15), we have

$\displaystyle f^{(k)}(z)=\frac{p_{1}(z)e^{P(z)}+p_{2}(z)e^{-P(z)}}{2},$
(2.17)

where,

$\displaystyle p_{1}(z)$ $\displaystyle=$ $\displaystyle
Q_{1}(z)\left[P^{\prime}(z)^{k}+M_{1,k}\left(P^{\prime},P^{\prime\prime},\ldots,P^{(k)}\right)\right]+Q_{1}^{\prime}(z)M_{2,k-1}\left(P^{\prime},P^{\prime\prime},\ldots,P^{(k-1)}\right)$
(2.18)
$\displaystyle+\ldots+Q_{1}^{(k-1)}(z)M_{k,1}(P^{\prime})+Q_{1}^{(k)}(z),$
$\displaystyle p_{2}(z)$ $\displaystyle=$ $\displaystyle
Q_{2}(z)\left[(-1)^{k}P^{\prime}(z)^{k}+N_{1,k}\left(P^{\prime},P^{\prime\prime},\ldots,P^{(k)}\right)\right]+(-1)^{k-1}Q_{2}^{\prime}(z)N_{2,k-1}\left(P^{\prime},P^{\prime\prime},\right.$
(2.19) $\displaystyle\left.\ldots,P^{(k-1)}\right)+\ldots-
Q_{2}^{(k-1)}(z)N_{k,1}(P^{\prime})+Q_{2}^{(k)}(z),$

where $M_{j,k-j+1}(N_{j,k-j+1})$ $(j=1,2)$ are differential polynomials of
$P^{\prime}$ with degree $k-1$. $M_{j,k-j+1}(N_{j,k-j+1})$ are differential
polynomial of $P^{\prime}$ with degree $k-j+1$ $(j=3,4,\ldots,k)$. It follows
that $p_{1}(z)$, $p_{2}(z)$ are polynomials with degree $k_{1}+k(k_{3}-1)\geq
k_{1}$ and $k_{2}+k(k_{3}-1)\geq k_{2}$, respectively. Using (2.16) and (2.17)
we have

$\displaystyle\frac{Q_{1}(z)e^{P(z)}-Q_{2}(z)e^{-P(z)}}{2iR(z)}$ (2.20)
$\displaystyle=$
$\displaystyle\frac{p_{1}(z+c)e^{P(z+c)}+p_{2}(z+c)e^{-P(z+c)}}{2}-\frac{p_{1}(z)e^{P(z)}+p_{2}(z)e^{-P(z)}}{2}$
$\displaystyle=$
$\displaystyle\frac{(p_{1}(z+c)e^{\Delta_{c}{P(z)}}-p_{1}(z))e^{P(z)}+(p_{2}(z+c)e^{-\Delta_{c}{P(z)}}-p_{2}(z))e^{-P(z)}}{2}.$

Then (2.20) can be written as

$\displaystyle\left[p_{1}(z+c)e^{\Delta_{c}{P(z)}}-p_{1}(z)-\frac{Q_{1}(z)}{iR(z)}\right]e^{P(z)}$
$\displaystyle+\left[p_{2}(z+c)e^{-\Delta_{c}{P(z)}}-p_{2}(z)+\frac{Q_{2}(z)}{iR(z)}\right]e^{-P(z)}=0.$
(2.21)

Applying Lemma 2.1 on (2), we have

$\displaystyle
p_{1}(z+c)e^{\Delta_{c}{P(z)}}-p_{1}(z)-\frac{Q_{1}(z)}{iR(z)}=0$ (2.22)

and

$\displaystyle
p_{2}(z+c)e^{-\Delta_{c}{P(z)}}-p_{2}(z)+\frac{Q_{2}(z)}{iR(z)}=0.$ (2.23)

Now we show that $P(z)$ is a one-degree polynomial. If not, suppose that
$\deg(P(z))\geq 2$. Then applying Lemma 2.1 on (2.22) and (2.23), we have
$p_{1}(z+c)=0$ and $p_{2}(z+c)=0$, a contradiction, which implies that $P(z)$
is a one-degree polynomial, say,

$\displaystyle P(z)=az+b,\;a,b\in\mathbb{C}.$ (2.24)

In view of (2.24), (2) yields

$\displaystyle\left[p_{1}(z+c)e^{ac}-p_{1}(z)-\frac{Q_{1}(z)}{iR(z)}\right]e^{P(z)}+\left[p_{2}(z+c)e^{-ac}-p_{2}(z)+\frac{Q_{2}(z)}{iR(z)}\right]e^{-P(z)}=0.$
(2.25)

Using (2.24), the expressions of $p_{1}(z)$ and $p_{2}(z)$ given by (2.18) and
(2.19) reduec to

$\displaystyle p_{1}(z)=\sum\limits_{l=0}^{k}{k\choose
l}a^{k-l}Q_{1}^{(l)}(z)\;\;and\;\;p_{2}(z)=\sum\limits_{l=0}^{k}{k\choose
l}(-a)^{k-l}Q_{2}^{(l)}(z),$ (2.26)

respectively. Now we have to consider the following two cases:
Case 1: If $e^{ac}\neq 1$. Then considering the degrees of $p_{1}$, $Q_{1}$
and $p_{2}$, $Q_{2}$ of the equation (2.25), one can conclude that $R(z)$ is
constant, say, $R$. Thus, using (2.24) and (2.26), (2.22) and (2.23) become

$\displaystyle iR\sum\limits_{l=0}^{k}{k\choose
l}a^{k-l}\left[e^{ac}Q_{1}^{(l)}(z+c)-Q_{1}^{(l)}(z)\right]=Q_{1}(z)$ (2.27)

and

$\displaystyle iR\sum\limits_{l=0}^{k}{k\choose
l}(-a)^{k-l}\left[Q_{2}^{(l)}(z)-e^{-ac}Q_{2}^{(l)}(z+c)\right]=Q_{2}(z).$
(2.28)

Considering the highest degree on both sides of (2.27) and (2.28), we have

$\displaystyle iRa^{k}(e^{ac}-1)=1\;and\;iR(-a)^{k}(1-e^{-ac})=1,$

which implies $e^{ac}=(-1)^{k}$. Since $e^{ac}\neq 1$, $k$ must be odd and

$\displaystyle a=\frac{(2m+1)\pi i}{c},$ (2.29)

$m$ is any integer. So,

$\displaystyle 2iRa^{k}+1=0.$ (2.30)

Eliminating $e^{\Delta_{c}P(z)}$ from (2.22) and (2.23), we have

$\displaystyle\frac{iRp_{1}(z)+Q_{1}(z)}{iRp_{1}(z+c)}=\frac{iRp_{2}(z+c)}{iRp_{2}(z)-Q_{2}(z)}$
(2.31) $\displaystyle\implies$
$\displaystyle(p_{1}(z)p_{2}(z)-p_{1}(z+c)p_{2}(z+c))R^{2}+(p_{1}(z)Q_{2}(z)-p_{2}(z)Q_{1}(z))iR$
$\displaystyle+Q_{1}(z)Q_{2}(z)=0.$

First suppose $Q(z)$ be constant, then $Q_{1}(z)$ and $Q_{2}(z)$ are
constants, then from (2.26) and (2.31) we have, $(2iRa^{k}+1)Q_{1}Q_{2}=0$,
which is possible from (2.30).

Suppose that $Q(z)$ be one-degree polynomial and let
$Q_{1}(z)=\alpha_{1}z+\alpha_{0}$, $\alpha_{1}\neq 0$ and $Q_{2}(z)$ be
constant. Using $e^{ac}=-1$, (2.29), (2.30), we have from (2.27) that

$\displaystyle
iR[a^{k}\\{-(\alpha_{1}(z+c)+\alpha_{0})-(\alpha_{1}z+\alpha_{0})\\}-2ka^{k-1}\alpha_{1}]=\alpha_{1}z+\alpha_{0}$
$\displaystyle\implies$ $\displaystyle
iR[a^{k}\\{2(\alpha_{1}z+\alpha_{0})+\alpha_{1}c\\}+2ka^{k-1}\alpha_{1}]+(\alpha_{1}z+\alpha_{0})=0$
$\displaystyle\implies$
$\displaystyle(2iRa^{k}+1)(\alpha_{1}z+\alpha_{0})+iRa^{k-1}\alpha_{1}(ac+2k)=0$
$\displaystyle\implies$ $\displaystyle 2k=-(2m+1)\pi i,$

a contradiction. Similarly, considering $Q_{1}(z)$ and $Q_{2}(z)$ respectively
as constant and one-degree polynomial, we can get a contradiction in a similar
way. So, $Q(z)$ cannot be one-degree polynomial.

Next suppose that $Q(z)$ is a polynomial of degree $2$. If $Q_{1}(z)$ and
$Q_{2}(z)$ both are one-degree polynomials, then from previous argument, we
get a contradiction. Now let $Q_{1}(z)=\beta_{2}z^{2}+\beta_{1}z+\beta_{0}$,
$\beta_{2}\neq 0$ and $Q_{2}(z)$ be a constant. Using $e^{ac}=-1$, (2.29),
(2.30), we have from (2.27) that

$\displaystyle
iR[a^{k}\\{-(\beta_{2}(z+c)^{2}+\beta_{1}(z+c)+\beta_{0})-(\beta_{2}z^{2}+\beta_{1}z+\beta_{0})\\}$
$\displaystyle+ka^{k-1}\\{-(2\beta_{2}(z+c)+\beta_{1})-(2\beta_{2}z+\beta_{1})\\}-2k(k-1)a^{k-2}\beta_{2}]$
$\displaystyle=\beta_{2}z^{2}+\beta_{1}z+\beta_{0}$ $\displaystyle\implies$
$\displaystyle(2iRa^{k}+1)(\beta_{2}z^{2}+\beta_{1}z+\beta_{0})+2iRa^{k-1}\beta_{2}(ac+2k)z+P_{0}(z)=0,$

where $P_{0}(z)$ is a polynomial of degree $0$. Now, comparing the coefficient
of $z$ from both side of the above equation, again we have $2k=-(2m+1)\pi i$,
a contradiction. Similarly, considering $Q_{1}(z)$ and $Q_{2}(z)$ respectively
as constant and two-degree polynomial, again we get a contradiction. Thus
$Q(z)$ cannot be a second degree polynomial.

Now suppose that $Q(z)$ is a polynomial of degree $n\geq 3$. Let
$Q_{1}(z)=\gamma_{j}z^{j}+\gamma_{j-1}z^{j-1}+\ldots+\gamma_{0}$,
$\gamma_{j}\neq 0$ and $Q_{2}(z)$ be of degree $n-j$, $3\leq j\leq n$ . Again
using $e^{ac}=-1$, (2.29), (2.30), we have from (2.27) that

$\displaystyle(2iRa^{k}+1)Q_{1}(z)+niRa^{k-1}\gamma_{n}(ac+2k)z^{j-1}+P_{j-2}(z)=0,$

$P_{j-2}(z)$ is a polynomial of degree $j-2$. Here, comparing the coefficient
of $z^{j-1}$ from both sides of the above equation, again we have
$2k=-(2m+1)\pi i$, a contradiction. Similarly, $Q_{2}(z)$ cannot be a
polynomial of degree $n-j$, as in that case also we get a contradiction. So,
$Q(z)$ cannot be a polynomial of degree $n\geq 3$.

Thus $Q(z)$ must be constant, say $Q$. Therefore, we must have the form of the
solution is

$\displaystyle f(z)=\frac{Q_{1}e^{az+b}+Q_{2}e^{-az-b}}{2},$

such that $Q_{1}Q_{2}=Q$.
Case 2: If $e^{ac}=1$, then (2.25) becomes

$\displaystyle\left[p_{1}(z+c)-p_{1}(z)-\frac{Q_{1}(z)}{iR(z)}\right]e^{P(z)}+\left[p_{2}(z+c)-p_{2}(z)+\frac{Q_{2}(z)}{iR(z)}\right]e^{-P(z)}=0.$
(2.32)

Applying Lemma 2.1 and using (2.26) on (2.32), we have

$\displaystyle iR(z)\sum\limits_{l=0}^{k}{k\choose
l}a^{k-l}\left[Q_{1}^{(l)}(z+c)-Q_{1}^{(l)}(z)\right]=Q_{1}(z)$ (2.33)

and

$\displaystyle iR(z)\sum\limits_{l=0}^{k}{k\choose
l}(-a)^{k-l}\left[Q_{2}^{(l)}(z)-Q_{2}^{(l)}(z+c)\right]=Q_{2}(z).$ (2.34)

Note that the highest degree of $\sum\limits_{l=0}^{k}{k\choose
l}a^{k-l}[Q_{1}^{(l)}(z+c)-Q_{1}^{(l)}(z)]$ and
$\sum\limits_{l=0}^{k}{k\choose l}(-a)^{k-l}[Q_{2}^{(l)}(z)-Q_{2}^{(l)}(z+c)]$
are $k_{1}-1$ and $k_{2}-1$ respectively. Also we see that none of $Q_{1}(z)$
and $Q_{2}(z)$ be constants, i.e., $\deg{Q(z)}\geq 2$. Comparing the total
degree of the equations (2.33) and (2.34), we can deduce that $\deg{R(z)}=1$.
Moreover,

$\displaystyle
R(z)=\frac{Q_{1}(z)}{P(Q_{1})}=\frac{Q_{2}(z)}{(-1)^{k-l+1}P(Q_{2})}$

such that

$P(x)=i\sum\limits_{l=0}^{k}{k\choose
l}a^{k-l}\left[x^{(l)}(z+c)-x^{(l)}(z)\right].$

∎

######  Proof of Theorem 2.2.

Assume that $f(z)$ is a finite order transcendental entire solution of (2.7),
then

$\displaystyle[f(z)+iR(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))][f(z)-iR(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))]=1.$
(2.35)

Thus both of $f(z)+iR(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))$ and
$f(z)-iR(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))$ have no zeros. Combining (2.35), with
Hadamard factorization theorem, we assume that

$f(z)+iR(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))=e^{P(z)}$

and

$f(z)-iR(z)(Af^{(m)}(z+c)+Bf^{(n)}(z))=e^{-P(z)},$

where $P(z)$ is a non-constant polynomial, otherwise $f(z)$ will be constant.
Thus we have,

$\displaystyle f(z)=\frac{e^{P(z)}+e^{-P(z)}}{2}$ (2.36)

and

$\displaystyle Af^{(m)}(z+c)+Bf^{(n)}(z)=\frac{e^{P(z)}-e^{-P(z)}}{2iR(z)}.$
(2.37)

Using (2.36) in (2.37), we have

$\displaystyle\frac{e^{P(z)}-e^{-P(z)}}{2iR(z)}$
$\displaystyle=\frac{A}{2}\left[p_{1}(z+c)e^{P(z+c)}+p_{2}(z+c)e^{-P(z+c)}\right]+\frac{B}{2}\left[q_{1}(z)e^{P(z)}+q_{2}(z)e^{-P(z)}\right],$
(2.38)

where

$\displaystyle
p_{1}(z+c)=P^{\prime}(z+c)^{m}+M_{1,m-1}\left(P^{\prime}(z+c),P^{\prime\prime}(z+c),\ldots,P^{(m)}(z+c)\right),$
$\displaystyle
p_{2}(z+c)=(-1)^{m}P^{\prime}(z+c)^{m}+M_{2,m-1}\left(P^{\prime}(z+c),P^{\prime\prime}(z+c),\ldots,P^{(m)}(z+c)\right),$
$\displaystyle
q_{1}(z)=P^{\prime}(z)^{n}+N_{1,n-1}\left(P^{\prime},P^{\prime\prime},\ldots,P^{(n)}\right),$
$\displaystyle
q_{2}(z)=(-1)^{n}P^{\prime}(z)^{n}+N_{2,n-1}\left(P^{\prime},P^{\prime\prime},\ldots,P^{(n)}\right),$

where $M_{i,m-1}$, $N_{i,n-1}$ $(i=1,2)$ are differential polynomials of
$P^{\prime}$ with degree $m-1$, $n-1$. Then (2) can be written as

$\displaystyle\left[iAR(z)p_{1}(z+c)e^{\Delta_{c}{P(z)}}+iBR(z)q_{1}(z)-1\right]e^{P(z)}$
$\displaystyle+\left[iAR(z)p_{2}(z+c)e^{-\Delta_{c}{P(z)}}+iBR(z)q_{2}(z)+1\right]e^{-P(z)}=0.$
(2.39)

Applying Lemma 2.1 on (2), we have

$\displaystyle iAR(z)p_{1}(z+c)e^{\Delta_{c}{P(z)}}+iBR(z)q_{1}(z)-1=0$ (2.40)

and

$\displaystyle iAR(z)p_{2}(z+c)e^{-\Delta_{c}{P(z)}}+iBR(z)q_{2}(z)+1=0.$
(2.41)

Now we show that $P(z)$ is a one-degree polynomial. If not, suppose that
$\deg(P(z))\geq 2$. Applying Lemma 2.1 on (2.40) and (2.41), we have
$p_{1}(z+c)=0$ and $p_{2}(z+c)=0$, a contradiction, which concludes that
$P(z)$ is a one-degree polynomial, say,

$\displaystyle P(z)=az+b.$ (2.42)

Using (2.42), we have

$\displaystyle
p_{1}(z)=a^{m},\;p_{2}(z)=(-a)^{m},\;\;q_{1}(z)=a^{n},\;\;q_{2}(z)=(-a)^{n}.$
(2.43)

Using (2.42) and (2.43) in (2), we have

$\displaystyle\left[iAa^{m}R(z)e^{ac}+iBa^{n}R(z)-1\right]e^{az+b}$
$\displaystyle+\left[iA(-a)^{m}R(z)e^{-ac}+iB(-a)^{n}R(z)+1\right]e^{-az-b}=0.$
(2.44)

Again applying Lemma 2.1 on (2), we have

$\displaystyle
iAa^{m}R(z)e^{ac}+iBa^{n}R(z)-1=0\;and\;iA(-a)^{m}R(z)e^{-ac}+iB(-a)^{n}R(z)+1=0.$
(2.45)

From here, we can conclude that $R(z)$ is constant, say, $R$. So, (2.45)
becomes

$\displaystyle iR(a^{m}e^{ac}A+a^{n}B)=1$ (2.46)

and

$\displaystyle iR((-a)^{m}e^{-ac}A+(-a)^{n}B)=-1.$ (2.47)

Case I: Let $m$, $n$ be even. Then (2.47) becomes

$\displaystyle iR(a^{m}e^{-ac}A+a^{n}B)=-1.$ (2.48)

Now eliminating $e^{ac}$ from (2.46) and (2.48) we get,
$R^{2}(a^{2m}A^{2}-a^{2n}B^{2})=1$, which implies,
$\displaystyle\frac{a^{m}A}{a^{n}B}\neq\pm 1$, i.e.,
$a^{m-n}\neq\pm\frac{B}{A}$.
Again from (2.46) and (2.48) we get,
$e^{ac}=\displaystyle\frac{-a^{n}B\pm\sqrt{(a^{n}B)^{2}-(a^{m}A)^{2}}}{a^{m}A}$.
Also, $e^{ac}\not\in\left\\{\pm
1,\left(\displaystyle\frac{a^{m}A}{a^{n}B}\right)^{\pm 1}\right\\}$.

Case II: Let $m$, $n$ be odd. Then (2.47) becomes

$\displaystyle iR(a^{m}e^{-ac}A+a^{n}B)=1.$ (2.49)

So, in this case, from (2.46) and (2.49) we get, $e^{ac}=\pm 1$. Also,
$a^{m-n}\neq\mp\frac{B}{A}$. Then $R=\displaystyle\frac{-i}{a^{n}B\pm
a^{m}A}$.

Case III: Let $m$ be even and $n$ be odd. Then (2.47) becomes

$\displaystyle iR(a^{m}e^{-ac}A-a^{n}B)=-1.$ (2.50)

So, in this case, from (2.46) and (2.50) we get, $e^{ac}=\pm i$ and
$a^{m-n}\neq\pm i\frac{B}{A}$. Then $R=\displaystyle\frac{-i}{a^{n}B\pm
ia^{m}A}$.

Case IV: Let $m$ be odd and $n$ be even. Then (2.47) becomes

$\displaystyle iR(a^{m}e^{-ac}A-a^{n}B)=1.$ (2.51)

Now eliminating $e^{ac}$ from (2.46) and (2.51) we get,
$R^{2}(a^{2n}B^{2}+a^{2m}A^{2})=-1$, which implies,
$\displaystyle\frac{a^{m}A}{a^{n}B}\neq\pm i$. Again from (2.46) and (2.51) we
get,
$e^{ac}=\displaystyle\frac{-a^{n}B\pm\sqrt{(a^{n}B)^{2}+(a^{m}A)^{2}}}{a^{m}A}$.
Also, $e^{ac}\not\in\left\\{\pm
1,\displaystyle\frac{a^{m}A}{a^{n}B},\displaystyle-\frac{a^{n}B}{a^{m}A}\right\\}$.

∎

######  Proof of Theorem 2.3.

Assume that $f(z)$ is a finite order transcendental entire solution of (2.8),
then

$\displaystyle[f(z)+iL_{c}(z,f)][f(z)-iL_{c}(z,f)]=1.$ (2.52)

Proceeding in the same way as done in the previous theorem, from (2.52), we
have,

$\displaystyle f(z)=\frac{e^{P(z)}+e^{-P(z)}}{2}$ (2.53)

and

$\displaystyle L_{c}(z,f)=\frac{e^{P(z)}-e^{-P(z)}}{2i}.$ (2.54)

From (1.1), (2.53) and (2.54), we have

$\displaystyle\frac{e^{P(z)}-e^{-P(z)}}{2i}=\sum\limits_{j=0}^{\tau}a_{j}\frac{e^{P(z+jc)}+e^{-P(z+jc)}}{2}$
(2.55) $\displaystyle\implies$
$\displaystyle\sum\limits_{j=1}^{\tau}a_{j}\left(e^{P(z+jc)}+e^{-P(z+jc)}\right)=-(a_{0}+i)e^{P(z)}-(a_{0}-i)e^{-P(z)}.$

Then (2.55) can be written as

$\displaystyle\left(a_{0}+i+\sum\limits_{j=1}^{\tau}a_{j}e^{\Delta_{jc}P(z)}\right)e^{P(z)}+\left(a_{0}-i+\sum\limits_{j=1}^{\tau}a_{j}e^{-\Delta_{jc}{P(z)}}\right)e^{-P(z)}=0.$
(2.56)

Applying Lemma 2.1 on (2.56), we have

$\displaystyle a_{0}+i+\sum\limits_{j=1}^{\tau}a_{j}e^{\Delta_{jc}P(z)}=0$
(2.57)

and

$\displaystyle a_{0}-i+\sum\limits_{j=1}^{\tau}a_{j}e^{-\Delta_{jc}{P(z)}}=0.$
(2.58)

Now we show that $P(z)$ is a one-degree polynomial. On the contrary, suppose
that $\deg(P(z))\geq 2$. Applying Lemma 2.1 on (2.57) and (2.58), we have
$a_{j}=0$ for all $1\leq j\leq\tau$, which in view of (2.8) implies that
$f(z)$ is constant, a contradiction. So, $P(z)$ is a one-degree polynomial,
say, $P(z)=az+b$. Then using (2.57) and (2.58), the relation between $a$, $c$
and $a_{j}$, $0\leq j\leq\tau$, can be determined by (2.11). Also from (2.11),
it is clear that $e^{ac}\neq\pm 1$. If $\tau=1$, $a_{0}\neq\pm a_{1}$ is
required. ∎

##  3\. Non-linear $c$-Shift Equations

For the existence of solutions of non-linear $c$-shift equation, in 2011, Qi
[17] obtained the following theorems:

######  Theorem C.

[17] Let $q(z)$, $p(z)$ be polynomials and let $n$, $m$ be distinct positive
integers. Then the equation

$\displaystyle f^{m}(z)+q(z)f(z+c)^{n}=p(z)$ (3.1)

has no transcendental entire solutions of finite order.

In 2015, Qi-Liu-Yang [18] obtained the meromorphic variant of Theorem C and
improved this as follows:

######  Theorem D.

[18] Let $f(z)$ be a transcendental meromorphic function with finite order,
$m$ and $n$ be two positive integers such that $m\geq n+4$, $p(z)$ be a
meromorphic function satisfying
$\overline{N}\left(r,\frac{1}{p(z)}\right)=S(r,f)$ and $q(z)$ be a non-zero
meromorphic function satisfying that $T(r,q(z))=S(r,f)$. Then, $f(z)$ is not a
solution of equation

$\displaystyle f^{m}(z)+q(z)f(z+c)^{n}=p(z).$ (3.2)

######  Theorem E.

[18] Let $f(z)$ be a transcendental entire function with finite order, $m$ and
$n$ be two positive integers such that $m\geq n+2$, $p(z)$ be a meromorphic
function satisfying $\overline{N}\left(r,\frac{1}{p(z)}\right)=S(r,f)$ and
$q(z)$ be a non-zero meromorphic function satisfying that $T(r,q(z))=S(r,f)$.
Then $f(z)$ is not a solution of equation (3.2).

In this paper we extend Theorems D-E at the expense of replacing $f(z+c)$ by
$L_{c}(z,f)$.

######  Theorem 3.1.

Let $f(z)$ be a transcendental meromorphic function with finite order, $m$ and
$n$ be two positive integers such that $m\geq(\tau+1)(n+2)+2$, $p(z)$ be a
meromorphic function satisfying
$\overline{N}\left(r,\frac{1}{p(z)}\right)=S(r,f)$ and $q(z)$ be a non-zero
meromorphic function satisfying that $T(r,q(z))=S(r,f)$. Then, $f(z)$ is not a
solution of the non-linear $c$-shift equation

$\displaystyle f^{m}(z)+q(z)(L_{c}(z,f))^{n}=p(z).$ (3.3)

######  Corollary 3.1.

Let $f(z)$ be a transcendental entire function with finite order, $m$ and $n$
be two positive integers such that $m\geq n+2$, $p(z)$ be a meromorphic
function satisfying $\overline{N}\left(r,\frac{1}{p(z)}\right)=S(r,f)$ and
$q(z)$ be a non-zero meromorphic function satisfying that $T(r,q(z))=S(r,f)$.
Then, $f(z)$ is not a solution of the non-linear $c$-shift equation (3.3).

The next examples show that if the condition $m\geq n+2$ is omitted then the
equation (3.3) can admit a transcendental entire solution.

First considering $n=1$ and $m=2$ we have the following examples.

######  Example 3.1.

For an odd integer $s$, the function $f(z)=e^{\frac{s\pi iz}{c}}+z$ is a
solution of the equation $f^{2}(z)-zL_{c}(z,f)=e^{\frac{2s\pi iz}{c}}$, for
$k\geq 2$, provided that the coefficients of $L_{c}(z,f)$ satisfy the
following simultaneous equations:

$\displaystyle\left\\{\begin{array}[]{clcr}&a_{0}-a_{1}+a_{2}-a_{3}+\ldots+(-1)^{k}a_{k}&=2,\\\
&a_{0}+a_{1}+a_{2}+a_{3}+a_{4}+\ldots+a_{k}&=1,\\\
&a_{1}+2a_{2}+3a_{3}+4a_{4}+\ldots+ka_{k}&=0.\end{array}\right.$

Next considering $m=n=1$ we have the following example.

######  Example 3.2.

The function $f(z)=ze^{\frac{\pi iz}{c}}$ satisfies the equation
$f(z)+\frac{1}{z+1}L_{c}(z,f)=\frac{z(z+2)}{z+1}e^{\frac{\pi iz}{c}}$ where
the coefficients of $L_{c}(z,f)$ is chosen such that they satisfy
simultaneously the equations

$\displaystyle\left\\{\begin{array}[]{clcr}a_{0}-a_{1}+a_{2}-\ldots+(-1)^{k}a_{k}&=&1,\\\
-a_{1}+2a_{2}-3a_{3}+\ldots+k(-1)^{k}a_{k}&=&0.\end{array}\right.$

To proceed further we require the following lemmas:

######  Lemma 3.1.

[2, Lemma 5.1] Let $f(z)$ be a finite order meromorphic function and
$\varepsilon>0$, then
$T(r,f(z+c))=T(r,f(z))+o(r^{\sigma-1+\varepsilon})+O(\log r)$ and
$\sigma(f(z+c))=\sigma(f(z))$. Thus, if $f(z)$ is a transcendental meromorphic
function with finite order, then we know $T(r,f(z+c))=T(r,f)+S(r,f)$.

######  Lemma 3.2.

[5, Theorem 2.1] Let $f(z)$ be a meromorphic function with finite order, and
let $c\in\mathbb{C}$ and $\delta\in(0,1)$. Then
$m\left(r,\frac{f(z+c)}{f(z)}\right)+m\left(r,\frac{f(z)}{f(z+c)}\right)=o\left(\frac{T(r,f)}{r^{\delta}}\right)=S(r,f)$.

######  Lemma 3.3.

[7] Let $f$ be a non-constant meromorphic function of finite order and
$c\in\mathbb{C}$. Then

$\displaystyle N(r,\infty;f(z+c))\leq
N(r,\infty;f(z))+S(r,f),\;\;\;\overline{N}(r,\infty;f(z+c))\leq\overline{N}(r,\infty;f)+S(r,f).$

######  Proof of Theorem 3.1.

Suppose by contradiction that $f(z)$ is a transcendental meromorphic function
with finite order satisfying equation (3.3).

If $T(r,p(z))=S(r,f)$, then applying Lemma 3.1 to equation (3.3), we have

$\displaystyle m.T(r,f)$ $\displaystyle=$ $\displaystyle T(r,f^{m})$
$\displaystyle=$ $\displaystyle T(r,p(z)-q(z)(L_{c}(z,f))^{n})$
$\displaystyle=$ $\displaystyle T(r,L_{c}(z,f)^{n})+S(r,f)$
$\displaystyle\leq$ $\displaystyle(\tau+1)n.T(r,f)+S(r,f),$

which contradicts the assumption that $m\geq(\tau+1)(n+2)+2$.

If $T(r,p(z))\neq S(r,f)$, differentiating equation (3.3), we get

$\displaystyle(f^{m}(z))^{\prime}+(q(z)(L_{c}(z,f))^{n})^{\prime}=p^{\prime}(z).$
(3.6)

Next dividing (3.6) by (3.3) we have

$\displaystyle
p^{\prime}(z)[f^{m}(z)+q(z)(L_{c}(z,f))^{n}]=p(z)[(f^{m}(z))^{\prime}+(q(z)(L_{c}(z,f))^{n})^{\prime}]$
(3.7) $\displaystyle\implies$ $\displaystyle
f^{m}(z)=\displaystyle\frac{\frac{p^{\prime}(z)}{p(z)}q(z)(L_{c}(z,f))^{n}-(q(z)(L_{c}(z,f))^{n})^{\prime}}{\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-\frac{p^{\prime}(z)}{p(z)}}.$

First observe that
$\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-\frac{p^{\prime}(z)}{p(z)}$ cannot
vanish identically. Indeed, if
$\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-\frac{p^{\prime}(z)}{p(z)}\equiv 0$,
then we get $p(z)=\alpha f^{m}(z)$, where $\alpha$ is a non-zero constant.
Substituting the above equality to equation (3.3), we have
$q(z)(L_{c}(z,f))^{n}=(\alpha-1)f^{m}(z)$. From Lemma 3.1 and the above
equation, we immediately see as above that
$mT(r,f)\leq(\tau+1)nT(r,f)+S(r,f)$, which is a contradiction to
$m\geq(\tau+1)(n+2)+2$. From equation (3.7), we know

$\displaystyle mT(r,f)$ $\displaystyle=$ $\displaystyle T(r,f^{m})$ (3.8)
$\displaystyle\leq$ $\displaystyle
m\left(r,q(z)(L_{c}(z,f))^{n}\right)+m\left(r,\frac{p^{\prime}(z)}{p(z)}-\frac{(q(z)(L_{c}(z,f))^{n})^{\prime}}{q(z)(L_{c}(z,f))^{n}}\right)$
$\displaystyle+$ $\displaystyle
N\left(r,\frac{p^{\prime}(z)}{p(z)}\;q(z)(L_{c}(z,f))^{n}-(q(z)(L_{c}(z,f))^{n})^{\prime}\right)$
$\displaystyle+$ $\displaystyle
m\left(r,\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-\frac{p^{\prime}(z)}{p(z)}\right)+N\left(r,\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-{\frac{p^{\prime}(z)}{p(z)}}\right)+S(r,f).$

As Lemma 3.1 together with equation (3.3) implies that

$\displaystyle(m-(\tau+1)n)T(r,f)+S(r,f)\leq
T(r,p(z))\leq(m+(\tau+1)n)T(r,f)+S(r,f),$

we conclude that

$\displaystyle S(r,p(z))=S(r,f).$ (3.9)

Applying Lemmas 3.1, 3.2 and (3.9) to equation (3.8), we obtain that

$\displaystyle mT(r,f)$ $\displaystyle\leq$ $\displaystyle
nm(r,f)+N\left(r,\frac{p^{\prime}(z)}{p(z)}q(z)(L_{c}(z,f))^{n}-(q(z)(L_{c}(z,f))^{n})^{\prime}\right)$
(3.10)
$\displaystyle+N\left(r,\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-{\frac{p^{\prime}(z)}{p(z)}}\right)+S(r,f).$

Let

$\displaystyle
H(z)=\frac{p^{\prime}(z)}{p(z)}q(z)(L_{c}(z,f))^{n}-(q(z)(L_{c}(z,f))^{n})^{\prime}$
(3.11)

and

$\displaystyle
G(z)=\frac{(f^{m}(z))^{\prime}}{f^{m}(z)}-{\frac{p^{\prime}(z)}{p(z)}}.$
(3.12)

First of all, we deal with $N(r,H(z))$. From (3.3) and (3.11), we know the
poles of $H(z)$ are at the zeros of $p(z)$ and at the poles of $f(z)$,
$f(z+jc)$, $(j=1,2,\ldots,\tau)$ and $q(z)$. Poles of $p(z)$ will not
contribute towards the poles of $H(z)$ as from the equation (3.3) we know that
the poles of $p(z)$ should be at the poles of $f(z)$, $f(z+jc)$,
$(j=1,2,\ldots,\tau)$ and $q(z)$. We note that $T(r,q(z))=S(r,f)$.

If $z_{0}$ is a zero of $p(z)$ then by (3.11), $z_{0}$ is at most a simple
pole of $H(z)$. If $z_{0}$ is a pole of $f(z)$ of multiplicity $t$ but not a
pole of $f(z+jc)$, $j=1,\ldots,\tau$, then $z_{0}$ will be a pole of $H(z)$ of
multiplicity at most $tn+1$. Next suppose $z_{1}$ be any pole of $f(z)$ of
multiplicity $t_{0}$ and a pole of at least one $f(z+jc)$,
$j=1,2,\ldots,\tau$, of multiplicity $t_{j}\geq 0$. Then $z_{1}$ may or may
not be a pole of $L_{c}(z,f)$. From the above arguments and our assumption, we
conclude that

$\displaystyle N(r,H)$ $\displaystyle\leq$
$\displaystyle\overline{N}\left(r,\frac{1}{p(z)}\right)+N(r,(L_{c}(z,f))^{n})+\overline{N}(r,L_{c}(z,f))+S(r,f)$
(3.13) $\displaystyle\leq$ $\displaystyle
nN(r,L_{c}(z,f))+(\tau+1)\overline{N}(r,f)+S(r,f).$

Next, we turn our attention towards the poles of $G(z)$. We know from (3.3)
and (3.12) that the poles of $G(z)$ are at the zeros of $p(z)$ and $f(z)$ and
at the poles of $f(z)$, $f(z+jc)$, $j=1,2,\ldots,\tau$. If $z_{0}$ is a zero
of $p(z)$, zero of $f(z)$, or pole of $f(z+jc)$, $j=1,2,\ldots,\tau$, then by
(3.12) we know $z_{0}$ will be at most a simple pole of $G(z)$. If $z_{0}$ is
a pole of $f(z)$ but not a pole of $f(z+jc)$, $j=1,2,\ldots,\tau$, then by the
Laurent expansion of $G(z)$ at $z_{0}$, we obtain that $G(z)$ is analytic at
$z_{0}$. Therefore, from our assumption and the discussions above, we know

$\displaystyle N(r,G)$ $\displaystyle\leq$
$\displaystyle\overline{N}\left(r,\frac{1}{p(z)}\right)+\overline{N}(r,L_{c}(z,f))+\overline{N}\left(r,\frac{1}{f}\right)+S(r,f)$
(3.14) $\displaystyle\leq$
$\displaystyle\overline{N}(r,L_{c}(z,f))+\overline{N}\left(r,\frac{1}{f}\right)+S(r,f).$

Using Lemma 3.3, from equations (3.10), (3.13) and (3.14) we have

$\displaystyle mT(r,f)$ $\displaystyle\leq$ $\displaystyle
nm(r,f)+nN(r,L_{c}(z,f))+(\tau+1)\overline{N}(r,f)+\overline{N}(r,L_{c}(z,f))+\overline{N}\left(r,\frac{1}{f}\right)+S(r,f)$
$\displaystyle\leq$ $\displaystyle
nm(r,f)+n(\tau+1)N(r,f)+(\tau+1)\overline{N}(r,f)+(\tau+1)\overline{N}(r,f)+\overline{N}\left(r,\frac{1}{f}\right)+S(r,f)$
$\displaystyle\leq$ $\displaystyle\\{(\tau+1)(n+2)+1\\}T(r,f)+S(r,f),$

which contradicts the assumption that $m\geq(\tau+1)(n+2)+2$. This completes
the proof of the theorem. ∎

Acknowledgement

The authors wish to thank the referee for his/her valuable suggestions towards
the improvement of the paper. The second author is thankful to University
Grant Commission (UGC), Govt. of India for financial suport under UGC-Ref.
No.: 1174/(CSIR-UGC NET DEC. 2017) dated 21/01/2019.

## References

  * [1] M. F. Chen and Z. S. Gao, Entire solutions of differential-difference equation and Fermat type $q$-difference-differential equations, Commun. Korean Math. Soc., 30(4)(2015), 447-456.
  * [2] Y. M. Chiang and S. J. Feng, On the Nevanlinna characteristic of $f(z+\eta)$ and difference equations in the Complex plane, Ramanujan J., 16(1)(2006), 105-129.
  * [3] F. Gross, On the equation $f^{n}(z)+g^{n}(z)=1$, Bull. Amer. Math. Soc., 72(1966), 86-88.
  * [4] F. Gross, On the equation $f^{n}(z)+g^{n}(z)=h^{n}(z)$, Amer. Math. Monthly, 73(1966), 1093–1096.
  * [5] R. G. Halburd and R. J. Korhonen, Nevanlinna theory for the difference operator, Ann. Acad. Sci. Fenn. Math., 31(2006), 463-478.
  * [6] W. K. Hayman, Meromorphic Functions, Clarendon Press, Oxford, 1964.
  * [7] J. Heittokangas, R. Korhonen, I. Laine, J. Rieppo and J. L. Zhang, Value sharing results for shifts of meromorphic function and sufficient conditions for periodicity, J. Math. Anal. Appl., 355(2009), 352-363.
  * [8] I. Laine, Nevanlinna Theory and Complex Differential Equations, Walter de Gruyter, Berlin-New York, 1993.
  * [9] B. Q. Li, On certain non-linear differential equations in complex domains, Arch. Math., 91(2008), 344-353.
  * [10] K. Liu, Meromorphic functions sharing a set with applications to difference equations, J. Math. Anal. Appl., 359(2009), 384-393.
  * [11] K. Liu and T. B. Cao, Entire solutions of Fermat type $q$-difference-differential equations, Electron. J. Diff. Equ., 2013(59)(2013), 1-10.
  * [12] K. Liu, T. B. Cao and H. Z. Cao, Entire solutions of Fermat type differential-difference equations, Arch. Math., 99(2012), 147-155.
  * [13] K. Liu and X. J. Dong, Fermat type differential and difference equations, Electron. J. Diff. Equ., 2015(159)(2015), pp. 1-10.
  * [14] K. Liu and L. Yang, On entire solutions of some differential-difference equations, Cmput. Methods. Funct. Theory, 13(2013), 433-447.
  * [15] K. Liu and L. Z. Yang, A note on meromorphic solutions of Fermat types equations, An. Stiint. Univ. Al. I. Cuza Lasi Mat. (N. S.), 1(2016), 317–325.
  * [16] P. Montel, Lecons sur les familles de nomales fonctions analytiques et leurs applications, Gauthier-Viuars Paris, 1927, 135-136.
  * [17] X. Qi, Value distribution and uniqueness of difference polynomials and entire solutions of difference equations, Ann. Polon. Math., (2011), 129-142.
  * [18] X. Qi, Y. Liu and L. Yang, On meromorphic solutions of certain type of difference equations, Bull. Irn. Math. Soc., 41(1)(2015), 281-289.
  * [19] J. F. Tang and L. W. Liao, The transcendental meromorphic solutions of a certain type of non-linear differential equations, J. Math. Anal. Appl., 334(2007), 517-527.
  * [20] H. Wang, H. Y. Xu and J. Tu, The existence and forms of solutions for some Fermat-type differential-difference equations, AIMS Mathematics, 5(1)(2019), 685-700.
  * [21] C. C. Yang, A generalization of a theorem of P. Montel on entire functions, Proc. Amer. Math. Sci., 26(1970), 332-334.
  * [22] C. C. Yang and P. Li, On the transcendental solutions of a certain type of non-linear differential equations, Arch. Math., 82(2004), 442–448.
  * [23] C. C. Yang and H. X. Yi, Uniqueness Theory of Meromorphic Functions, Science Press and Kluwer Acad. Publ., Beijing, 2003.
  * [24] X. Zhang and L. W. Liao, On a certain type of non-linear differential equations admitting transcendental meromorphic solutions, Science China, 56(10)(2013), 2025-2034.


# Guidance on Individualized Treatment Rule Estimation in High Dimensions

Philippe Boileau
Graduate Group in Biostatistics and
Center for Computational Biology, and
University of California, Berkeley
<EMAIL_ADDRESS>
&Ning Leng
Genentech Inc.
<EMAIL_ADDRESS>
&Sandrine Dudoit
Department of Statistics,
Division of Biostatistics, and
Center for Computational Biology,
University of California, Berkeley
<EMAIL_ADDRESS>

###### Abstract

Individualized treatment rules, cornerstones of precision medicine, inform
patient treatment decisions with the goal of optimizing patient outcomes.
These rules are generally unknown functions of patients’ pre-treatment
covariates, meaning they must be estimated from clinical or observational
study data. Myriad methods have been developed to learn these rules, and these
procedures are demonstrably successful in traditional asymptotic settings with
moderate number of covariates. The finite-sample performance of these methods
in high-dimensional covariate settings, which are increasingly the norm in
modern clinical trials, has not been well characterized, however. We perform a
comprehensive comparison of state-of-the-art individualized treatment rule
estimators, assessing performance on the basis of the estimators’ accuracy,
interpretability, and computational efficacy. Sixteen data-generating
processes with continuous outcomes and binary treatment assignments are
considered, reflecting a diversity of randomized and observational studies. We
summarize our findings and provide succinct advice to practitioners needing to
estimate individualized treatment rules in high dimensions. All code is made
publicly available, facilitating modifications and extensions to our
simulation study. A novel pre-treatment covariate filtering procedure is also
proposed and is shown to improve estimators’ accuracy and interpretability.

_Keywords_ heterogeneous treatment effects $\cdot$ observational studies
$\cdot$ precision medicine $\cdot$ randomized control trials

##  1 Introduction

Broadly, clinical trials are performed to evaluate the efficacy and safety of
novel therapies relative to the standard of care. Efficacy is generally
assessed by performing inference on a pre-defined marginal parameter measuring
the treatment effect, like the average treatment effect. The rise of
personalized medicine is however partially at odds with this approach,
especially when the patient population is large and diverse. Seemingly similar
patients’ response to a given treatment can vary widely, which is often
attributable to the existence of unknown patient subpopulations. Identifying
these patient subgroups is therefore required to maximize patient outcomes and
optimize trial success rates.

These subpopulations are generally defined as an unknown function of patient
characteristics, like age, sex-at-birth, or genetic mutations. Certain pre-
treatment covariates are therefore said to modify the treatment effect. These
variables are referred to as treatment effect modifiers (TEMs). TEMs in turn
may be used to define individualized treatment rules (ITRs) that assign
patients to treatment regimens based on their covariates.

Much effort has therefore been dedicated to developing statistical methods for
learning ITRs, and in particular optimal ITRs (see Robins, 2004; Qian and
Murphy, 2011; Zhang et al., 2012a, b; Tian et al., 2014; Luedtke and van der
Laan, 2016; Künzel et al., 2019, to name but a few). Here, an optimal rule
refers to that which optimizes the average clinical outcome in the patient
population. Note too that ITRs may be inferred from clinical trial and
observational study data alike.

Beyond informing patient treatment decisions, estimated ITRs may provide some
insight on the biological functioning of therapies by classifying covariates
as TEMs. TEM classification generally requires that the rules be inferred
using “interpretable” methods (see Murdoch et al., 2019, for a discussion on
interpretability in machine learning). Examples include penalized linear
modeling procedures and decision trees, methods capable of communicating which
covariates are believed to interact with the treatment in order to predict
outcomes. Importantly, rules estimated by interpretable ITR estimators can be
vetted by domain experts, possibly increasing their trustworthiness from the
perspective of patients and healthcare providers alike.

While many approaches can successfully estimate the optimal ITR in randomized
settings where the number of pre-treatment covariates is small relative to the
number of patients, it is not necessarily the case when there are more
covariates than observations. The challenge is even greater when estimating
this rule from observational study data due to possible confounding. As in
many high-dimensional settings, simplifying assumptions about the complexity
of the estimand must be made to ensure consistent estimation and therefore
reliable predictive performance. Examples include sparsity, such as
subpopulation-specific expression of a particular biomarker (Qian and Murphy,
2011; Tian et al., 2014; Chen et al., 2017; Zhao et al., 2017; Bahamyirou et
al., 2022), or that the outcome is a (partially) linear function of the
covariates and treatment (Tian et al., 2014; Chen et al., 2017). To the best
of our knowledge, there neither exists a comprehensive comparison of popular
ITR estimators in high-dimensional settings or an analysis of these methods’
sensitivity to violations of their assumptions.

Recent developments for uncovering TEMs in high-dimensional data have been
proposed by Boileau et al. (2022, 2023) that might partially alleviate the
difficulty of estimating optimal ITRs in these settings. These authors define
a formal statistical framework for treatment effect modifier variable
importance parameters (TEM-VIPs), permitting the rigorous classification of
individual pre-treatment covariates as TEMs. In the same way that variables
can be filtered prior to fitting traditional regression methods to improve
interpretability and predictive performance, these authors posit that their
TEM-VIP inference procedures may be used with any of the existing ITR
estimators to a similar effect. As a first step, pre-treatment covariates
would be classified as TEMs or non-effect modifiers. Inferred TEMs — in
addition to all other variables known to be clinically relevant — would then
be used by an ITR estimator to learn the treatment rule. This filtering
strategy may be particularly useful when combined with black-box models, which
generally provide greater predictive performance than more rigid parametric
methods at the expense of increased opacity.

Given the lack of practical guidance on ITR estimation in high dimensions and
the potential for the aforementioned TEM-VIPs framework to improve ITR
estimation, we review popular ITR estimators and benchmark them with and
without the use of TEM-VIP-based filtering. To accomplish the latter, a
comprehensive simulation study is performed using a variety of data-generating
processes (DGPs) that are representative of randomized control trials (RCTs)
and observational studies with continuous outcomes, binary treatment
variables, and high-dimensional covariates. ITR estimators are then compared
on the basis of accuracy, interpretability, and computational efficiency.
Recommendations based on these metrics are then provided.

##  2 Problem Formulation

We begin by formally outlining the statistical problem. Consider $n$
independent and identically distributed (i.i.d.) random vectors
$X_{i}=(W_{i},A_{i},Y_{i}^{(0)},Y_{i}^{(1)})\sim P_{X,0}\in\mathcal{M}_{X}$
for, $i=1,\ldots,n$, representing the full, partially-unobserved data. Indices
are omitted for convenience where possible throughout the remainder of the
text. In the full-data random vector, $X$, $W$ is a vector of $p$ pre-
treatment covariates, $A$ is a binary indicator of treatment assignment, and
$Y^{(0)}$ and $Y^{(1)}$ are random variables representing observations’
potential outcomes under assignment to the control and treatment groups,
respectively (Rubin, 1974). We assume that $p$ is approximately equal to or
larger than $n$. Further, $Y^{(0)}$ and $Y^{(1)}$ are either continuous or
binary and it is assumed that large values of $Y$ are beneficial”; this has no
bearing on the definition of the estimand, the optimal ITR. The true full-data
DGP $P_{X,0}$ is typically unknown and is a member of the nonparametric
statistical model $\mathcal{M}_{X}$, the collection of all possible full-data
DGPs.

Only one of the potential outcomes is generally observed for any given
participant; the full data are censored by the treatment assignment mechanism.
However, this characterization of the full-data random vectors permits a
rigorous definition of the optimal ITR and related quantities.

In particular, an intermediate parameter for performing inference about the
optimal ITR is the conditional average treatment effect (CATE):

$\Gamma_{X,0}(W)\equiv\mathbb{E}_{P_{X,0}}\left[Y^{(1)}-Y^{(0)}\big{|}W\right]\;.$
(1)

This functional corresponds to the expected difference of individuals’
potential outcomes conditional on pre-treatment covariates. For any given
participant with $W=w$, a $\Gamma_{X,0}(w)$ larger than zero indicates that,
on average, participants with the same pre-treatment covariates benefit more
from being assigned to the treatment group than to the control group. The
opposite is true when $\Gamma_{X,0}(w)$ is less than zero.

With the CATE in hand, we define the optimal ITR as follows:

$\Psi_{X,0}(W)\equiv I(\Gamma_{X,0}(W)>0)\;.$ (2)

That is, the optimal ITR is an indicator function that equals one when the
CATE is larger than zero and zero otherwise. Proposition 1 follows:

######  Proposition 1.

Let $\mathbf{\Psi}_{X}$ be the collection of all possible ITRs for the full-
data model, i.e., the set of all functions from the covariate space into
$\\{0,1\\}$. Then

$\Psi_{X,0}=\operatorname*{arg\,max}_{\Psi\in\mathbf{\Psi}_{X}}\;\mathbb{E}_{P_{X,0}}\left[Y^{(\Psi(W))}\right]\;.$

###### Proof.

$\displaystyle\Psi_{X,0}$
$\displaystyle=\operatorname*{arg\,max}_{\Psi\in\mathbf{\Psi}_{X}}\;\mathbb{E}_{P_{X,0}}\left[Y^{(\Psi(W))}\right]$
$\displaystyle=\operatorname*{arg\,max}_{\Psi\in\mathbf{\Psi}_{X}}\;\mathbb{E}_{P_{X,0}}\left[\mathbb{E}_{\
P_{X,0}}\left[Y^{(1)}\Psi(W)+Y^{(0)}(1-\Psi(W))\big{|}W\right]\right]$
$\displaystyle=\operatorname*{arg\,max}_{\Psi\in\mathbf{\Psi}_{X}}\;\mathbb{E}_{P_{X,0}}\left[\Psi(W)\mathbb{E}_{P_{X,0}}\left[Y^{(1)}-Y^{(0)}\big{|}W\right]\right]$
$\displaystyle=\operatorname*{arg\,max}_{\Psi\in\mathbf{\Psi}_{X}}\;\mathbb{E}_{P_{X,0}}\left[\Psi(W)\Gamma_{X,0}(W)\right].$

∎

$\Psi_{X,0}$ is therefore optimal in the sense that it maximizes the expected
outcome under an ITR from among all possible ITRs. Note too that it is not
necessarily unique.

Now, as previously mentioned, only one of the potential outcomes is observed
for a given random unit. In place of $\\{X_{i}\\}_{i=1}^{n}$, we observe $n$
i.i.d. random vectors $O=(W,A,Y)\sim P_{0}\in\mathcal{M}$, where $W$ and $A$
are defined as above, and $Y=AY^{(1)}+(1-A)Y^{(0)}$. $P_{0}$ is defined as the
unknown observed-data DGP, and is fully determined by $P_{X,0}$. This
treatment assignment mechanism is denoted by
$\pi_{0}(W)\equiv\mathbb{P}_{P_{0}}[A=1|W]$ throughout the text. Note that
this treatment assignment mechanism is known in RCTs and is unknown in
observational studies. Additionally, we refer to the conditional outcome
$\mathbb{E}_{P_{0}}[Y|A,W]$ by $\mu_{0}(W,A)$.

Regrettably, performing inference about the parameters of Equations (1) and
(2) using the observed data in place of the full data is generally not
possible without making additional assumptions. Identifiability conditions
linking parameters of the observed data DGP to the desired parameters of the
full data DGP are provided in Proposition 2.

######  Proposition 2.

Under the assumptions of no unmeasured confounding between the treatment and
the outcome conditional on $W$, i.e., $Y^{(a)}\perp A\mid W$ for
$a\in\\{0,1\\}$, and that the probability of receiving treatment conditional
on $W$ is bounded away from zero, i.e., $\pi_{0}(W)>0$:

$\begin{split}\Gamma_{0}(W)&\equiv\mu_{0}(W,1)-\mu_{0}(W,0)\\\
&=\Gamma_{X,0}(W)\;,\end{split}$

and

$\begin{split}\Psi_{0}(W)&\equiv I(\Gamma_{0}(W)>0)\\\
&=\Psi_{X,0}(W)\;.\end{split}$

###### Proof.

This result follows immediately from the fact that
$\mu_{0}(W,a)=\mathbb{E}_{P_{X,0}}[Y^{(a)}|W]$ for $a\in\\{0,1\\}$ under the
stated assumptions. ∎

Proposition 2 establishes that the optimal ITR is estimable from the observed
data in the continuous and binary outcome setting. Its assumptions are
satisfied in RCTs: randomization ensures that (1) there are no
treatment–outcome confounders, (2) the treatment assignment mechanism is
known, and (3) the probability of receiving treatment, conditional on
covariates, is generally bounded away from zero. In studies using
observational data, these assumptions are such that the data may be considered
as if generated by an RCT but with an unknown treatment assignment mechanism.

Assuming Proposition 2 is satisfied, we take as target of inference the
optimal ITR, $\Psi_{0}(W)$, which, as a preliminary step, requires the
estimation of the CATE, $\Gamma_{0}(W)$. Possible estimators are presented and
discussed in Section 3.

##  3 Methods

Myriad methods have been developed for estimating the CATE and therefore the
optimal ITR. We review a subset of these approaches capable of performing
inference about the latter in the high-dimensional, binary treatment setting.
We begin with the simplest parametric approach and end with state-of-the-art
nonparametric methods. Throughout, we represent estimators of
$\Gamma_{0}(W),\Psi_{0}(W),\mu_{0}(W,A)$, and $\pi_{0}(W)$ fit using the
empirical distribution $P_{n}$ by $\Gamma_{n}(W),\Psi_{n}(W),\mu_{n}(W,A)$,
and $\pi_{n}(W)$, respectively.

###  3.1 Plug-In Estimators

The most straightforward approach for estimating $\Gamma_{0}$ is to assume
that the conditional expected outcome admits a parametric form, like that of a
generalized linear model (GLM). For example, we might assume that
$\mu_{0}(W,A)=l^{-1}(\alpha_{0}+\beta_{0}A+\gamma_{0}^{\top}W+\delta_{0}^{\top}(AW))$
for some link function $l(\cdot)$. Letting the outcome be continuous and using
the identity link, $l(x)=x$, $\mu_{0}(W,A)$ could be estimated using penalized
regression methods, like the LASSO (Tibshirani, 1996) or the elastic net (Zou
and Hastie, 2005), to produce coefficient estimates
$(\alpha_{n},\beta_{n},\gamma_{n},\delta_{n})$. It then follows that the plug-
in CATE estimator under this linear model is given by

$\begin{split}\Gamma_{n}(W)&\equiv\mu_{n}(W,1)-\mu_{n}(W,0)\\\
&=\beta_{n}+\delta_{n}^{\top}W\;.\end{split}$ (3)

This plug-in estimator can then be used to construct an optimal ITR estimator
under this linear model:

$\Psi_{n}(W)\equiv\text{I}(\Gamma_{n}(W)>0).$ (4)

Analogous estimators can be constructed for scenarios with binary outcomes
using penalized logistic regression to estimate the expected conditional
outcome.

While computationally efficient implementations can be painlessly constructed
in standard software, like the R language and environment for statistical
computing (R Core Team, 2021), and resulting estimators are generally
interpretable, this approach may be hampered by its reliance on strong
parametric assumptions about the functional form of the CATE. When the
parametric form of the conditional expected outcome is misspecified, the
estimator of Equation (3) could be biased, and so too could be the
corresponding ITR estimator of Equation (4). Further, when the covariates are
correlated, the covariates classified as TEMs by penalized regression methods
like the LASSO—variables whose corresponding treatment interaction coefficient
estimates in $\delta_{n}$ are non-zero—may be unreliable (Zhao and Yu, 2006).

We might instead estimate $\mu_{0}(W,A)$ using more flexible machine learning
procedures, like Random Forests (Breiman, 2001) or XGBoost (Chen and Guestrin,
2016). In fact, $\mu_{0}(W,0)$ and $\mu_{0}(W,1)$ can even be estimated
individually on observations in control and treatment conditions,
respectively, using possibly different estimators. Estimators of the CATE and
the accompanying optimal ITR are then constructed similarly to those of
Equations (3) and (4).

These flexible plug-in estimators share similar limitations to their
parametric counterparts, however: if the estimators of $\mu_{0}(W,0)$ and
$\mu_{0}(W,1)$ are not consistent, then neither will their associated CATE
estimator. If estimated using black-box methods, then they may also be less
interpretable and computationally efficient than plug-in estimators obtained
from parametric modeling procedures.

###  3.2 (Augmented) Modified Covariates Estimators

In part motivated by the limitations of the parametric plug-in estimator,
(Tian et al., 2014) developed a semiparametric framework for estimating the
CATE in RCTs. Assume for now that the outcome is a continuous variable and let

$\mu_{0}(W,A)=\alpha_{0}^{\top}f(W)+\beta_{0}^{\top}h(W)\frac{2A-1}{2}\;,$ (5)

where $f$ and $h$ are some functions of the covariates. For simplicity of
presentation, we let $f(W)=h(W)=W$. It follows that
$\Gamma_{0}(W)=\mathbb{E}_{P_{0}}[2Y(2A-1)|W]=\beta_{0}^{\top}W$. Dubbed the
modified outcome method, this approach permits a straightforward estimation of
the CATE using standard ordinary least squares or penalized methods like the
LASSO and elastic net. As (Tian et al., 2014) note, however, this framework
does not easily generalize to DGPs with other kinds of outcome variables.

They instead propose the more general modified covariates method, which relies
on GLMs or the proportional hazards model. Again, first assuming a continuous
outcome and positing the following working model,

$\mu_{0}(W,A)=\beta_{0}^{\top}h(W)\frac{2A-1}{2}\;,$ (6)

Tian et al. (2014) show that the $\beta_{0}$s of Equations (6) and (5) share
identical statistical interpretation and can be estimated in the same manner.
We again let $h(W)=W$. The resulting CATE estimator in this model therefore
takes the form of

$\Gamma_{n}(W)\equiv\beta_{n}^{\top}W\;,$

where $\beta_{n}$ is again estimated via (penalized) linear methods. The
accompanying optimal ITR estimator take the same form as the estimator of
Equation (4).

GLMs can be used to apply this approach to other kinds of outcomes. When $Y$
is binary, the conditional outcome regression can be represented by the
following logistic regression model:

$\mu_{0}(W,A)=\frac{\exp(\beta_{0}^{\top}h(W)(2A-1)/2)}{1+\exp(\beta_{0}^{\top}h(W)(2A-1)/2)}\;.$
(7)

When this working model is well-specified,

$\Gamma_{0}(W)=\frac{\exp(\beta_{0}^{\top}h(W)/2)-1}{\exp(\beta_{0}^{\top}h(W)/2)+1}.$

Estimates of $\beta_{0}$ are again obtained by fitting the working model for
the conditional expected outcome using standard (penalized) regression
approaches. An optimal ITR estimator under this model is constructed as in
Equation (4).

While the modified covariates framework is more flexible than the parametric
approaches based on GLMs, it is generally inefficient. Tian et al. (2014)
therefore proposed an “augmented” modified covariates framework that can be
used to construct CATE estimators with generally smaller variance but which
are asymptotically equivalent. In the continuous outcome scenario with
$h(W)=W$, the resulting CATE estimator is equivalent to $\Gamma_{n}(W)$ of
Equation (3). This simplicity is not shared with nonlinear working models like
that of Equation (7), however (Tian et al., 2014).

Chen et al. (2017) generalized the (augmented) modified covariates methodology
through a loss-based estimation procedure relying on loss functions
incorporating the propensity score. For a continuous outcome, the non-
augmented CATE estimator is defined as

$\Gamma_{n}(W)\equiv\operatorname*{arg\,min}_{\Gamma\in\boldsymbol{\Gamma}_{n}}\;\sum_{i=1}^{n}\frac{(Y_{i}-(2A_{i}-1)\Gamma(W_{i}))^{2}}{(2A_{i}-1)\pi_{0}(W_{i})+1-A_{i}}\;,$

where $\boldsymbol{\Gamma}_{n}$ is the set of possible CATE estimators and
$\pi_{0}(W)$ must be replaced by $\pi_{n}(W)$ in observational studies. This
generalized framework also permits more flexible estimation of $\mu_{0}(W,A)$,
like with XGBoost. More efficient estimation of the CATE and the ITR is
therefore made possible using this modified covariates model in observational
studies, assuming that $\mu_{n}(W,A)$ and $\pi_{n}(W)$ are consistent. Readers
interested in the augmented version of this weighted CATE estimator are
invited to review Chen et al. (2017).

Like the parametric plug-in estimators, these approaches result in biased CATE
estimators when their models are misspecified. This may result in a biased
optimal ITR estimator, too. When the difference in expected outcomes is
estimated using more complex methods, like XGBoost, or penalized GLMs are used
but covariates are correlated (Zhao and Yu, 2006), the interpretability of the
resulting optimal ITR estimate is reduced, too.

A computationally efficient implementations of this methodology is made
available in the personalized R package (Huling and Yu, 2021). We use it
throughout the simulation studies in Sections 4.

###  3.3 Nonparametric Approaches

Much work in the nonparametric estimation literature has been dedicated to
developing estimators of the CATE and the optimal ITR; see Robins (2004);
Robins et al. (2008); Zhang et al. (2012a, b); van der Laan and Luedtke
(2015); Luedtke and van der Laan (2016); Zhao et al. (2017); Bahamyirou et al.
(2022); Kennedy (2022). These estimators, typically relying on nuisance
parameters like $\mu_{0}(W,A)$ and $\pi_{0}(W)$, are generally double robust.
That is, these ITR estimators are consistent estimators under mild regularity
conditions and so long as either $\mu_{n}(W,A)$ or $\pi_{n}(W)$ is consistent.
Given that the treatment assignment rule $\pi_{0}(W)$ is typically known in a
clinical trial, these nonparametric estimators are guaranteed to be consistent
in these settings.

We briefly present the approach inspired by Luedtke and van der Laan (2016).
It relies on the Augmented-Inverse Probability Weighted (AIPW) transform,
defined as

$T(O;\mu_{n},\pi_{n})=\frac{2A-1}{A\pi_{n}(W)+(1-A)(1-\pi_{n}(W))}(Y-\mu_{n}(W,A))+\mu_{n}(W,1)-\mu_{n}(W,0)\;,$

and sometimes referred to as the pseudo-outcome difference.

A nonparametric, “AIPW-based” ITR estimator can then be defined as follows:

$\Gamma_{n}(W)\equiv\operatorname*{arg\,min}_{\Gamma\in\boldsymbol{\Gamma}_{n}}\;\sum_{i=1}^{n}\left(\Gamma(W_{i})-T(O_{i};\mu_{n},\pi_{n})\right)^{2}\;.$

That is, $\Gamma_{n}(W)$ is the squared AIPW-transform error risk minimizer
from among the set of possible CATE estimators. $\Gamma_{n}(W)$ can be
constructed by regressing the pseudo-outcome differences given by
$T(O;\mu_{n},\pi_{n})$ on the covariates $W$. This estimator is a consistent
estimator of $\Gamma_{0}(W)$ when either $\mu_{n}(W,A)$ or $\pi_{n}(W)$ are
consistent (Luedtke and van der Laan, 2016; Bahamyirou et al., 2022). The
accompanying optimal ITR estimator is constructed similarly to the ITR
estimators of Equation (4).

This nonparametric approach makes minimal assumptions about the DGP and is
made all the more attractive by the possible use of flexible machine learning
algorithms to both estimate nuisance parameters and the CATE, curbing the risk
of model misspecification and potentially translating to increased finite
sample precision. We recommend using an estimator based on the Super Learner
framework of van der Laan et al. (2007). Building upon the theory of cross-
validated loss-based estimation (Dudoit and van der Laan, 2003; van der Laan
and Dudoit, 2003; Dudoit and van der Laan, 2005; van der Vaart et al., 2006),
a Super Learner constructs a convex combination of estimators from a pre-
specified library that minimizes the cross-validated risk of a pre-defined
loss. This collection of estimators can be made up of modern machine learning
algorithms, avoiding the need to make strong parametric assumptions about the
target parameter. The resulting estimator converges in probability to the
oracle asymptotically. The oracle corresponds to the estimator that would be
selected for the given dataset if $P_{0}$ were in fact known. The Super
Learner is conveniently implemented in the sl3 R package (Coyle et al., 2021),
which we use in the subsequent simulation study.

This flexibility of this estimation procedure comes at the cost of
computational efficiency, however. Whether the aforementioned asymptotic
properties lead to improved performance in the high-dimensional settings
considered here remains uncertain, too.

The resulting CATE estimator’s level of interpretability also largely depends
on the method used to regress pseudo-outcome differences on the covariates. A
compromise between flexibility and finite sample precision might be achieved
by fitting the pseudo-outcome differences as a function of the covariates
using penalized GLMs (Zhao et al., 2017; Bahamyirou et al., 2022). If
interpretability is not a concern, then a Super Learner composed of flexible
machine learning algorithms will likely provide the best finite sample
performance.

An alternative but related nonparametric approach is that of Causal Random
Forests (Wager and Athey, 2018; Athey et al., 2019). Like the Random Forests
algorithm proposed by Breiman (2001), a large number of decision trees are
constructed on random subsamples of the data. Whereas traditional Random
Forests generate branches in any given tree by maximizing the variance of the
outcome within subgroups formed across random selections of covariates, the
causal version attempts to maximize the variance of the estimated treatment
effects. The CATE is then constructed using a weighted Robinson’s residual-on-
residual regression (Robinson, 1988). This methodology is implemented in the
grf R package (Tibshirani et al., 2022), which we use in the simulation study
of Section 4. Similar to the previously described nonparametric CATE
estimator, $\mu_{n}(W,A)$ and $\pi_{n}(W)$ must be fit. This is typically done
using traditional Random Forests. If either estimator is consistent, then so
too is the CATE estimator given by the Causal Random Forests under mild
conditions on the DGP. Reliably recovering TEMs is however not possible in
high dimensions.

###  3.4 Feature Selection using Treatment Effect Modifier Variable Importance
Parameters

Recent work by Boileau et al. (2022, 2023) discussed and demonstrated how
CATE-based methods aiming to uncover TEMs in high-dimensional RCTs and
observational studies generally fail to reliably uncover true TEMs. This has
been reported elsewhere as well (for example, see the simulation studies of
Tian et al., 2014; Bahamyirou et al., 2022). This suggests that the previously
reviewed CATE estimators may needlessly incorporate uninformative covariates
in their estimation procedures, potentially decreasing their accuracy and
diminishing their interpretability.

Motivated by the need for reliable TEM discovery methodology applicable to
high-dimensional data, Boileau et al. (2022, 2023) proposed a nonparametric
framework for defining TEM-VIPs based on any pathwise differentiable treatment
effect and for subsequently constructing regular asymptotically linear (RAL)
estimators of these parameters. This framework provides formal statistical
guarantees, like Type I error rate control.

Of particular relevance is the absolute TEM-VIP based on the average treatment
effect (ATE), and therefore related to the CATE. Assuming that the expected
value of $\mu_{0}(W,1)-\mu_{0}(W,0)$ conditional on any given $W_{j}$ is
linear in $W_{j}$, this TEM-VIP is defined as the simple linear regression
coefficient obtained by regressing the difference in conditional expected
potential outcomes on a single covariate. For the $j^{\text{th}}$ pre-
treatment covariate $W_{j}$, assuming without loss of generality to be
centered at zero, this parameter is given by

$\frac{\mathbb{E}_{P_{0}}\left[\left(\mu_{0}(W,1)-\mu_{0}(W,0)\right)W_{j}\right]}{\mathbb{E}_{P_{0}}[W_{j}^{2}]}.$

Though the true relationship between this expected outcome difference and the
covariate conditional on said covariate is unlikely to be linear, this
parameter is an informative measure of treatment effect modification: any
value greater than zero is indicative of effect modification, and the
magnitude and sign of the TEM-VIP summarize, respectively, the strength and
direction of the effect.

Boileau et al. (2022, 2023) provide one-step, estimating equation, and
targeted maximum likelihood (TML) estimators of this TEM-VIP that possess
$\mu_{0}(W,A)$ and $\pi_{0}(W)$ as nuisance parameters, and prove that these
estimators are double robust. In observational study scenarios, assuming
$\mu_{n}(W,A)$ and $\pi_{n}(W)$ converge to $\mu_{0}(W,A)$ and $\pi_{0}(W)$,
respectively, at the appropriate nonparametric rate, $o_{P}(n^{-1/4})$, these
TEM-VIP estimators are shown to be RAL. In RCTs, these estimators are
guaranteed to be RAL by virtue of $\pi_{0}(W)$ being known. Recall that RAL
estimators’ sampling distributions are asymptotically Gaussian, permitting
computationally efficient hypothesis testing. Inference procedures for this
TEM-VIP are implemented in the unihtee R package. We use this software in the
simulation studies.

These TEM-VIP estimators can be used in a simple two-stage CATE estimation
procedure. First, the TEM-VIPs are estimated for each covariate, and a
hypothesis test about their modification status is performed. Covariates
classified as TEMs based on a predefined null hypothesis and level of
significance are then used to estimate the CATE.

As previously mentioned, decreasing the number of covariates used by CATE
estimators in the second stage may decrease said estimators’ variance and
increase their accuracy. The two-stage procedure may also be more
computationally efficient than fitting an estimator using the entire set of
covariates. Ranking the estimated TEM-VIPs also adds a layer of
interpretability to the resulting CATE estimate, regardless of the CATE
estimation strategy used in the second stage.

##  4 Simulation Study

###  4.1 Data-Generating Processes and Simulation Details

Recall that, in a continuous outcome setting, we have access to realizations
of $O=(W,A,Y)$, where $W$ is a $p$-dimensional vector of pre-treatment
covariates (and possible confounders), $A$ is a binary treatment indicator,
and $Y$ is the observed continuous outcome. We consider generative models
based on the template below:

$\displaystyle W$ $\displaystyle\sim N(0,\Sigma)$ $\displaystyle A|W$
$\displaystyle\sim\text{Bernoulli}(\pi(W))$ $\displaystyle Y^{(A)}|W$
$\displaystyle\sim N(\mu(W,A),1)$ $\displaystyle Y$
$\displaystyle=AY^{(1)}+(1-A)Y^{(0)}\;.$

Here, $\Sigma$ is some $p\times p$ covariance matrix and $\pi(W)$ and
$\mu(W,A)$ are generic propensity score and conditional outcome functions,
respectively.

Setting $p=500$, we define 16 DGPs using every possible combination of the
following factors:

$\begin{split}\Sigma_{1}&=I_{500\times 500}\\\ \Sigma_{2}&=\text{Block
diagonal}\\\
\end{split}\quad\times\quad\begin{split}\pi_{1}(W)&=\frac{1}{2}\\\
\pi_{2}(W)&=\text{logit}^{-1}\left(\frac{W_{1}+W_{2}+W_{3}+W_{4}}{5}\right)\\\
\end{split}\quad\times\quad\begin{split}\mu_{1}(A,W)&=A+W^{\top}\beta+AW^{\top}\gamma^{(10)}\\\
\mu_{2}(A,W)&=A+W^{\top}\beta+AW^{\top}\gamma^{(50)}\\\
\mu_{3}(A,W)&=W^{\top}\beta+2\;\text{arctan}\left\\{AW^{\top}\gamma^{(10)}\right\\}\\\
\mu_{4}(A,W)&=W^{\top}\beta+2\;\text{arctan}\left\\{AW^{\top}\gamma^{(50)}\right\\}\\\
\end{split}$

where $\beta_{1}=\ldots=\beta_{5}=2$, $\beta_{6}=\ldots=\beta_{500}=0$,
$\gamma^{(10)}_{1}=\ldots=\gamma^{(10)}_{10}=2$,
$\gamma^{(10)}_{11}=\ldots=\gamma^{(10)}_{500}=0$,
$\gamma^{(50)}_{1}=\ldots=\gamma^{(50)}_{50}=1/2$, and
$\gamma^{(50)}_{51}=\ldots=\gamma^{(50)}_{500}=0$. $\Sigma_{2}$ is constructed
by randomly generating $50$ positive definite square symmetric matrices with
diagonal elements equal to one. Details on the generation of these matrices
are provided in the accompanying code.

We note that DGPs using $\pi_{2}$ mimic an RCT, and $\pi_{2}$ is treated as
known in these DGPs. We also highlight that $\gamma^{(10)}$ and
$\gamma^{(50)}$ permit the evaluation of methods in DGPs where the number of
TEMs is sparse and non-sparse, respectively. The effect sizes associated with
each TEM also vary in terms of the conditional expected outcome function:
$\mu_{1}$ produces TEMs with the largest effect sizes, followed by $\mu_{3}$,
$\mu_{2}$, and $\mu_{4}$ in decreasing order.

One hundred learning datasets made up of $n=250,500$ and $1,\\!000$
observations are generated for each DGP. The potential outcomes $Y^{(0)}$ and
$Y^{(1)}$ are unobserved in these datasets. Accompanying each learning set is
a test set made up of $n^{\prime}=100$ observations. These data contain the
potential outcomes for evaluation purposes. The learning datasets are used to
fit the considered CATE and optimal ITR estimators which are subsequently
assessed using the test sets.

###  4.2 Estimators

The CATE-based optimal ITR estimators considered in this simulation study are
listed in Table 1. Details on their fitting procedures are provided, as is a
column indicating whether they provide built-in TEM classification. Indeed,
only estimators relying solely on LASSO linear regressions to estimate the
conditional expected outcome possess this ability: covariates with non-zero
estimated treatment-covariate interaction coefficients are categorized as
TEMs. Note too that all methods employing the LASSO select the regularization
hyperparameter by minimizing the 10-fold cross-validated mean squared error.

CATE Estimator |  Details |  Built-in TEM Classification
---|---|---
LASSO-based Plug-in |  $\mu(W,A)$ is estimated using the LASSO. Linear model coefficients for the treatment indicator, the pre-treatment covariates, and the treatment-covariates interactions are considered. |  Yes
XGBoost-based Plug-in |  XGBoost (Chen and Guestrin, 2016) is used to estimate $\mu(W,A)$. |  No
Modified covariates with LASSO GLMs |  $\mu(W,A)$ is estimated with a LASSO linear regression. $\pi(W)$ is estimated in non-randomized DGPs with a logistic LASSO regression. |  Yes
Augmented modified covariates with LASSO GLMs |  $\mu(W,A)$ is estimated with a LASSO linear regression. $\pi(W)$ is estimated in non-randomized DGPs with a logistic LASSO regression. |  Yes
Modified covariates with XGBoost conditional outcome estimator |  $\mu(W,A)$ is estimated with XGBoost. $\pi(W)$ is estimated in non-randomized DGPs with a logistic LASSO regression. |  No
Augmented modified covariates with XGBoost conditional outcome estimator |  $\mu(W,A)$ is estimated with XGBoost. $\pi(W)$ is estimated in non-randomized DGPs with a logistic LASSO regression. |  No
AIPW-based estimator with LASSO pseudo-outcome regression |  A Super Learner comprised of penalized linear regression (LASSO, ridge, elastic net), Random Forests, and XGBoost is used to estimate $\mu(W,A)$. A Super Learner comprised of penalized logistic regressions (LASSO, ridge, elastic net), Random Forests, and XGBoost is used to estimate $\pi(W)$ when necessary. The difference in predicted pseudo-outcomes is regressed on the covariates using a LASSO linear regression. |  Yes
AIPW-based estimator with Super Learner pseudo-outcome regression |  A Super Learner comprised of penalized linear regression (LASSO, ridge, elastic net), Random Forests, and XGBoost is used to estimate $\mu(W,A)$. A Super Learner comprised of penalized logistic regressions (LASSO, ridge, elastic net), Random Forests, and XGBoost is used to estimate $\pi(W)$ when necessary. The difference in predicted pseudo-outcomes is regressed on the covariates using a Super Learner identical to that used to estimate the conditional expected outcomes. |  No
Causal Random Forests |  Random Forests are used to estimate $\mu(W,A)$ and, when necessary, $\pi(W)$. |  No
Table 1: CATE-based ITR estimators benchmarked in simulation study.

A second “filtered” version each estimator outlined in Table 1, which is
paired with the TEM-VIP framework of Boileau et al. (2022, 2023), is also
included in this benchmark. As previously described in Section 3, joint
inference about TEM-VIPs based on the CATE is performed to identify TEMs on
each learning dataset replicate. These predicted TEMs are then used to train
the individual CATE-based ITR estimators.

We use the nonparametric double-robust one-step estimator to perform inference
about the absolute ATE-based TEM-VIPs (Boileau et al., 2022, 2023). Recall
that in observational studies, this filtering procedure requires that two
nuisance parameters be estimated: $\mu(W,A)$ and $\pi(W)$. They are fit in the
observational study DGPs considered here with the same Super Learners used by
the AIPW-based estimators listed in Table 1. In the RCT-like DGPs, $\pi(W)$
need not be estimated, and $\mu(W,A)$ is estimated with a LASSO regression
that includes main terms for the treatment and covariates, as well as all
possible treatment-covariate interaction terms. Again, regardless of whether a
DGP mimics an observational study or an RCT, this one-step estimator is
asymptotically Gaussian under mild conditions. Wald-type confidence intervals
are constructed for each pre-treatment covariate’s TEM-VIP, and any variable
with an Benjamini-Hochberg FDR-adjusted $p$-value (Benjamini and Hochberg,
1995) lower than $5\%$ is classified as a TEM.

###  4.3 Estimator Performance Metrics

Estimators are evaluated using three metrics: accuracy of the resulting ITR-
based classification on new data, interpretability, and computational
efficiency.

#### Accuracy:

For a generic ITR estimate $\psi_{b}(W)$ produced by fitting a generic ITR
estimator $\Psi_{n}(W)$ to the $b^{\text{th}}$ learning dataset, we compute
the mean outcome over the $b^{\text{th}}$ test set based on the resulting
estimate’s classification as

$\frac{1}{n^{\prime}}\sum_{i=1}^{n^{\prime}}Y_{i}^{(\psi_{b}(W_{i}))}\;,$

where we recall that each test set is made up of $n^{\prime}=100$
observations. $\Psi_{n}(W)$’s performance across the $B=100$ test set
replicates is then summarized as

$\frac{1}{B}\sum_{b=1}^{B}\left(\frac{1}{n^{\prime}}\sum_{i=1}^{n^{\prime}}Y_{i}^{(\psi_{b}(W_{i}))}\right)\;.$

This metric, corresponding to the empirical mean of the potential outcomes
under the ITR $\Psi_{n}(W)$, is therefore an estimator of
$\mathbb{E}_{P_{X,0}}[Y^{(\Psi_{X,0}(W))}]$, assuming the conditions of
Proposition 2 are satisfied. ITR estimators that produce the largest value are
deemed the most “accurate”—they accurately assign new observations to the
correct treatment group.

#### Interpretability:

The non-zero entries in $\gamma^{(10)}$ and $\gamma^{(50)}$ correspond to the
TEMs of their respective DGPs, and are therefore the driving force behind
heterogeneous treatment effects. For any given learning set replicate, the
false discovery, true negative, and true positive proportions with respect to
TEM status are computed for all estimators capable of classifying pre-
treatment covariates. These proportions are then averaged across replicates to
produce the empirical false discovery rate (FDR), true negative rate (TNR),
and true positive rate (TPR). Optimal ITR estimators are deemed interpretable
if the empirical false discovery rate (FDR) approximately achieves the nominal
Type I error rate of $5\%$ while producing TNRs and TPRs near 100%.

#### Computational Efficiency:

The time to fit each estimator, using a single core, to each learning dataset
is recorded and considered a surrogate for overall computational efficiency.
The mean fit time is then computed across sample sizes and DGPs. Smaller times
indicate improved computational efficiency. We note that some of the
estimators implementations permit parallelization, which generally decreases
the required time. We nevertheless fit all estimators in serial to ensure a
comparison that is relevant and informative for an audience with access to a
diverse set of computing environments.

###  4.4 Results

Figure 1:  Accuracy: Mean ITR outcome computed over the 100 test set
replicates vs. sample size, for each combination of DGP. The dotted line
corresponds to the mean outcome under the optimal ITR defined in Proposition
1, which is approximated using a Monte Carlo procedure for each DGP.

#### Accuracy:

No estimator uniformly dominates the others (Figure 1). In the DGPs with
sparse TEMs—regardless of whether treatment is randomized or of the
covariates’ correlation structure—the filtered and non-filtered versions of
the LASSO-based plug-in, AIPW-based, and augmented modified covariates
estimators are the top performers. These estimators produce empirical means
that are approximately equal to the population mean under the optimal ITR in
the linear model DGPs at all sample sizes and in the non-linear DGPs when
$n=1,\\!000$. When $n=250$, the non-filtered versions of these estimators
generally perform marginally better than their filtered counterparts (Figure
A1). The reverse is true as sample size increases.

In the non-sparse settings, the non-filtered plug-in LASSO, AIPW-based, and
LASSO-based augmented modified covariates estimators produce the best
empirical mean outcomes on average (Figures 1, A1). The empirical means
obtained under these estimators are near the population mean under the optimal
ITRs as sample size increases. The performance of the filtered versions of
these estimators’ generally converges with increasing sample size to that of
their non-filtered counterparts. An exception to this trend is observed in the
DGPs with non-linear conditional expected outcomes and whose covariates’
covariance matrix is $I_{500\times 500}$.

The XGBoost-based plug-in, modified covariates, and Causal Random Forest
estimators are generally less accurate than other estimators considered here,
though their performance is improved by the covariate filtering procedure.

Figure 2:  Interpretability, FDR: Empirical FDR computed over the 100 learning
set replicates vs. sample size for each combination of DGP. The dotted line
corresponds to desired nominal Type I error rate of $5\%$.

#### Interpretability:

As illustrated in Boileau et al. (2022, 2023), only ITR estimators employing
the TEM-VIP-based filtering procedure reliably control the FDR near the
nominal $5\%$ level regardless of the treatment assignment mechanism (Figure
2). Indeed, these estimators produce empirical FDR equal to $\approx 5\%$ in
all DGPs for sample sizes equal to or larger than $500$. All other estimators
capable of classifying TEMs produce empirical FDRs ranging from $40\%$ to
$80\%$. However, we find that this FDR control sometimes results in a worse
empirical TPR in non-sparse settings (Figure A2). With respect to the
empirical TNR, only the filtered estimators produce near-perfect results, and
this regardless of sample size (Figure A3).

Figure 3:  Computational Efficiency: Fit time computed over the 100 learning
set replicates vs. sample size for each combination of DGP. The dotted line
corresponds to the minimum mean fit time observed in each DGP.

#### Computational Efficiency:

In RCT DGPs, the non-filtered LASSO-based plug-in estimator generally has the
lowest mean fit time across all sample sizes and DGPs, followed closely by the
non-filtered (augmented) modified covariates estimators using the LASSO and
the non-filtered XGBoost plug-in estimator (Figure 3). The non-filtered AIPW-
based estimators are the slowest estimators to produce a rule in all DGPs and
sample sizes. All estimators using the TEM-VIP-based filtering procedure
exhibit similar mean fit times that are several orders of magnitude larger
than that of the non-filtered LASSO-based plug-in estimator. These filtered
estimators’ mean fit times are sometimes similar, however, to their non-
filtered counterparts when $n=1,\\!000$, as is the case with Causal Random
Forests and the augmented modified covariates estimators using XGBoost.
Further, the filtered AIPW-based estimators have lower mean fit times than
their non-filtered versions.

Similar patterns are observed in the observational study DGPs’ simulations.
The non-filtered plug-in estimator using the LASSO is again the most
computationally efficient procedure, followed by the non-filtered plug-in
estimator using XGBoost (Figure 3). The non-filtered (augmented) modified
covariates estimators and non-filtered the Causal Random Forests exhibit
similar mean fit times which are approximately two orders of magnitude larger
than the fastest estimator. The nonparametric AIPW-based estimators are again
the slowest of the non-filtered estimators, though they produce estimated
rules more quickly than any of the estimators relying on the TEM-VIP-based
filtering strategy.

##  5 Discussion

The simulation results suggest that no ITR estimator uniformly dominates all
others with respect to accuracy, interpretability, and computational
efficiency in the continuous outcome, binary treatment setting. We posit that
the same is true in binary outcome, binary treatment settings. A trade-off
must therefore be made when selecting a method with which to learn a treatment
rule. Given that the primary goal of learning an ITR is to optimize patient
benefit—that is, to learn the optimal ITR—we hold that only accurate
estimators should be considered. Selecting a procedure therefore comes down to
favoring either interpretability or computational efficiency as a secondary
goal.

If more importance is placed on interpretability than computational
efficiency, the TEM-VIP-based filtering procedure is a virtual necessity. None
of the non-filtered estimators considered here reliably control the empirical
FDR or produce near-perfect empirical TPRs (Figures 2, A3). Even when
empirical TPR optimization is required, filtered estimators generally perform
as well as non-filtered estimators (Figure A2). In sparse TEM scenarios, the
most accurate estimators are generally filtered (Figure A1). In non-sparse TEM
settings, the filtered LASSO-based plug-in and AIPW-based estimators are
typically only marginally worse than the most accurate non-filtered
estimators.

However, the use of this filtering procedure may be constrained by available
computing power. In RCTs, this filtering procedure makes fitting ITR
estimators slower than their non-filtered versions in small sample sizes
(Figure 3)—with the exception of the AIPW-based ITR estimators. However, this
gap decreases in larger sample sizes: the mean fit times for the filtered
estimators are approximately constant, while the non-filtered estimators’ fit
times follow an approximately exponential function of sample size. In
observational studies, where more computationally intensive nuisance parameter
estimators are generally required to ensure reliable TEM classification, the
filtered estimators are generally several orders of magnitude slower than
their non-filtered counterparts.

When there is limited time, effort, or computing resources available to learn
an ITR, the simulation study results suggest that the non-filtered LASSO-based
plug-in estimator or the non-filtered augmented modified covariates approach
using XGBoost should be used for analyzing data collected in an RCT. In
observational studies, the LASSO-based plug-in estimator is the most
computationally efficient of the sufficiently accurate estimators. We again
emphasize, however, that the computational efficiency of these procedures
comes at the cost interpretability: these estimators cannot reliably
distinguish TEMs from non-modifiers.

##  6 Recommendations

We summarize our discussion of the simulation results with concise
recommendations about the choice of ITR estimators based on desired operating
characteristics.

#### Accurate and interpretable:

The filtered LASSO-based plug-in and AIPW-based estimators generally produce
accurate ITR estimates while reliably recovering TEMs. These estimators are
computationally intensive, however. When computational resources are
available, parallelized implementations of these estimators will decrease the
time required to produce a treatment rule estimate.

#### Accurate and computationally efficient:

The LASSO-based plug-in estimator is among the most accurate in our simulation
study, providing empirical evidence that it is robust to model
misspecification while being exceptionally computationally efficient. However,
as previously mentioned, this estimator’s built-in feature selection
capabilities should not be used for TEM recovery.

We emphasize that while this simulation study is comprehensive, practitioners
may nevertheless benefit from performing bespoke benchmarks of their own. The
code used to generate the simulations results of Section 4, found on GitHub,
may serve as foundation for extensions and modifications.

## Code

Code used to generate the results is available in the
PhilBoileau/pub_guidance-on-ITR-estimation-in-HD GitHub repository.

## Appendix

Figure A1:  Accuracy: Mean ITR outcome ranks computed over the 100 test set
replicates vs. sample size for each combination of DGP. The dotted line
corresponds to the desired rank, one. Figure A2:  Interpretability, TPR:
Empirical TPR computed over the 100 learning set replicates vs. sample size
for each combination of DGP. The dotted line corresponds to desired TPR of
$100\%$. Figure A3:  Interpretability, TNR: Empirical TNR computed over the
100 learning set replicates vs. sample size for each combination of DGP. The
dotted line corresponds to desired TNR of $100\%$.

## References

  * Athey et al. [2019] S. Athey, J. Tibshirani, and S. Wager.  Generalized random forests.  _The Annals of Statistics_ , 47(2):1148–1178, 4 2019.
  * Bahamyirou et al. [2022] A. Bahamyirou, M. E. Schnitzer, E. H. Kennedy, L. Blais, and Y. Yang.  Doubly robust adaptive lasso for effect modifier discovery.  _The International Journal of Biostatistics_ , 2022.  doi: doi:10.1515/ijb-2020-0073.  URL https://doi.org/10.1515/ijb-2020-0073.
  * Benjamini and Hochberg [1995] Y. Benjamini and Y. Hochberg.  Controlling the false discovery rate: A practical and powerful approach to multiple testing.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 57(1):289–300, 1995.  ISSN 00359246.  URL http://www.jstor.org/stable/2346101.
  * Boileau et al. [2022] P. Boileau, N. T. Qi, M. J. van der Laan, S. Dudoit, and N. Leng.  A flexible approach for predictive biomarker discovery.  _Biostatistics_ , 07 2022.  ISSN 1465-4644.  doi: 10.1093/biostatistics/kxac029.  URL https://doi.org/10.1093/biostatistics/kxac029.  kxac029.
  * Boileau et al. [2023] P. Boileau, N. Leng, N. S. Hejazi, M. van der Laan, and S. Dudoit.  A nonparametric framework for treatment effect modifier discovery in high dimensions, 2023.
  * Breiman [2001] L. Breiman.  Random forests.  _Machine Learning_ , 45(1):5–32, 2001.
  * Chen et al. [2017] S. Chen, L. Tian, T. Cai, and M. Yu.  A general statistical framework for subgroup identification and comparative treatment scoring.  _Biometrics_ , 73(4):1199–1209, 2017.  doi: https://doi.org/10.1111/biom.12676.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12676.
  * Chen and Guestrin [2016] T. Chen and C. Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939785.  URL http://doi.acm.org/10.1145/2939672.2939785.
  * Coyle et al. [2021] J. R. Coyle, N. S. Hejazi, I. Malenica, R. V. Phillips, and O. Sofrygin.  _sl3: Modern Pipelines for Machine Learning and Super Learning_ , 2021.  URL https://doi.org/10.5281/zenodo.1342293.  R package version 1.4.2.
  * Dudoit and van der Laan [2003] S. Dudoit and M. J. van der Laan.  Asymptotics of cross-validated risk estimation in estimator selection and performance assessment.  Working Paper 126, University of California, Berkeley, Berkeley, 2003\.  URL https://biostats.bepress.com/ucbbiostat/paper126.
  * Dudoit and van der Laan [2005] S. Dudoit and M. J. van der Laan.  Asymptotics of cross-validated risk estimation in estimator selection and performance assessment.  _Statistical Methodology_ , 2(2):131–154, 2005\.
  * Huling and Yu [2021] J. D. Huling and M. Yu.  Subgroup identification using the personalized package.  _Journal of Statistical Software_ , 98(5):1–60, 2021.  doi: 10.18637/jss.v098.i05.  URL https://www.jstatsoft.org/index.php/jss/article/view/v098i05.
  * Kennedy [2022] E. H. Kennedy.  Towards optimal doubly robust estimation of heterogeneous causal effects, 2022.
  * Künzel et al. [2019] S. R. Künzel, J. S. Sekhon, P. J. Bickel, and B. Yu.  Metalearners for estimating heterogeneous treatment effects using machine learning.  _Proceedings of the National Academy of Sciences_ , 116(10):4156–4165, 2019.
  * Luedtke and van der Laan [2016] A. R. Luedtke and M. J. van der Laan.  Super-learning of an optimal dynamic treatment rule.  _The International Journal of Biostatistics_ , 12(1):305–332, 2016.  doi: doi:10.1515/ijb-2015-0052.  URL https://doi.org/10.1515/ijb-2015-0052.
  * Murdoch et al. [2019] W. J. Murdoch, C. Singh, K. Kumbier, R. Abbasi-Asl, and B. Yu.  Definitions, methods, and applications in interpretable machine learning.  _Proceedings of the National Academy of Sciences_ , 116(44):22071–22080, 2019.
  * Qian and Murphy [2011] M. Qian and S. A. Murphy.  Performance guarantees for individualized treatment rules.  _The Annals of Statistics_ , 39(2):1180 – 1210, 2011.  doi: 10.1214/10-AOS864.  URL https://doi.org/10.1214/10-AOS864.
  * R Core Team [2021] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2021.  URL https://www.R-project.org/.
  * Robins et al. [2008] J. Robins, L. Orellana, and A. Rotnitzky.  Estimation and extrapolation of optimal treatment and testing strategies.  _Statistics in Medicine_ , 27(23):4678–4721, 2008\.  doi: https://doi.org/10.1002/sim.3301.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.3301.
  * Robins [2004] J. M. Robins.  _Optimal Structural Nested Models for Optimal Sequential Decisions_ , pages 189–326.  Springer New York, New York, NY, 2004.
  * Robinson [1988] P. M. Robinson.  Root-n-consistent semiparametric regression.  _Econometrica_ , 56:931–954, 1988.
  * Rubin [1974] D. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688–701, 1974.
  * Tian et al. [2014] L. Tian, A. A. Alizadeh, A. J. Gentles, and R. Tibshirani.  A simple method for estimating interactions between a treatment and a large number of covariates.  _Journal of the American Statistical Association_ , 109(508):1517–1532, 2014.  doi: 10.1080/01621459.2014.951443.  URL https://doi.org/10.1080/01621459.2014.951443.  PMID: 25729117.
  * Tibshirani et al. [2022] J. Tibshirani, S. Athey, E. Sverdrup, and S. Wager.  _grf: Generalized Random Forests_ , 2022.  URL https://CRAN.R-project.org/package=grf.  R package version 2.2.1.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58(1):267–288, 1996.  ISSN 00359246.  URL http://www.jstor.org/stable/2346178.
  * van der Laan and Dudoit [2003] M. J. van der Laan and S. Dudoit.  Unified cross-validation methodology for selection among estimators and a general cross-validated adaptive epsilon-net estimator: Finite sample oracle inequalities and examples.  Working Paper 130, University of California, Berkeley, Berkeley, 2003\.  URL https://biostats.bepress.com/ucbbiostat/paper130/.
  * van der Laan and Luedtke [2015] M. J. van der Laan and A. R. Luedtke.  Targeted learning of the mean outcome under an optimal dynamic treatment rule.  _Journal of Causal Inference_ , 3(1):61–95, 2015\.  doi: doi:10.1515/jci-2013-0022.  URL https://doi.org/10.1515/jci-2013-0022.
  * van der Laan et al. [2007] M. J. van der Laan, E. C. Polley, and A. E. Hubbard.  Super learner.  _Statistical Applications in Genetics and Molecular Biology_ , 6(1), 2007.  doi: doi:10.2202/1544-6115.1309.  URL https://doi.org/10.2202/1544-6115.1309.
  * van der Vaart et al. [2006] A. W. van der Vaart, S. Dudoit, and M. J. van der Laan.  Oracle inequalities for multi-fold cross validation.  _Statistics and Decisions_ , 24:351–371, 2006.
  * Wager and Athey [2018] S. Wager and S. Athey.  Estimation and inference of heterogeneous treatment effects using Random Forests.  _Journal of the American Statistical Association_ , 113(523):1228–1242, 2018.  doi: 10.1080/01621459.2017.1319839.  URL https://doi.org/10.1080/01621459.2017.1319839.
  * Zhang et al. [2012a] B. Zhang, A. A. Tsiatis, M. Davidian, M. Zhang, and E. Laber.  Estimating optimal treatment regimes from a classification perspective.  _Stat_ , 1(1):103–114, 2012a.  doi: https://doi.org/10.1002/sta.411.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sta.411.
  * Zhang et al. [2012b] B. Zhang, A. A. Tsiatis, E. B. Laber, and M. Davidian.  A robust method for estimating optimal treatment regimes.  _Biometrics_ , 68(4):1010–1018, 2012b.  doi: https://doi.org/10.1111/j.1541-0420.2012.01763.x.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1541-0420.2012.01763.x.
  * Zhao and Yu [2006] P. Zhao and B. Yu.  On model selection consistency of lasso.  _Journal of Machine Learning Research_ , 7(90):2541–2563, 2006.  URL http://jmlr.org/papers/v7/zhao06a.html.
  * Zhao et al. [2017] Q. Zhao, D. S. Small, and A. Ertefaie.  Selective inference for effect modification via the lasso, 2017.  URL https://arxiv.org/abs/1705.08020.
  * Zou and Hastie [2005] H. Zou and T. Hastie.  Regularization and variable selection via the elastic net.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 67(2):301–320, 2005.