[R] Help with lme
R. Michael Weylandt
michael.weylandt at gmail.com
Fri Nov 9 18:53:00 CET 2012
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.
More information about the R-help
mailing list