Help with lme
On Fri, Nov 9, 2012 at 5:27 PM, Bert Gunter <gunter.berton at gene.com> wrote:
Well, you've posted to the wrong list! First off, you're almost always better off posting to an R SIG list when one exists in your area of concern,as it does here: r-sig-mixed-models. Second, this appears to be primarily a statistics question, and R-help is not a statistics help list (though, I admit, the intersection is nonempty, and that may be the case here). For statistics primarily questions,you should generally post to a statistics help list like stats.stackexchange.com. Otherwise, you're fine. :-) Cheers, Bert On Fri, Nov 9, 2012 at 9:15 AM, Annie Hoen <anniehoen at gmail.com> wrote:
Hi,
This is my first time posting to the list so please forgive any breeches
of etiquette!
I am new to mixed-effects modeling.
This is my dataset:
subject treatment day replicate outcome
1 1 1 1 0.0
1 1 4 1 0.0
1 1 8 1 14.5
1 1 8 2 15.4
2 1 2 1 0.0
2 1 4 1 0.0
2 1 7 1 12.1
2 1 7 2 11.9
3 1 2 1 0.0
3 1 4 1 0.0
3 1 7 1 0.0
4 1 2 1 4.2
4 1 2 2 5.0
4 1 4 1 8.5
4 1 4 2 10.0
4 1 6 1 16.4
4 1 6 2 18.1
5 1 2 1 0.0
5 1 4 1 0.0
5 1 7 1 0.0
6 2 2 1 0.0
6 2 4 1 9.1
6 2 4 2 9.7
6 2 7 1 12.6
6 2 7 2 10.3
7 2 1 1 3.3
7 2 1 2 4.8
7 2 4 1 6.2
7 2 4 2 6.4
7 2 7 1 12.9
7 2 7 2 13.1
8 2 2 1 0.0
8 2 4 1 0.0
8 2 8 1 0.0
9 2 2 1 2.7
9 2 2 2 3.2
9 2 4 1 5.6
9 2 4 2 5.4
9 2 8 1 14.9
9 2 8 2 14.8
10 2 1 1 0.0
10 2 4 1 10.7
10 2 4 2 11.0
10 2 7 1 13.7
10 2 7 2 12.9
11 2 1 1 0.0
11 2 4 1 0.0
11 2 7 1 0.0
12 2 1 1 0.0
12 2 4 1 0.0
12 2 7 1 0.0
It can be made using this:
subject=c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 6,
6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10,
10, 11, 11, 11, 12, 12, 12)
treatment=c(rep(1, 20), rep(2, 31))
day=c(1, 4, 8, 8, 2, 4, 7, 7, 2, 4, 7, 2, 2, 4, 4, 6, 6, 2, 4, 7, 2, 4, 4,
7, 7, 1, 1, 4, 4, 7, 7, 2, 4, 8, 2, 2, 4, 4, 8, 8, 1, 4, 4, 7, 7, 1, 4, 7,
1, 4, 7)
replicate=c(1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 1, 1, 1,
1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1,
1, 1, 1, 1, 1)
outcome=c(0, 0, 14.5, 15.4, 0, 0, 12.1, 11.9, 0, 0, 0, 4.2, 5.0, 8.5,
10.0, 16.4, 18.1, 0, 0, 0, 0, 9.1, 9.7, 12.6, 10.3, 3.3, 4.8, 6.2, 6.4,
12.9, 13.1, 0,0,0, 2.7,3.2, 5.6, 5.4, 14.9, 14.8, 0, 10.7, 11.0, 13.7,
12.9, 0, 0, 0, 0, 0, 0)
data<-data.frame(cbind(subject, treatment, day, replicate, outcome))
Seconding Bert's advice to repost to R-SIG-Mixed-models (it's a great mailing list with brilliant statisticians) and commending a really good reproducible example. The only thing I'd note is that you generally don't want to use constructs like: data.frame(cbind(....)) cbind() will create a matrix object before data.frame() is called. The disadvantage of that is that a matrix can have only one sort of data, so cbind() will force it all to be the same (usually either numeric or character). You can loose information here, but, more importantly data.frame() doesn't know what went into cbind() so it doesn't bother to convert back to the original data type. Then you're left with a data.frame of all one data type, which more or less defeats the point of data.frames in the first place. :-) It's better to just call data.frame() directly with something like data.frame(subject, treatment, day, replicate, outcome) You'll also get the added bonus of R figuring out names for columns here. Of course, in your case, it doesn't hurt because your data really is all of the same type. But good practice and all that! Cheers and welcome! Michael
I have two groups of subjects, each given a different treatment. The outcome (growth) was observed post-treatment on each of three days. If there was growth, it was measured in duplicate measurements. There are uneven numbers of subjects in my two treatment groups. Also, the outcome was always observed on 3 days, but the exact day of observation is not always consistent. I just want to know the effect of treatment on outcome. This is the model I've run: model <- lme(outcome ~ treatment * day, random = list(subject = pdDiag(~ day)), data = data) summary(model) Can any experts out there let me know if I'm doing this right? Thanks!!!
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
-- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.