Existence of radial limits of products of certain power series and $1-x$
As for probabilistic approaches... let's just say they're not my cup of tea Well, that certainly has to be fixed if one "thinks of himself as an analyst". Let me give you a quick overview of the things most relevant to your current project as a series of exercises.
1) Read https://en.wikipedia.org/wiki/Hoeffding%27s_inequality
2) Look up the proof of the Bernstein inequality $$ \|P'\|_\infty\le Cn\|P\|_\infty $$ for trigonometric polynomials $P(t)=\sum_{k=-n}^n c_ke^{ikt}$ of degree $n$. You do not need the sharp value of $C$.
3) Derive from it that $\|P\|_\infty\le C\max_{k=0}^{100n}|P(e^{(2\pi i k)/(100n)}|$
4) Using Hoeffding and 3), show that if $\xi_j$ are mean zero independent Bernoulli that are $p_j-1$ with probability $p_j$ and $p_j$ with probability $(1-p_j)$, then $$ P(\|\frac 1N\sum_{j=1}^N\xi_jz^j\|_\infty\ge N^{-0.1})\le N^{-0.1} $$ for large $N$.
5) Using Borel-Cantelli and the fact that Chesaro means of a bounded sequence can be read from any subsequence $N_k$ with $N_{k+1}/N_k\to 1$ as $k\to\infty$, show that with probability $1$, we have $$ \lim_{N\to\infty}\frac 1N \sum_{j=1}^N\xi_jz^j=0 $$ for all $z\in\mathbb T$ simultaneously.
6) Consider now the random set $V$ to which each particular integer $j$ belongs independently with probability $a_j\in[0,1]$. Show that with probability $1$, for every $t$, $$ \lim_{N\to\infty}\frac 1N\left|\sum_{v\in V, v\le N}z^v-\sum_{v=1}^N a_vz^v\right|=0 $$
Thus, it does not matter if you consider sums over sets or arbitrary sums with coefficients in $[0,1]$, but with the latter you have way more freedom in various constructions. For instance,
7) Construct a series that gives you a non-zero limit at $0$, $z$, $\bar z$ but $0$ limit everywhere else.
8) Construct a series that gives you non-zero limits at any prescribed self-conjugate set of "rational points" containing $t=0$ and zero limits everywhere else.
Etc.
There are lots of different questions going on here, but here are a few observations.
First of all, if you change $V$ by any set with (upper) density zero, then $c_V(t)$ doesn't change for any value of $t$. So for any $T_V$ there are uncountably many possible choices of $V$, and so almost all of them will be transcendental functions. (For an explicit example where $T_V = [0,1)$, let $V_S$ consist of the set of all integers minus any subset $S$ consisting of powers of $2$. The set of possible $S$ is uncountable.)
Second, it's also relatively easy to construct $V$ such that $T_V$ is missing infinitely many rational numbers. Let's start out with an example which is just missing $1/3$ and $2/3$. The idea is to take $V$ to consist of a string of numbers which are $1$ mod $3$, then a much longer string of numbers which is $2$ mod $3$, then a much much longer string which is $1$ mod $3$ and so on. That is,
$$\varsigma(x) = \frac{x}{1-x^3} + x^{3a_1} \frac{x^2 - x}{1-x^3} + x^{3 a_2} \frac{x - x^2}{1 - x^3} + x^{3 a_3} \frac{x^2 - x}{1-x^3} + \ldots $$
Note that
$$\varsigma(x) = \frac{x}{1-x^3} \left(1 + (x-1) x^{3a_1} - (x-1) x^{3 a_2} + (x-1) x^{3 a_3} + \ldots \right)$$
Here we imagine that the $a_n$ are growing really really fast, (say $x = n!^3$ or something). Let $t = 1/3$ and let $\omega = e^{2 \pi i t}$. Then $c_V(1/3)$ exists if and only if
$$\lim_{x \rightarrow 1^{-}} \left(1 + (\omega x-1) x^{3a_1} - (\omega x-1) x^{3 a_2} + (\omega x-1) x^{3 a_3} + \ldots \right)$$
exists. The claim is that the value of this sum will vary between $1$ and $\omega$ as $x$ approaches $1$. The reason is that when $a_n$ increases sufficiently fast, one can choose $x$ so that the first $n$ values of $x^{3 a_n}$ are very close to one and the others decrease super-exponentially. For example, choose a value of $x$ such that
$$1 - \frac{1}{n^2 a_n} < x < 1 - \frac{n}{a_{n+1}} < 1 - \frac{n^{k}}{a_{n+k}}.$$
Then for $m \le n$,
$$x^{3 a_m} \ge x^{3 a_n} \ge \left(1 - \frac{1}{n^2 a_n}\right)^{3 a_n} \sim 1 + O\left(\frac{1}{n^2}\right),$$
but for $m = n+k > n$,
$$x^{3 a_m} \le \left(1 - \frac{3 n^{k}}{a_{n+k}}\right)^{3 a_{n+k}} \sim e^{-3 n^k}$$
which are thus negligible. Thus the value of the sum is approximately
$$1 + (\omega - 1) - (\omega -1) + (\omega -1) - \ldots $$
and so will vary between $1$ and $\omega$ and not converge.
On the other hand, it's not so hard to see that there will be no convergence issues for all other rational numbers except for $t = 1/3$ and $t = 2/3$. The point is that the $1-x$ term kills of any power series with terms that are so spread out.
But now consider:
$$\varsigma_W(x):=\varsigma_V(x) + \varsigma_V(x^3) + \varsigma_V(x^9) + \ldots $$
Since $V$ consisted only of terms $1$ and $2$ modulo $3$, there exists a corresponding set $W$. But now the analysis above shows that $T_W$ avoids the set $a/3^n$ for every $(a,3) = 1$ and $n \ge 1$. The point is that these values are screwed up exactly by $\sigma_V(x^{3^{n-1}})$. So $T_W$ is missing infinitely many rational numbers.
Some preliminary analysis suggests one can construct functions missing all rational numbers except zero but the analysis could be a bit painful. (The idea is similar to the above, but to slowly increase the modulus one is considering from $2!$ to $3!$ to $4!$ etc.)
Other observations:
Certainly a $V$ chosen randomly will, almost always (in the appropriate sense) have $S \subset T_V$ for any finite set $S$. I might even guess that for a random $V$ the set $T_V$ contains all rational numbers, did you try to prove that?