[R] trouble with model.tables SE means

David Winsemius dwinsemius at comcast.net
Wed Dec 23 18:05:16 CET 2009 

On Dec 23, 2009, at 11:45 AM, Jon Prince wrote:

> David and Peter,
>
> Thanks so much for all of your help, I think I understand R much  
> better as a result. Omitting the Error() term in my aov does indeed  
> allow me to get SE means, so I guess that was the issue. I suppose I  
> can go back and calculate the SE values for each p*t entry (averaged  
> across subject) from Matlab.
>
> By the way, the reason I separated the p and t objects derived from  
> rate_data in my initial example was for readability.

No need; see below.

> If I do it as below then I get the same result.
>
> In any case, again many thanks for your help. Cheers and Happy  
> Holidays.
>
> Jon
>
>
> rate_data=read.table("/Users/jonprince/Desktop/attend_pitch_3p3t.txt")
>  rate_data$p=factor(rate_data$V3)
>  rate_data$t=factor(rate_data$V2)
>  rate_data$subj=factor(rate_data$V1)
>  rate_data$rate=rate_data$V4
>  fm=aov(rate_data$rate ~ rate_data$p*rate_data$t + Error(rate_data 
> $subj/(rate_data$p*rate_data$t)),rate_data)

Try instead:
rm(p); rm(rate); rm(t); rm(subject)
# my memory of what your free-standing variables were named

fm=aov( rate ~ p*t + Error( subj/(p*t) ), data= rate_data)

R regression functions generally let you specify the column names as  
long as you supply the data.frame to a "data=" argument.

>
>
> Peter Alspach wrote:
>>
>> Tena korua David and Jon
>>
>> Without an Error() in the model, you get standard errors for the  
>> effects
>> and standard errors of the difference for the means.  With fully
>> balanced data, as in the example, these are directly comparable  
>> (compare
>> model.tables(npk.aov, se=T)*sqrt(2) with model.tables(npk.aov,
>> type='means', se=T)$se).  For data with uneven replication this  
>> will not
>> be the case.
>>
>> The standard errors of the differences for the means is not yet
>> implemented for aovlist (which is returned with an Error() in the  
>> model
>> as in npk.aovE).  I imagine this is because of the issues that  
>> arise in
>> comparing means from different strata.
>>
>> HTH .....
>>
>> Peter Alspach
>>
>>
>>> -----Original Message-----
>>> From: r-help-bounces at r-project.org
>>> [mailto:r-help-bounces at r-project.org] On Behalf Of David Winsemius
>>> Sent: Wednesday, 23 December 2009 4:13 p.m.
>>> To: Jon Prince
>>> Cc: r-help at r-project.org
>>> Subject: Re: [R] trouble with model.tables SE means
>>>
>>>
>>> On Dec 22, 2009, at 8:52 PM, Jon Prince wrote:
>>>
>>>
>>>> David Winsemius wrote:
>>>>
>>>>> On Dec 22, 2009, at 5:19 PM, Jon Prince wrote:
>>>>>
>>>>>> David Winsemius wrote:
>>>>>>
>>>>>>> On Dec 22, 2009, at 4:22 PM, Jon Prince wrote:
>>>>>>>
>>>>>>>
>>>>>>>> Hi, I'm new to R, with some experience with Matlab and SPSS.
>>>>>>>> I've figured out how to run my repeated measures anova and am
>>>>>>>> getting the right numbers for my effects (comparing
>>>>>>>>
>>> with results
>>>
>>>>>>>> from other software), but am having trouble with the
>>>>>>>>
>>> model.tables
>>>
>>>>>>>> function. Specifically, using:
>>>>>>>>
>>>>>>>> prints the means, but then won't do the SE values,
>>>>>>>>
>>> instead giving:
>>>
>>>>>>>> Warning message:
>>>>>>>> In model.tables.aovlist(fm, "means", se = TRUE) :
>>>>>>>> SEs for type 'means' are not yet implemented"
>>>>>>>>
>>>>>>>> Asking for SEs for "effects" works fine, but that's not what I
>>>>>>>> want. I searched the help for this issue and one other
>>>>>>>>
>>> person has
>>>
>>>>>>>> had this problem last year
>>>>>>>>
>>>>>>>>
>>> (http://markmail.org/message/k5yxxqcfiihvzvtp?q=list:r-project+mod
>>>
>>>>>>>> el%2Etables ), but the person helping them was unable
>>>>>>>>
>>> to replicate
>>>
>>>>>>>> it, inferring that it was an out-of-date version. My version  
>>>>>>>> is:
>>>>>>>>
>>>>>>>> R version 2.10.1 (2009-12-14)
>>>>>>>>
>>>>>>>> I only downloaded it the other day, and therefore
>>>>>>>>
>>> cannot have an
>>>
>>>>>>>> outdated version. How can I fix this error and get my SE  
>>>>>>>> values?
>>>>>>>> Apologies if I have not provided sufficient information, and
>>>>>>>> thanks in advance for your help.
>>>>>>>>
>>>>>>> When I look at the output of the first model.tables call copied
>>>>>>> from the help page, I see a list element that holds "se" values.
>>>>>>> Try:
>>>>>>>
>>>>>>> model.tables(fm,"means",se=TRUE)$se
>>>>>>>
>>>>>>>
>>>>>> Thanks for the rapid reply! Unfortunately adding the $se returns
>>>>>> NULL, and repeats the same warning message ("...not yet
>>>>>> implemented"). If you're not experiencing the issue, is
>>>>>>
>>> it possible
>>>
>>>>>> for me to replace the relevant code/source file with what
>>>>>>
>>> you have
>>>
>>>>>> (or would that require recompiling)? Could this be an OS
>>>>>>
>>> issue? I'm
>>>
>>>>>> running Mac OSX 10.6.2.
>>>>>>
>>>>>> By the way, I "replied all" on this message, but let me
>>>>>>
>>> know if that
>>>
>>>>>> is not the preferred convention. Cheers,
>>>>>>
>>>>> Reply all. That way people can correct my mistakes and general
>>>>> cluelessness. I'm running MacOSX 10.5.8 so it would seem
>>>>>
>>> less likely
>>>
>>>>> that is the explanation.
>>>>>
>>>>> 1) Did you run the example in the help pages?
>>>>>
>>>>> 2) When I look at :
>>>>>
>>>>>
>>>>>> methods(model.tables)
>>>>>>
>>>>> [1] model.tables.aov*     model.tables.aovlist*
>>>>>
>>>>> ... I see both an "aov" method and an "aovlist" method. Is it
>>>>> possible that there is something about the object that you are
>>>>> working on that makes it an aovlist at thus invokes a different
>>>>> function than what the help page invokes?
>>>>>
>>>>> The code would not require complination... it's available with
>>>>> getAnywhere() and I do not see any calls to compiled or .Internal
>>>>> subroutines. Tell me what happens with the above questions first.
>>>>>
>>>>>
>>>> Sorry for the delayed response, I've been trying to work out the
>>>> discrepancy between my data and the example data in terms
>>>>
>>> of getting
>>>
>>>> the SE means.
>>>>
>>>> 1.) Upon trying the example, it appears to work on my
>>>>
>>> machine (i.e., I
>>>
>>>> get the SE means). Since then I've been spending my time trying to
>>>> figure out why it doesn't do the same for my data. I can't
>>>>
>>> figure it
>>>
>>>> out.
>>>>
>>>> 2.) Given my lack of experience, it is quite possible that
>>>>
>>> I've loused
>>>
>>>> up an object somewhere, but I don't know how to track that down. I
>>>> have put my code below, and can send the data file too if you're
>>>> willing to take a look. Both the p and t factors have three levels,
>>>> combined factorially. There are four replications of each p- t
>>>> combination, and the rating data vary from 1 to 7.
>>>>
>>>> ##load data
>>>>   rate_data<-read.table("/Users/jonprince/Desktop/
>>>> attend_pitch_3p3t.txt")
>>>> ##set factors
>>>>   p<-factor(rate_data$V3)
>>>>   t<-factor(rate_data$V2)
>>>>   subject<-factor(rate_data$V1)
>>>> ##set data
>>>>   rate<-rate_data$V4
>>>> ##run anova
>>>>   fm<-aov(rate ~ p*t + Error(subject/(p*t)),rate_data)
>>>>
>>> But, but, but, ... none of those objects (p,t,subject, rate)
>>> are part of rate_data!
>>>
>>>
>>>> ##get summary
>>>>   summary(fm)
>>>> ##tables
>>>>   model.tables(fm,"means",se=TRUE)
>>>>
>>> Without more information, I would be flailing around, and
>>> even then I am not particularly experienced with the aov
>>> framework. My guess is that you will get better advice on the
>>> R-mixed-models-SIG. I would offer more information to that
>>> group, at a minimum the output of str on rate_data (once you
>>> actually create the variables IN rate_data.)
>>>
>>> --
>>> David
>>>
>>>> Thanks again, I really appreciate the help.
>>>>
>>>> Jon
>>>>
>>>>
>>>> --
>>>> Jon Prince
>>>> Postdoctoral Research Associate
>>>>
>>>>
>>>>
>>> David Winsemius, MD
>>> Heritage Laboratories
>>> West Hartford, CT
>>>
>>> ______________________________________________
>>> 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.
>>>
>>>
>>
>>
>>
>
> -- 
> Jon Prince
> Postdoctoral Research Associate
> Department of Psychology
> 356 Park Hall
> University at Buffalo, SUNY
> Buffalo, NY 14260
> office (716) 645 0235
> lab (716) 645 0225
> home (716) 839 1315
> jonprinc at buffalo.edu
> http://www.acsu.buffalo.edu/~jonprinc/

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

More information about the R-help mailing list