10
$\begingroup$

So, it seems to me that the weights function in lm gives observations more weight the larger the associated observation's 'weight' value, while the lme function in lme does precisely the opposite. This can be verified with a simple simulation. 

#make 3 vectors- c is used as an uninformative random effect for the lme model
a<-c(1:10)
b<-c(2,4,6,8,10,100,14,16,18,20)
c<-c(1,1,1,1,1,1,1,1,1,1)

If you were now to run a model where you weight the observations based on the inverse of the dependent variable in lm, you can only generate the exact same result in nlme if you weight by just the dependent variable, without taking the inverse.

summary(lm(b~a,weights=1/b))
summary(lme(b~a,random=~1|c,weights=~b))

You can flip this and see the converse is true- specifying weights=b in lm requires weights=1/b to get a matching lme result.

So, I understand this much, I just want validation on one thing and to ask a question about another.

  1. If I want to weight my data based on the inverse of the dependent variable, is it fine to just code weights=~(dependent variable) within lme?
  2. Why is lme written to handle weights completely differently than lm? What is the purpose of this other than to generate confusion?

Any insight would be appreciated!

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ The answer to 2. is that they were written by very different people to do very different things. lm() needed to compatible with S and various books, nlme didn't, and it aimed to be more flexible, allowing the heterogeneity to be modelled more flexibly than lm allows. $\endgroup$
    – Gavin Simpson
    Commented Sep 18, 2013 at 22:20

1 Answer 1

Reset to default
15
$\begingroup$

Q1

In lme the notation weights = ~ b would result in the varFixed variance function being used with sole argument b. This function would add to the model a variance function $s^2(v)$ that has the form $s^2(v) = |v|$, where $v$ takes the values of the vector argument b.

Hence, you should use weights = ~ I(1/b) in lme() to have the variance of $\varepsilon_i = 1/b_i$.

In lm what you pass weights seems is the exact opposite; weights is inversely proportional to the variance.

I'm not 100% sure what you mean by weight my data, but if you mean provide the heterogeneous variance of the observations, then I think you want weights = ~ I(1/b).

Q2

My gut feeling (you'd have to ask the respective authors of the two functions) is that this is beacuse lm() and lme() were written by very different people to do very different things. lm() needed (was desired to be) to be compatible with S and various books, nlme didn't, and it aimed to be more flexible, allowing the heterogeneity to be modelled more flexibly than lm allows through the use of variance functions via the varFunc infrastructure.

Improve this answer
$\endgroup$
4
  • $\begingroup$ This is clear enough. By 'weight my data' I mean I want the model fit to consider that that large residuals should be expected from large observations, and fit something akin to percentage least squares, rather than ordinary least squares. ALSO- I deleted the cross post on stack overflow, sorry! $\endgroup$
    – colin
    Commented Sep 19, 2013 at 17:08
  • $\begingroup$ You might want to look at other variance functions in nlme then. What you are doing is saying that the variances of your observations are exactly the (absolute) value of b. It would seem better to just say that the variance increased with b. varPower() for example would have the variance as $\hat{\sigma}^2 \times |b|^{2\delta}$ with $\delta$ estimated a model parameter. This is OK if b doesn't take 0 values. If it can take 0 values, then the varExp() function may be better, there the variance is $var(\varepsilon_i) = \hat{\sigma}^2 \times e^{2\delta \times b_i}$. $\endgroup$
    – Gavin Simpson
    Commented Sep 19, 2013 at 17:29
  • $\begingroup$ In lm(), note the wording that the variance is proportional to the inverse of weights. In the lme code we discussed, b is the variance. Following your explanation, I don't think you really want that... Also note that if the variance increases with mean response, then a GLMM may be appropriate and the lme4 package would be suitable as it can model the mean-variance relationship directly, rather than via modification to the covariance matrix - which is what the lme code is doing. $\endgroup$
    – Gavin Simpson
    Commented Sep 19, 2013 at 17:33
  • $\begingroup$ Finally, sorry if I sounded grumpy on Stack Overflow. It wasn't intentional. I just forgot that you can't vote to close as OT & migrate to Cross Validated. You have to leave a comment as to why but I'd already left the first comment. Don't choose a SE site for your question based on the number of eyes that will see it. Choose the most appropriate venue. There is nothing wrong with promoting your question on Cross Validated to get more eyes, you could even post the link to it in the R public Chat Room on Stack Overflow. Crossposting or posting OT questions dilutes the resource if we have too many, hence close votes etc. $\endgroup$
    – Gavin Simpson
    Commented Sep 19, 2013 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.