Amanda Li
2014-Jan-13 22:43 UTC
[R] likelihood ratio test for mean difference assuming unequal variance
Hello, I am so sorry, but I have been struggling with the code for the entire day. I have a very simple dataset that looks like this: response=c(45,47,24,35,47,56,29) sub=c("A","A","B","B","C","C","C"£© time=c(1,2,1,2,1,2,3) gdata=cbind(response,sub,time) Namely, for three subjects, each has 2 or 3 observations. Assuming that each subject has a different variance, I want to test whether the mean for the three subjects are equal. I tried: library(nlme) weight <- varIdent(form = ~ 1 | sub ) lme1 <- lme(response~sub, weights = weight, data = gdata) However, it shows error: Error in getGroups.data.frame(dataMix, groups) : invalid formula for groups Can anyone help me? Eventually want to try this: lme1 <- lme(response~sub, weights = weight, data = gdata) lme2 <- lme(response~1, weights = weight, data = gdata) lrtest(lme1,lme2) Am I on the right track? Or can anyone think of a better method? Thank you very much! Best regards, Amanda [[alternative HTML version deleted]]
Ben Bolker
2014-Jan-14 02:34 UTC
[R] likelihood ratio test for mean difference assuming unequal variance
Amanda Li <amandali <at> uchicago.edu> writes:> > Hello, > > I am so sorry, but I have been struggling with > the code for the entire day. > > I have a very simple dataset that looks like this: > response=c(45,47,24,35,47,56,29) > sub=c("A","A","B","B","C","C","C"?? > time=c(1,2,1,2,1,2,3) > gdata=cbind(response,sub,time) > > Namely, for three subjects, each has 2 or 3 observations. > Assuming that each subject has a different > variance, I want to test whether > the mean for the three subjects are equal. > > I tried: > library(nlme) > weight <- varIdent(form = ~ 1 | sub ) > lme1 <- lme(response~sub, weights = weight, data = gdata) > > However, it shows error: > Error in getGroups.data.frame(dataMix, groups) : > invalid formula for groups > > Can anyone help me? Eventually want to try this: > lme1 <- lme(response~sub, weights = weight, data = gdata) > lme2 <- lme(response~1, weights = weight, data = gdata) > lrtest(lme1,lme2) > > Am I on the right track? Or can anyone think of a better method? > > Thank you very much! > > Best regards, > AmandaA variety of issues: (1) you should use data.frame() rather than cbind() to combine your data: cbind() makes it into a matrix, which converts everything to character vectors. (2) you should use gls() rather than nlme() -- you don't have a random effect. (3) I really hope this is a toy example and that your real data are larger, because estimating means and variance for three groups from 7 data points is *extremely* optimistic. (If you look at intervals(lme1) you'll see that the confidence intervals are very wide ...) (4) I'm guessing lrtest() is from the lmtest package, but you didn't tell us that ... Ben Bolker