Hi Stata experts!

I have a situation that I've been unable to Google my way to a resolution.

I have a three period interrupted time series study and would like to find the change in level between the first period and third period using a mixed model with restricted maximum likelihood method and lincom. I can do this part of things without any hassle at all.

My dilemma comes when we want to use the small sample approximation, e.g. Satterthwaite. Again, doing things in a straightforward manner is no problem.

In some instances, the Satterthwaite degrees of freedom estimations drop quite low. We would like to set a lower-bound for the degrees of freedom.

In the example below, (made up data, but sort of represents what I'd like) we have degrees of freedom that dip below 2 for some parameters. Let's just say for the sake of the example that I would like to set a lower limit to the degrees of freedom to 2.

Currently the dofs are
| observa~n | lnsig_e | r_atr1
| level_0 slope_0 level_1 slope_1 level_2 slope_2 | _cons | _cons
-------------+------------------------------------------------------------------+-----------+-----------
r1 | 1.731149 1.836793 2.258847 1.567554 2.513234 2.298275 | . | .


And I would like them to change to 2, 2, 2.258847, 2, 2.513234, 2.298275. i.e. leave them alone if greater than 2, change to 2 if less than 2.

In practice sometimes the values are VERY different, I understand in this particular example there's not a huge difference between them, but in practice there can be ranges from below 1 to much, much higher with some of the datasets I'm working with!

When I use lincom I can set all degrees of freedom to a particular number OR use the small option, which passes the degrees of freedom from the mixed Satterthwaite model (or so I understand it to work?). I guess what I'm wondering is if I can set the dof to a particular matrix of my choosing?

Cheers for any help you can provide me with!
Below is the code I've used to generate the example problem (the graph shows what I'm looking at, this example is after the difference between the red line and green line at time point 61 (start of the green line).

Thank you kindly,
Simon.

Code:
clear 
graph drop _all 

set seed 1

set obs 90

// set parameters for level and slope changes
local level_0 0
local slope_0 0
local level_1 0
local slope_1 .1
local level_2 -3
local slope_2 -.1

// set autocorrelation parameter
local rho 0.3

gen time = _n

// create error
gen error = rnormal(0,1)
// autocorrelated error
replace error = error + `rho'*error[_n-1] in 2/90

// indicator variables for segments
gen indicator_0 = 1

gen indicator_1 = 1
replace indicator_1 = 0 if time < 31

gen indicator_2 = 1
replace indicator_2 = 0 if time < 61

// segments
gen level_0 = indicator_0
gen slope_0 = indicator_0 * time
gen level_1 = indicator_1
gen slope_1 = indicator_1 * (time-30-1)
gen level_2 = indicator_2
gen slope_2 = indicator_2 * (time-60-1)

// combine to single observation
gen observation = `level_0'*indicator_0 + `slope_0'*time*indicator_0 + ///
                  `level_1'*indicator_1 + `slope_1'*time*indicator_1 + ///
                  `level_2'*indicator_2 + `slope_2'*time*indicator_2 + ///
                  error

// normal dof
mixed observation level_0 slope_0 level_1 slope_1 level_2 slope_2 , nocons  res(ar 1, t(time)) var reml iter(1000)
local normal_level_0_beta = _b[level_0]
local normal_level_0_se = _se[level_0]

// find the level change from the counterfactual (extrapolated from first time period)    
lincom level_0 + slope_0*`=slope_0[31]' + level_1 + slope_1*`=slope_1[61]' + level_2 + slope_2*`=slope_2[61]'

// analysis with adjusted dof
mixed observation level_0 slope_0 level_1 slope_1 level_2 slope_2 , nocons  res(ar 1, t(time)) var reml iter(1000) dfmethod(satt)
local ttest_level_0_beta = _b[level_0]
local ttest_level_0_se = _se[level_0]
matrix dofs = e(df)
matlist dofs
local ttest_level_0_dof = dofs[1,1]

// compare the two upper confidence limits
local normal_ciu = `normal_level_0_beta' - invnormal(0.025)*`normal_level_0_se' 
local ttest_ciu = `ttest_level_0_beta' + invttail(`ttest_level_0_dof',0.025)*`ttest_level_0_se' 
display "normal ciu: `normal_ciu' ttest ciu: `ttest_ciu'"


// find the level change from the counterfactual (extrapolated from first time period)    
lincom level_0 + slope_0*`=slope_0[31]' + level_1 + slope_1*`=slope_1[61]' + level_2 + slope_2*`=slope_2[61]', small

// graph
gen counterfactual = _b[level_0] + _b[slope_0]*time
twoway scatter observation time ///
    || line counterfactual time ///
    || lfit observation time if time > 60