Time dependent delay differential equation

A faster and more straightforward approach is to use NDSolve as follows. Begin by noting that the first segment of the solution can be computed by

xd[t_?NumericQ] := 1.5;
s1 = NDSolve[[{x'[t] == x[t] (1 - xd[t]), x[0] == 1.5}, x[t], {t, 0, t1] // Values;

where t1 - (2 + Sin[t1]) == 0. With s1 determined, it becomes possible to compute the next section with

xd[t_?NumericQ] := s1[[0]][t - (2 + Sin[t])]

and integrating from t1 to t2, where t2 - (2 + Sin[t2]) == t1. In all, 109 steps are required to reach t = 200, calculated by

step = Rest@NestList[t /. FindRoot[t - (2 + Sin[t]) == #, {t, Max[#, 2]}] &, 0, 109]
(* {2.5542, 3.88062, 4.89775, 7.89684, ..., 196.712, 198.321, 199.334, 202.268} *)

Of course, executing NDSolve 109 times is both slow and cumbersome, requiring that the 109 solution segments be spliced together. Using NDSolve Components, however, simplifies the computation greatly. It is initialized with

xd[t_?NumericQ] := 1.5;
ndss = First[NDSolve`ProcessEquations[{x'[t] == x[t] (1 - xd[t]), x[0] == 1.5}, x[t], t]];
NDSolve`Iterate[ndss, step[[1]]];
s = First@NDSolve`ProcessSolutions[ndss] // Values;
xd[t_?NumericQ] := s[[0]][t - (2 + Sin[t])]

and completed by iterating through the remaining values of step

Do[NDSolve`Iterate[ndss, step[[i]]];
   s = First@NDSolve`ProcessSolutions[ndss] // Values;, {i, 2, 109}]

The iteration requires about 1/40 the time and 1/50 the memory of the approach used to obtain the second plot in my earlier answer. Plotting the final expression for s in Red and overlaying it on that second plot yields.

Agreement is excellent except at the tips of some of the highest peaks in the curve. Perhaps, this is due to the higher order interpolation used by NDSolve. (The earlier answer employs linear interpolation.) The key observation is that the two solutions do not drift apart as t increases.

The method described here should generalize to most ODEs with time-varying delays, provided that the minimum size of step elements is not too small.

Edited for clarity and improved accuracy.

Here is a simple solution, which perhaps can be improved. First, replace x[t] by Exp[y[t]], to obtain

y'[t] == 1 - Exp[y[t - d[t]]]

which guarantees that x[t] > 0 after discretization, and also is a bit simpler to discretize. Then, the natural discretization would be

f = 1/2 + (dl + a*Sin[(n - 1/2) dt])/dt
y[n] = y[n - 1] + (1 - Exp[y[n - f]]) dt

except that f is not an integer. Hence, interpolation is needed. For instance,

Clear[y]; dl = 2.; tl = 60; dt = 1/400; y0 = Log[1.5]; a = 0.;
Table[y[n] = y0, {n, -3/dt, 0}];
y[n_] := y[n] = (f = 1/2 + (dl + a*Sin[(n - 1/2) dt])/dt; y[n - 1] + 
    (1 - Exp[y[n - Floor[f]] (1 - Mod[f, 1]) + y[n - Ceiling[f]] Mod[f, 1]]) dt);
ListPlot[Table[Exp[y[n]], {n, 0, tl/dt}], PlotRange -> All, Joined -> True, 
    DataRange -> tl, ImageSize -> Large, LabelStyle -> {15, Bold, Black}]

which is the same result as

NDSolveValue[{x'[t] == x[t] (1 - x[t - 2]), x[t /; t <= 0] == 1.5}, x[t], {t, 0, 60}];
Plot[%, {t, 0, 60}, PlotRange -> All, ImageSize -> Large, LabelStyle -> {15, Bold, Black}]

as desired. Note that we chose a delay of 2.0 rather than 1.4 in the Question, because 2.0 is the average value of d. Now set

Clear[y]; dl = 2.; tl = 200; dt = 1/1000; y0 = Log[1.5]; a = 1.;

to reflect the sinusoidal variation in d. In addition, using a smaller time step is helpful for good accuracy, and a longer domain is desirable to show variations in the solution pattern. The result is

The irregularity is not necessarily surprising, and probably represents beating between the oscillations shown in the first plot and those in d. It also is possible that the solution is mildly chaotic.