# Finding a trend in variance

I have written a little function that samples from a normal distribution with linearly changing variance

rnormltv <- function(n,mu,rsd) {
vc <- vector(mode="numeric", length=n)
i <- 1
for (x in seq(from=rsd[1], to=rsd[2], length=n)) {
vc[i]=rnorm(1,mu,x)
i <- i+1
}
return(vc)
}


Suppose we obtain 1000 samples, with a decreasing variance from 9 to 1.

set.seed(1)
x <- rnormltv(1000,0,c(9,1))
plot(x)


The goal is to estimate the parameters 9 and 1 given this data, the mean of zero, assuming a normal distribution and a linear trend in variance. Conceptually I can see how this can be done with maximum likelihood, but how can it by done in practice, preferably in R ?

Note: This is not homework.

-
If you can see how this is done with ML estimation, then just write a function that takes in the data and returns the ML estimates. –  Max Jan 17 '12 at 14:55
Thanks, but I don't know how to do that in R. Wouldn't the likelihood function would be unwieldy with n=1000 ? I fiddled around a little in Maple and didn't have much success there either. –  P Sellaz Jan 17 '12 at 15:10
The likelihood is a function of just two parameters: what's "unwieldy" about that? The likelihood equations for a linear trend with intercept $a$ and slope $b$ are (1) $\sum_t\left(1/(a+bt) - x_t^2/(a+bt)^3\right)=0$ and (2) $\sum_t\left(t/(a+bt) - t x_t^2/(a+bt)^3\right)=0$. Since this is not homework, and presumably addresses real data, you should be very cautious about these strong assumptions. Before proceeding, consider an exploratory analysis that includes robust checking of the linearity-of-variance and zero-mean assumptions. –  whuber Jan 17 '12 at 15:24
Thanks for your comment. It seems to be trivial to you, but it's unwieldy to me and I don't know how to solve such equations in R. –  P Sellaz Jan 17 '12 at 15:41
Use any of the many solvers available or, for greater convenience, use one of the dedicated ML "wrappers," such as MaxLik (which I haven't tried but looks like just the thing). The only code you have to write is the expression of the log likelihood function. I did not mean to imply the problem is trivial; my point is that having 1000 (or 1,000,000 for that matter) data values is not a complication in ML estimation. –  whuber Jan 17 '12 at 15:50

Here is a similar solution in R.

Please note that the code you provided generates random normal variables with a linear trend in sthe standard deviation, not in the variance. So I used this, you can easily modify the code for other kind of trends.

The only caveat is to compute directly the log likelihood, which will behave better with big sample, where the likelihood is very near to 0

set.seed(1)
x <- rnormltv(1000,0,c(9,1))
f <- function(a,b,x) sum( dnorm(x, mean=0, sd = a+(b-a)*seq(0,1,length=length(x)), log=TRUE) )

nlm( function(p) -f(p[1], p[2], x), p = c(8,1) )
$minimum [1] 2923.871$estimate
[1] 9.120767 1.126542

$gradient [1] 4.048507e-05 2.179800e-05$code
[1] 1

$iterations [1] 9  So the estimate is$a = 9.12$and$b=1.13$. Please follow whuber advice to check the robustness of this estimation to your assumptions (zero mean, linear trend in standard deviation). If you need this for a particular kind of data analysis, always have a look on what other people did with the same kind of data. - Thank you. I'll mark this as the accepted answer soon. First, if you don't mind I'd like to ask a related question which refers to the points @whuber and yourself made. I asked about this because it arose as a consequence of other modelling I am doing. The /real/ problem is described in this open question: stats.stackexchange.com/questions/19911/… It is thought that there has been a near-linear change in randoms effects variance, and that is actually what I'm eventually hoping to model. Cont...... – P Sellaz Jan 18 '12 at 12:36 ...Cont... At this point I am searching for published papers, in any field/discipline, where a trend in variance has been modeled at level 1 in a multilevel model. So far I have not found any, but I am pretty sure that there should be some. If anyone can point me in the right direction, or provide links to published research I would be very grateful. – P Sellaz Jan 18 '12 at 12:38 I know some papers on ad hoc "variance stabilizing transformations" applied to variables where the variance depends on the mean value, but that’s not the same. In you hospital data, which variable seems to have such a trend? If you want to test the existence of such a trend, you don’t really need to model it, you’d better use a test that detects non-constant variance. "Test for homoscedasticity" should be the keyword. – Elvis Jan 18 '12 at 13:45 Hi again. Sorry that I'm not explaining it very well. It's not that the variance of any particular variable has a trend. It's the variance of the hospital-level effects (random effects) that has a trend. It is hypothesized that the variation in hospital-level effects has decreased linearly over the period of the data, due to more /consistency/ in clinical diagnosis and treatment (not necessarily /better/ treatment/outcomes, though that is also likely, but not the point). – P Sellaz Jan 18 '12 at 14:54 I think I got it. I think you should ask a specific question on comparison of estimated variances of random effects, explaining how you conducted the analysis (describe variables, models used for analysis, if possible R commands used...). I don’t even know how to perform a logistic regression with random effects, so I can’t help more! – Elvis Jan 18 '12 at 16:45 show 1 more comment You mentioned initially trying Maple. I did, too, and here's what I got (disclosure: I work for them): # To generate a test sample: decayingSample := proc(N :: nonnegint, a :: positive, b :: positive,$)
local i, sd, sample, result;
result := Vector(N, 'datatype=float');
sample := Statistics:-Sample(Normal(0, sd));
for i to N do
sd := a + (i-1)*(b-a)/(N-1);
result[i] := sample(1)[1];
end do;
return result;
end proc:

# aa, bb are the values for our test sample; a and b will represent
# the unknown symbolic parameters to be determined.
NN, aa, bb := 1000, 9, 1:
ds := decayingSample(NN, aa, bb):
with(Statistics):
PointPlot(ds);

# Form the symbolic log likelihood for the ith entry.
ll := LogLikelihood(Normal(0, normal(a + (i-1)*(b-a)/(NN-1))),
x, samplesize=1) assuming 'real';
# Now add these guys for i = 1 .. NN:
expr := evalf(add('eval(ll, ['i' = i, x[1] = ds[i]])', i = 1..NN)):
# Create an initial guess for a and b:
initialGuess := [a = StandardDeviation(ds[1 .. iquo(NN, 2)]),
b = StandardDeviation(ds[iquo(NN, 2) .. NN])];
# Optionally:
# infolevel[Optimization] := 5: # for more detailed info
Optimization:-Maximize(expr, {a >= 0, b >= 0}, 'initialpoint' = initialGuess);


This lead to an estimate of $a = 8.72, b = 1.07$; trying it on a second sample gave $a = 9.20, b = 0.94$.

-
Hi Erik and thanks. You know, I almost posted the question on mapleprimes, not SE ! I ran your code and also got good results, albeit with a "Warning, undefined value encountered" –  P Sellaz Jan 17 '12 at 19:24
Yeah, such messages just mean that one of the points that Optimization evaluates expr at leads to an undefined value. You can ignore them. –  Erik P. Jan 17 '12 at 19:32