It turns out that both approaches (my initial attempt and another based on suggestions in the comment section) in my original post give the same answer. I will outline both methods here for a complete answer to the question.
Here, $\text{Gamma}(n,\theta)$ means the gamma density $f(y)=\frac{\theta^n}{\Gamma(n)}e^{-\theta y}y^{n-1}\mathbf1_{y>0}$ where $\theta,n>0$, and $\text{Exp}(\theta)$ denotes an exponential distribution with mean $1/\theta$, ($\theta>0$). Clearly, $\text{Exp}(\theta)\equiv\text{Gamma}(1,\theta)$.
Since $T=\sum_{i=1}^n\ln X_i$ is complete sufficient for $\theta$ and $\mathbb E(X_1)=\frac{\theta}{1+\theta}$, by the Lehmann-Scheffe theorem $\mathbb E(X_1\mid T)$ is the UMVUE of $\frac{\theta}{1+\theta}$.
So we have to find this conditional expectation.
We note that $X_i\stackrel{\text{i.i.d}}{\sim}\text{Beta}(\theta,1)\implies-\ln X_i\stackrel{\text{i.i.d}}{\sim}\text{Exp}(\theta)\implies-T\sim\text{Gamma}(n,\theta)$.
Method I:
Let $U=-\ln X_1$ and $V=-\sum_{i=2}^n\ln X_i$, so that $U$ and $V$ are independent. Indeed, $U\sim\text{Exp}(\theta)$ and $V\sim\text{Gamma}(n-1,\theta)$, implying $U+V\sim\text{Gamma}(n,\theta)$.
So, $\mathbb E(X_1\mid \sum_{i=1}^n\ln X_i=t)=\mathbb E(e^{-U}\mid U+V=-t)$.
Now we find the conditional distribution of $U\mid U+V$.
For $n>1$ and $\theta>0$, joint density of $(U,V)$ is
\begin{align}f_{U,V}(u,v)&=\theta e^{-\theta u}\mathbf1_{u>0}\frac{\theta^{n-1}}{\Gamma(n-1)}e^{-\theta v}v^{n-2}\mathbf1_{v>0}\\&=\frac{\theta^n}{\Gamma(n-1)}e^{-\theta(u+v)}v^{n-2}\mathbf1_{u,v>0}\end{align}
Changing variables, it is immediate that the joint density of $(U,U+V)$ is $$f_{U,U+V}(u,z)=\frac{\theta^n}{\Gamma(n-1)}e^{-\theta z}(z-u)^{n-2}\mathbf1_{0<u<z}$$
Let $f_{U+V}(\cdot)$ be the density of $U+V$. Thus conditional density of $U\mid U+V=z$ is \begin{align}f_{U\mid U+V}(u\mid z)&=\frac{f_{U,U+V}(u,z)}{f_{U+V}(z)}\\&=\frac{(n-1)(z-u)^{n-2}}{z^{n-1}}\mathbf1_{0<u<z}\end{align}
Therefore, $\displaystyle\mathbb E(e^{-U}\mid U+V=z)=\frac{n-1}{z^{n-1}}\int_0^z e^{-u}(z-u)^{n-2}\,\mathrm{d}u$.
That is, the UMVUE of $\frac{\theta}{1+\theta}$ is $\displaystyle\mathbb E(X_1\mid T)=\frac{n-1}{(-T)^{n-1}}\int_0^{-T} e^{-u}(-T-u)^{n-2}\,\mathrm{d}u\tag{1}$
Method II:
As $T$ is a complete sufficient statistic for $\theta$, any unbiased estimator of $\frac{\theta}{1+\theta}$ which is a function of $T$ will be the UMVUE of $\frac{\theta}{1+\theta}$ by the Lehmann-Scheffe theorem. So we proceed to find the moments of $-T$, whose distribution is known to us. We have,
$$\mathbb E(-T)^r=\int_0^\infty y^r\theta^n\frac{e^{-\theta y}y^{n-1}}{\Gamma(n)}\,\mathrm{d}y=\frac{\Gamma(n+r)}{\theta^r\,\Gamma(n)},\qquad n+r>0$$
Using this equation we obtain unbiased estimators (and UMVUE's) of $1/\theta^r$ for every integer $r\ge1$.
Now for $\theta>1$, we have $\displaystyle\frac{\theta}{1+\theta}=\left(1+\frac{1}{\theta}\right)^{-1}=1-\frac{1}{\theta}+\frac{1}{\theta^2}-\frac{1}{\theta^3}+\cdots$
Combining the unbiased estimators of $1/\theta^r$ we obtain $$\mathbb E\left(1+\frac{T}{n}+\frac{T^2}{n(n+1)}+\frac{T^3}{n(n+1)(n+2)}+\cdots\right)=1-\frac{1}{\theta}+\frac{1}{\theta^2}-\frac{1}{\theta^3}+\cdots$$
That is, $$\mathbb E\left(\sum_{r=0}^\infty \frac{T^r}{n(n+1)...(n+r-1)}\right)=\frac{\theta}{1+\theta}$$
So assuming $\theta>1$, the UMVUE of $\frac{\theta}{1+\theta}$ is $g(T)=\displaystyle\sum_{r=0}^\infty \frac{T^r}{n(n+1)...(n+r-1)}\tag{2}$
I am not certain about the case $0<\theta<1$ in the second method.
According to Mathematica, equation $(1)$ has a closed form in terms of the incomplete gamma function. And in equation $(2)$, we can express the product $n(n+1)(n+2)...(n+r-1)$ in terms of the usual gamma function as $n(n+1)(n+2)...(n+r-1)=\frac{\Gamma(n+r)}{\Gamma(n)}$. This perhaps provides the apparent connection between $(1)$ and $(2)$.
Using Mathematica I could verify that $(1)$ and $(2)$ are indeed the same thing.
Your reasoning is mostly correct.
The joint density of the sample $(X_1,X_2,\ldots,X_n)$ is
\begin{align}
f_{\theta}(x_1,x_2,\ldots,x_n)&=\frac{\theta^n}{\left(\prod_{i=1}^n (1+x_i)\right)^{1+\theta}}\mathbf1_{x_1,x_2,\ldots,x_n>0}\qquad,\,\theta>0
\\\\\implies \ln f_{\theta}(x_1,x_2,\ldots,x_n)&=n\ln(\theta)-(1+\theta)\sum_{i=1}^n\ln(1+x_i)+\ln(\mathbf1_{\min_{1\le i\le n} x_i>0})
\\\\\implies\frac{\partial}{\partial \theta}\ln f_{\theta}(x_1,x_2,\ldots,x_n)&=\frac{n}{\theta}-\sum_{i=1}^n\ln(1+x_i)
\\\\&=-n\left(\frac{\sum_{i=1}^n\ln(1+x_i)}{n}-\frac{1}{\theta}\right)
\end{align}
Thus we have expressed the score function in the form
$$\frac{\partial}{\partial \theta}\ln f_{\theta}(x_1,x_2,\ldots,x_n)=k(\theta)\left(T(x_1,x_2,\ldots,x_n)-\frac{1}{\theta}\right)\tag{1}$$
, which is the equality condition in the Cramér-Rao inequality.
It is not difficult to verify that $$E(T)=\frac{1}{n}\sum_{i=1}^n\underbrace{E(\ln(1+X_i))}_{=1/\theta}=\frac{1}{\theta}\tag{2}$$
From $(1)$ and $(2)$ we can conclude that
- The statistic $T(X_1,X_2,\ldots,X_n)$ is an unbiased estimator of $1/\theta$.
- $T$ satisfies the equality condition of the Cramér-Rao inequality.
These two facts together imply that $T$ is the UMVUE of $1/\theta$.
The second bullet actually tells us that variance of $T$ attains the Cramér-Rao lower bound for $1/\theta$.
Indeed, as you have shown,
$$E_{\theta}\left[\frac{\partial^2}{\partial\theta^2}\ln f_{\theta}(X_1)\right]=-\frac{1}{\theta^2}$$
This implies that the information function for the whole sample is $$I(\theta)=-nE_{\theta}\left[\frac{\partial^2}{\partial\theta^2}\ln f_{\theta}(X_1)\right]=\frac{n}{\theta^2}$$
So the Cramér-Rao lower bound for $1/\theta$ and hence the variance of the UMVUE is
$$\operatorname{Var}(T)=\frac{\left[\frac{d}{d\theta}\left(\frac{1}{\theta}\right)\right]^2}{I(\theta)}=\frac{1}{n\theta^2}$$
Here we have exploited a corollary of the Cramér-Rao inequality, which says that for a family of distributions $f$ parametrised by $\theta$ (assuming regularity conditions of CR inequality to hold), if a statistic $T$ is unbiased for $g(\theta)$ for some function $g$ and if it satisfies the condition of equality in the CR inequality, namely $$\frac{\partial}{\partial\theta}\ln f_{\theta}(x)=k(\theta)\left(T(x)-g(\theta)\right)$$, then $T$ must be the UMVUE of $g(\theta)$. So this argument does not work in every problem.
Alternatively, using the Lehmann-Scheffe theorem you could say that $T=\frac{1}{n}\sum_{i=1}^{n} \ln(1+X_i)$ is the UMVUE of $1/\theta$ as it is unbiased for $1/\theta$ and is a complete sufficient statistic for the family of distributions. That $T$ is compete sufficient is clear from the structure of the joint density of the sample in terms of a one-parameter exponential family. But variance of $T$ might be a little tricky to find directly.
Best Answer
From your previous question, you already have the complete sufficient statistic:
$$T(\mathbf{X}) = \sum_{i=1}^n \ln(1+X_i).$$
The simplest way to find the UMVUE estimator for $\theta$ is to appeal to the Lehmann-Scheffé theorem, which says that any unbiased estimator of $\theta$ which is a function of $T$ is the unique UMVUE.
To find an estimator with these properties, let $T_i = \ln(1+X_i)$ and observe that $T_i \sim \text{Exp}(\theta)$ so that $T \sim \text{Gamma}(n,\theta)$. Hence, we can use the complete sufficient statistic to form the estimator:
$$\hat{\theta}(\mathbf{X}) = \frac{n-1}{T(\mathbf{X})} = \frac{n-1}{\sum_{i=1}^n \ln(1+X_i)} \sim (n-1) \cdot \text{Inv-Gamma}(n,\theta).$$
From the known moments of the inverse gamma distribution, we have $\mathbb{E}(\hat{\theta}) = \theta$ so we have found an unbiased estimator that is a function of the complete sufficient statistic. The Lehmann-Scheffé theorem ensures that our estimator is UMVUE for $\theta$. The variance of the UMVUE can easily be found by appeal to the moments of the inverse gamma distribution.