[R] ANOVA

[Angelo Scozzarella]

Hi,

How can I make an  ANOVA if I haven't got all data set but I  know the  
numbers of subjects for each group, the mean and di standard deviation  
for each group?


Thanks


Angelo Scozzarella

__
R-help@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.

[R] Anova

[Dani Valverde]

Hello,
I am trying to perform an ANOVA on a dataframe like this:

 ID Grup  ML  PC  
mI   Gly Glx
X373.txt 11 0.004685889  0.05994939 0.004589104   0.007144027  
0.02042549
X373bis.txt  21 0.004685889   0.0599494 0.004589104   0.007144028   
0.0204255
X376.txt 31 0.007604898  0.07600049  0.019352470.01678972  
0.03363403
X392.txt 41 0.004982747  0.02760991  0.018555960.01413439  
0.02013267
X483.txt 51 0.009131522  0.02636227  0.024852320.01493074  
0.02532194
X520.txt 61 0.008553895  0.06304886  0.024902610.01907819  
0.03356447
X671.txt 71 0.008004383  0.05243909   0.0229447  0.024266  
0.03553224
X673A.txt81  0.01306893  0.02799318  0.01155338   0.007668364  
0.02123958
X673B.txt91 0.006099385   0.0516078  0.030993240.02664487  
0.03766267
X674.txt101 0.005476397  0.04945545  0.014797420.01084873  
0.02307815
X681.txt111 0.009807893  0.01680487 0.004524913   0.004705414 
0.008969453
X725.txt121  0.01018343   0.0119391 0.007911063   0.008005791  
0.01014625
X746.txt131 0.008864319  0.02484126  0.011980480.01276371  
0.01042683
X772.txt141 0.005966057  0.07398342   0.02103990.01606057  
0.02932874
X814.txt151  0.01350947   0.0317243 0.003234958   0.002958675 
0.005194511
X835.txt161 0.005509359  0.02670237  0.074253230.04095533
0.031017
X847.txt171 0.004784228  0.05349206  0.011836790.01230316  
0.03073348
X851.txt181 0.004395889  0.03812635  0.019986530.02165763  
0.02074277
X865.txt191 0.001514055  0.07128447   0.01068870.00802192  
0.01478855
X878.txt201 0.003693918  0.03864821  0.015717970.01972568
0.022711
et2029.txt  211  0.05204285  0.03769878  0.031687990.03631842  
0.02732589
et2348.txt  221  0.01407517  0.01679014 0.004900163   0.004243802  
0.01052361


The table is truncated because it's very big, but it has 3 different 
groups (Grup variable) and it's called GBM.

I try to use the command:

aov(GBM$Grup~GBM$ML,data=GBM)

But I get this error:

/Error in storage.mode(y) <- "double" :
 invalid to change the storage mode of a factor
In addition: Warning message:
In model.response(mf, "numeric") :
 using type="numeric" with a factor response will be ignored/

How can I solve this?
Best,

Dani

--
Daniel Valverde Saubí

Grup de Biologia Molecular de Llevats
Facultat de Veterinària de la Universitat Autònoma de Barcelona
Edifici V, Campus UAB
08193 Cerdanyola del Vallès- SPAIN

Centro de Investigación Biomédica en Red
en Bioingeniería, Biomateriales y
Nanomedicina (CIBER-BBN)

Grup d'Aplicacions Biomèdiques de la RMN
Facultat de Biociències
Universitat Autònoma de Barcelona
Edifici Cs, Campus UAB
08193 Cerdanyola del Vallès- SPAIN
+34 93 5814126

__
R-help@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.

[R] Anova

[daniel jupiter]

Hi all,

I apologize for what might be a silly question.

I am interested in doing a one way anova.
This is not too hard in and of itself, either with anova, aov or oneway.test
.

However, I need to
1) get pvalues,
2) do a posthoc analysis with Tukey HSD,
3) and have (sometimes) an unbalanced design.

I just can't seem to put all the pieces together.

Any suggestions?

Thanks in advance,

Dan.

-- 
Daniel C. Jupiter, Ph.D.
Postdoctoral Research Associate
Department of Systems Biology and Translational Medicine
College of Medicine
Texas A&M Health Science Center
702 SW H.K. Dodgen Loop
Temple, TX 76504

979.997.2106 | Fax 254.742.7145
[EMAIL PROTECTED] | www.tamhsc.edu

[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA

[kayj]

Hi All,


I have about one hundred patients and all the patients had their glucose
measured on three different days. The days are all the same for all he
patients. So I have three measurement for each patient . I want to know
whether the day when the glucose was measured has an effect on the
measurements. I was thinking to use a single factor analysis of variance but
I am not sure how to do it in R. any other suggestion on dealing with this
problem is welcome.

Thanks,

-- 
View this message in context: 
http://www.nabble.com/ANOVA-tp22353919p22353919.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.

[R] Anova

[stephen sefick]

melt.updn <- structure(list(date = structure(c(11808, 11869, 11961, 11992,
12084, 12173, 12265, 12418, 12600, 12631, 12753, 12996, 13057,
13149, 11808, 11869, 11961, 11992, 12084, 12173, 12265, 12418,
12600, 12631, 12753, 12996, 13057, 13149), class = "Date"), variable =
structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("unrestored",
"restored"), class = "factor"), value = c(1.34057641541824, 0.918021774919366,
0.905654270934854, 0.305945104043220, 0.58298856330543, 1.36580645291274,
0.874195629894938, 0.87482377014642, 0.930267689669002, 0.41753134369356,
1.09248531450337, 1.72571397293738, 0.305751868168171, 0.584498524462223,
0.983300317501076, 1.27216569968585, 0.730578393573363, 0.88361473836175,
1.16501295544266, 2.08896500025784, 0.664286881841064, 1.03859387871079,
1.39172581649833, 0.323405269371357, 1.00207568577518, 1.54383416626015,
0.611261918697393, 0.848992483196744)), .Names = c("date", "variable",
"value"), row.names = c(NA, -28L), class = "data.frame")

aov(value~variable, data=melt.updn)

I am having problems making sure that I am doing the correct analysis.
 I am trying to see if there is a difference in the mean of the
restored segment versus the unrestored segment (variable in x).  These
are repeated measures on the same treatments through time.  Is there a
way to control for the differences in time steps?  Any ideas?
thanks for the help,



-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@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.

[R] ANOVA

[Georgina Sarah Humphreys]

I have the attached data set (csv) and I want to run an analysis of variance on 
the wingsize data (comparing infected vs non-infected) within and between 
experiments.

Can anyone help me with the command I should use?

Many thanks
Georgina



PhD Student
Division of Infection and Immunity
B5-29, GBRC
120 University Place
Glasgow
G12 8TA
Tel: 0141 330 5650
__
R-help@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.

[R] ANOVA error

[Dani Valverde]

Hello,
I have a dataframe named myMatrix with the structure

TreatmentTimeCrmIb ...

Being the treatment and time the predictors and Cr, mIb and so on the 
response variables. When I call


Cr.aov <- aov(Cr~Treatment, data=myMatrix)

I got this error:

Error in storage.mode(y) <- "double" :
 invalid to change the storage mode of a factor
In addition: Warning message:
In model.response(mf, "numeric") :
 using type="numeric" with a factor response will be ignored

Can anyone help me? I take the chance to ask another question. I would 
like to perform an ANOVA on different variables. Is the call 
aov(Cr+mIb+...~Treatment, data=myMatrix) o correct way of doing it? What 
I would like to have is a p value for the variance for each response 
variable.

Best,

Dani

--
Daniel Valverde Saubí

Grup de Biologia Molecular de Llevats
Facultat de Veterinària de la Universitat Autònoma de Barcelona
Edifici V, Campus UAB
08193 Cerdanyola del Vallès- SPAIN

Centro de Investigación Biomédica en Red
en Bioingeniería, Biomateriales y
Nanomedicina (CIBER-BBN)

Grup d'Aplicacions Biomèdiques de la RMN
Facultat de Biociències
Universitat Autònoma de Barcelona
Edifici Cs, Campus UAB
08193 Cerdanyola del Vallès- SPAIN
+34 93 5814126

__
R-help@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.

Re: [R] ANOVA

[stephen sefick]

I don't know of a way to do it without the data

On Sun, Aug 10, 2008 at 11:35 AM, Angelo Scozzarella
<[EMAIL PROTECTED]> wrote:
> Hi,
>
> How can I make an  ANOVA if I haven't got all data set but I  know the
> numbers of subjects for each group, the mean and di standard deviation for
> each group?
>
> Thanks
>
>
> Angelo Scozzarella
>
> __
> R-help@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.
>



-- 
Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@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.

Re: [R] ANOVA

[Chuck Cleland]

On 8/10/2008 11:35 AM, Angelo Scozzarella wrote:

Hi,

How can I make an  ANOVA if I haven't got all data set but I  know the 
numbers of subjects for each group, the mean and di standard deviation 
for each group?


  See anova.mean() in the HH package:

http://finzi.psych.upenn.edu/R/library/HH/html/anova.mean.html


Thanks

Angelo Scozzarella

__
R-help@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. 


--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

__
R-help@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.

[R] ANOVA help

[Gareth Campbell]

Hi,

I'm doing anova on a matrix of multivariate data where I want to assess the
effect of each column (element).

My matrix is 86 rows x 31 columns.  I've created a grouping factor of length
86 containing group assignments of 6 types.

Then I run:

x<- aov(matrix~grouping.factor)
summary(aov.fit.raw, test="Wilks")

This is working fine enough, but I'm getting different results to someone
who I'm comparing with - am I doing what I think I am doing here??

What I think I'm doing is - take the first column (element) and then look at
that element between the 6 groups and report F, Pvalue etc..., then move
onto the 2nd column and repeat.

Thanks



-- 
Gareth Campbell
PhD Candidate
The University of Auckland

P +649 815 3670
M +6421 256 3511
E [EMAIL PROTECTED]
[EMAIL PROTECTED]

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova

[Prof Brian Ripley]

You are trying to do ANOVA on a categorical response.


How can I solve this?


Use a more appropriate model.  Did you mean ML to be the response and Grup 
the group rather than v.v.?  Then try


aov(ML ~ Grup, data = GBM)

On Tue, 12 Aug 2008, Dani Valverde wrote:


Hello,
I am trying to perform an ANOVA on a dataframe like this:

ID Grup  ML  PC  mI 
Gly Glx
X373.txt 11 0.004685889  0.05994939 0.004589104   0.007144027 
0.02042549
X373bis.txt  21 0.004685889   0.0599494 0.004589104   0.007144028 
0.0204255
X376.txt 31 0.007604898  0.07600049  0.019352470.01678972 
0.03363403
X392.txt 41 0.004982747  0.02760991  0.018555960.01413439 
0.02013267
X483.txt 51 0.009131522  0.02636227  0.024852320.01493074 
0.02532194
X520.txt 61 0.008553895  0.06304886  0.024902610.01907819 
0.03356447
X671.txt 71 0.008004383  0.05243909   0.0229447  0.024266 
0.03553224
X673A.txt81  0.01306893  0.02799318  0.01155338   0.007668364 
0.02123958
X673B.txt91 0.006099385   0.0516078  0.030993240.02664487 
0.03766267
X674.txt101 0.005476397  0.04945545  0.014797420.01084873 
0.02307815
X681.txt111 0.009807893  0.01680487 0.004524913   0.004705414 
0.008969453
X725.txt121  0.01018343   0.0119391 0.007911063   0.008005791 
0.01014625
X746.txt131 0.008864319  0.02484126  0.011980480.01276371 
0.01042683
X772.txt141 0.005966057  0.07398342   0.02103990.01606057 
0.02932874
X814.txt151  0.01350947   0.0317243 0.003234958   0.002958675 
0.005194511
X835.txt161 0.005509359  0.02670237  0.074253230.04095533 
0.031017
X847.txt171 0.004784228  0.05349206  0.011836790.01230316 
0.03073348
X851.txt181 0.004395889  0.03812635  0.019986530.02165763 
0.02074277
X865.txt191 0.001514055  0.07128447   0.01068870.00802192 
0.01478855
X878.txt201 0.003693918  0.03864821  0.015717970.01972568 
0.022711
et2029.txt  211  0.05204285  0.03769878  0.031687990.03631842 
0.02732589
et2348.txt  221  0.01407517  0.01679014 0.004900163   0.004243802 
0.01052361


The table is truncated because it's very big, but it has 3 different groups 
(Grup variable) and it's called GBM.

I try to use the command:

aov(GBM$Grup~GBM$ML,data=GBM)

But I get this error:

/Error in storage.mode(y) <- "double" :
invalid to change the storage mode of a factor
In addition: Warning message:
In model.response(mf, "numeric") :
using type="numeric" with a factor response will be ignored/

How can I solve this?
Best,

Dani

--
Daniel Valverde Saubí

Grup de Biologia Molecular de Llevats
Facultat de Veterinària de la Universitat Autònoma de Barcelona
Edifici V, Campus UAB
08193 Cerdanyola del Vallès- SPAIN

Centro de Investigación Biomédica en Red
en Bioingeniería, Biomateriales y
Nanomedicina (CIBER-BBN)

Grup d'Aplicacions Biomèdiques de la RMN
Facultat de Biociències
Universitat Autònoma de Barcelona
Edifici Cs, Campus UAB
08193 Cerdanyola del Vallès- SPAIN
+34 93 5814126

__
R-help@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.



--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595__
R-help@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.

Re: [R] Anova

[milton ruser]

Dear Daniel,

Are you really interested on the analyzis of the effect of ML on Grup, our
on the effect of Grup on ML?

I reproduce your sample database, but "changed" some Grup values from 1 to 2
to allow an analysis of variance.

See the example below.

Vale!

miltinho astronaura
Brazil



GBM<-read.table(stdin(),head=T, sep=",")
ID,Grup,ML,PC,mI,Gly,Glx
X373.txt,1,1,0.004685889,0.05994939,0.004589104,0.007144027,0.02042549
X373bis.txt,2,1,0.004685889,0.0599494,0.004589104,0.007144028,0.0204255
X376.txt,3,1,0.007604898,0.07600049,0.01935247,0.01678972,0.03363403
X392.txt,4,1,0.004982747,0.02760991,0.01855596,0.01413439,0.02013267
X483.txt,5,1,0.009131522,0.02636227,0.02485232,0.01493074,0.02532194
X520.txt,6,1,0.008553895,0.06304886,0.02490261,0.01907819,0.03356447
X671.txt,7,1,0.008004383,0.05243909,0.0229447,0.024266,0.03553224
X673A.txt,8,1,0.01306893,0.02799318,0.01155338,0.007668364,0.02123958
X673B.txt,9,1,0.006099385,0.0516078,0.03099324,0.02664487,0.03766267
X674.txt,10,1,0.005476397,0.04945545,0.01479742,0.01084873,0.02307815
X681.txt,11,1,0.009807893,0.01680487,0.004524913,0.004705414,0.008969453
X725.txt,12,1,0.01018343,0.0119391,0.007911063,0.008005791,0.01014625
X746.txt,13,1,0.008864319,0.02484126,0.01198048,0.01276371,0.01042683
X772.txt,14,1,0.005966057,0.07398342,0.0210399,0.01606057,0.02932874
X814.txt,15,1,0.01350947,0.0317243,0.003234958,0.002958675,0.005194511
X835.txt,16,1,0.005509359,0.02670237,0.07425323,0.04095533,0.031017
X847.txt,17,1,0.004784228,0.05349206,0.01183679,0.01230316,0.03073348
X851.txt,18,1,0.004395889,0.03812635,0.01998653,0.02165763,0.02074277
X865.txt,19,1,0.001514055,0.07128447,0.0106887,0.00802192,0.01478855
X878.txt,20,1,0.003693918,0.03864821,0.01571797,0.01972568,0.022711
et2029.txt,21,1,0.05204285,0.03769878,0.03168799,0.03631842,0.02732589
et2348.txt,22,1,0.01407517,0.01679014,0.004900163,0.004243802,0.01052361

GBM[11:22,"Grup"]<-2
GBM

aov(Grup~ML,data=GBM)
aov(ML~Grup,data=GBM)

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova

[Gavin Simpson]

On Sat, 2008-03-15 at 12:48 -0500, daniel jupiter wrote:
> Hi all,
> 
> I apologize for what might be a silly question.
> 
> I am interested in doing a one way anova.
> This is not too hard in and of itself, either with anova, aov or oneway.test
> .
> 
> However, I need to
> 1) get pvalues,

if obj is the result of anova, aov, oneway.test, then

str(obj) ## for anova
str(summary(obj)) ## for aov
str(obj) ## for oneway.test

to find the names of the elements of obj that contain the p-values of
the various tests/models. In the first two you are looking for the
component "Pr(>F)" and the latter is obvious ("p.value")

For summary(aov) objects, the result is a list so this gets the p-value
you need:

obj[[1]]$`Pr(>F)` or obj[[1]][,5]]

for anova then this:

obj$`Pr(>F)` or obj[,5]

note the quoting of the component name using backticks.

For oneway.test

obj$p.value

> 2) do a posthoc analysis with Tukey HSD,

?TukeyHSD for the results of aov

> 3) and have (sometimes) an unbalanced design.

See ?lme in package nlme

> 
> I just can't seem to put all the pieces together.
> 
> Any suggestions?

I'm not sure what the problem is here - you don't say. All of what I say
above is documented in the relevant help pages for the various functions
and using str() is a basic tenet of using R and looking at returned
objects.

Ok, you might have needed help with getting the p-values for some of
those tests/models, but 2) and 3) are answered on ?aov

For what you describe, stick with aov for balanced designs if you want
to do TukeyHSD as there is a method for aov objects (otherwise) you'll
need to refit the model.

For unbalanced designs, check out lme and for that you may need to
get/borrow the book by Pinhiero and Bates, reference details of which
are given in item [7] on:

http://www.r-project.org/doc/bib/R-books.html

> 
> Thanks in advance,
> 
> Dan.
> 

HTH

G
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

__
R-help@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.

Re: [R] Anova

[Douglas Bates]

On Sun, Mar 16, 2008 at 6:39 AM, Gavin Simpson <[EMAIL PROTECTED]> wrote:
> On Sat, 2008-03-15 at 12:48 -0500, daniel jupiter wrote:
>  > Hi all,
>  >
>  > I apologize for what might be a silly question.
>  >
>  > I am interested in doing a one way anova.
>  > This is not too hard in and of itself, either with anova, aov or 
> oneway.test
>  > .
>  >
>  > However, I need to
>  > 1) get pvalues,
>
>  if obj is the result of anova, aov, oneway.test, then
>
>  str(obj) ## for anova
>  str(summary(obj)) ## for aov
>  str(obj) ## for oneway.test
>
>  to find the names of the elements of obj that contain the p-values of
>  the various tests/models. In the first two you are looking for the
>  component "Pr(>F)" and the latter is obvious ("p.value")
>
>  For summary(aov) objects, the result is a list so this gets the p-value
>  you need:
>
>  obj[[1]]$`Pr(>F)` or obj[[1]][,5]]
>
>  for anova then this:
>
>  obj$`Pr(>F)` or obj[,5]
>
>  note the quoting of the component name using backticks.
>
>  For oneway.test
>
>  obj$p.value
>
>
>  > 2) do a posthoc analysis with Tukey HSD,
>
>  ?TukeyHSD for the results of aov
>
>
>  > 3) and have (sometimes) an unbalanced design.
>
>  See ?lme in package nlme
>
>
>  >
>  > I just can't seem to put all the pieces together.
>  >
>  > Any suggestions?
>
>  I'm not sure what the problem is here - you don't say. All of what I say
>  above is documented in the relevant help pages for the various functions
>  and using str() is a basic tenet of using R and looking at returned
>  objects.
>
>  Ok, you might have needed help with getting the p-values for some of
>  those tests/models, but 2) and 3) are answered on ?aov
>
>  For what you describe, stick with aov for balanced designs if you want
>  to do TukeyHSD as there is a method for aov objects (otherwise) you'll
>  need to refit the model.
>
>  For unbalanced designs, check out lme and for that you may need to
>  get/borrow the book by Pinhiero and Bates, reference details of which
>  are given in item [7] on:
>
>  http://www.r-project.org/doc/bib/R-books.html

Thanks for the plug, Gavin, but lme and my book with Jose are for
linear *mixed-effects* models.  I think the question here is about
unbalanced fixed-effects models and those can be fit quite simply by
lm, aov and friends.

Users, quite reasonably, expect that model fitting in R is done using
the formulas that they saw in an introductory textbook.  In fact, that
is almost never the case.  Functions like lm and aov use computational
methods that do not depend on balance in the data.


>
>  >
>  > Thanks in advance,
>  >
>  > Dan.
>  >
>
>  HTH
>
>  G
>  --
>  %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>   Dr. Gavin Simpson [t] +44 (0)20 7679 0522
>   ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
>   Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
>   Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
>   UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
>  %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
>
>
>  __
>  R-help@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.
>

__
R-help@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.

[R] anova help

[Joe King]

Hi all, I am trying to run a two factor anova, but one of the factors is a
random factor, now I am also running in SPSS and it seems its dividing by
the wrong term to get the appropriate F term. here is my data. In SPSS the F
scores about double the ones in R, how can I specify one of my factors as a
random factor or change it to where it does the right model fitting? I am
using the lm command instead of glm. I am new to R so this might seem basic.

 

Joe King, M.A.
  j...@joepking.com

"Never give in, never give in, never; never; never; never - in nothing,
great or small, large or petty - never give in except to convictions of
honor and good sense" - Winston Churchill

 

"You have enemies? Good. That means you've stood up for something, sometime
in your life." - Winston Churchill

 


[[alternative HTML version deleted]]

__
R-help@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.

[R] Anova Help?

[Dar]

I am conducting an experiment where students are put into 5 total
groups (one is the control group).  They are given a task and then I'm
measuring if there are differences (A 5 X 1 between-subject design
will be used for the experiment).  I'm a little confused on the data
explanation (or should I say how do I explain what is being analyzed
versus just comparing values)  Any help?

__
R-help@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.

Re: [R] ANOVA

[Kingsford Jones]

On Thu, Mar 5, 2009 at 8:31 AM, kayj  wrote:
>
> Hi All,
>
>
> I have about one hundred patients and all the patients had their glucose
> measured on three different days. The days are all the same for all he
> patients. So I have three measurement for each patient . I want to know
> whether the day when the glucose was measured has an effect on the
> measurements. I was thinking to use a single factor analysis of variance but
> I am not sure how to do it in R.

Two options for one-way ANOVA are aov(response ~ factor), or
lm(response ~ factor) followed by anova(lm.object).

Butthis would probably be a bad idea.  I recommend Pinheiro and
Bates (2000), where there are many examples of modeling effects within
subjects over time.

Here's an example with a continuous time effect

library(nlme)
f1 <- lme(distance~age, data=Orthodont, random= ~1 + age|Subject,
weights=varPower())
plot(f1, distance~fitted(.)|Subject, abline=c(0,1))
intervals(f1)


hth,

Kingsford Jones


>any other suggestion on dealing with this
> problem is welcome.
>
> Thanks,
>
> --
> View this message in context: 
> http://www.nabble.com/ANOVA-tp22353919p22353919.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@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.
>

__
R-help@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.

Re: [R] Anova

[Gavin Simpson]

On Wed, 2009-05-13 at 12:43 -0400, stephen sefick wrote:
> melt.updn <- structure(list(date = structure(c(11808, 11869, 11961, 11992,
> 12084, 12173, 12265, 12418, 12600, 12631, 12753, 12996, 13057,
> 13149, 11808, 11869, 11961, 11992, 12084, 12173, 12265, 12418,
> 12600, 12631, 12753, 12996, 13057, 13149), class = "Date"), variable =
> structure(c(1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("unrestored",
> "restored"), class = "factor"), value = c(1.34057641541824, 0.918021774919366,
> 0.905654270934854, 0.305945104043220, 0.58298856330543, 1.36580645291274,
> 0.874195629894938, 0.87482377014642, 0.930267689669002, 0.41753134369356,
> 1.09248531450337, 1.72571397293738, 0.305751868168171, 0.584498524462223,
> 0.983300317501076, 1.27216569968585, 0.730578393573363, 0.88361473836175,
> 1.16501295544266, 2.08896500025784, 0.664286881841064, 1.03859387871079,
> 1.39172581649833, 0.323405269371357, 1.00207568577518, 1.54383416626015,
> 0.611261918697393, 0.848992483196744)), .Names = c("date", "variable",
> "value"), row.names = c(NA, -28L), class = "data.frame")
> 
> aov(value~variable, data=melt.updn)

You can think of this as a linear model and just use lm:

lm(value~variable, data=melt.updn)

> 
> I am having problems making sure that I am doing the correct analysis.
>  I am trying to see if there is a difference in the mean of the
> restored segment versus the unrestored segment (variable in x).  These
> are repeated measures on the same treatments through time.  Is there a
> way to control for the differences in time steps?  Any ideas?
> thanks for the help,

One option is to fit this model using generalised least squares:

## do some plotting to look at potential differences:

require(lattice)
xyplot(value ~ time | variable, data = melt.updn, 
   type = c("p","smooth"))
## so perhaps some evidence of trend,
## different in the two groups possibly
bwplot(value ~ variable, data = melt.updn)
## doesn't look like there is much difference though

require(nlme)
melt.updn$time <- rep(with(melt.updn[1:14,], date - date[1]) + 1, 2)
## include fixed time effect to account for any trend for example?
## use a CAR(1) structure allows for different separations in sampling times 
lmod <- gls(value ~ variable + time, data = melt.updn,
  corr = corCAR1(form=  ~ time | variable))
summary(lmod)
intervals(lmod) ## fitting problems with these dummy data
## test CAR(1) structure - do we need?
lmod2 <- gls(value ~ variable + time, data = melt.updn)
anova(lmod, lmod2) ## no need for the structure here
summary(lmod2) ## looks like no difference in un/restored
anova(lmod2)

Just a few thoughts, without knowing exactly your data and design it is
difficult to say more. With only two groups, it is difficult to more. I
also assume these are dummy data otherwise there really doesn't look
like there is any difference between the two groups of samples.

HTH

G
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

__
R-help@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.

Re: [R] ANOVA

[Tal Galili]

Hi Sarah
I can't see the file.
The general code can be seen in the example on
?aov
(assuming the design is balanced)







On Mon, Jun 29, 2009 at 3:07 PM, Georgina Sarah Humphreys <
g.humphrey...@research.gla.ac.uk> wrote:

> I have the attached data set (csv) and I want to run an analysis of
> variance on the wingsize data (comparing infected vs non-infected) within
> and between experiments.
>
> Can anyone help me with the command I should use?
>
> Many thanks
> Georgina
>
>
>
> PhD Student
> Division of Infection and Immunity
> B5-29, GBRC
> 120 University Place
> Glasgow
> G12 8TA
> Tel: 0141 330 5650
>
> __
> R-help@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.
>
>


-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
http://www.r-statistics.com/
http://www.talgalili.com
http://www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] ANOVA error

[Peter Alspach]

Dani

> I have a dataframe named myMatrix with the structure

If it really is a dataframe, why is it called myMatrix :-)  I guess you 
actually have a factor matrix.  From the message below, Cr would appear to be a 
factor.

There are some fundamental distinctions in R (such as between a matrix and a 
dataframe, or between factor and numeric) which can seem 'subtle' as first.  It 
is worth investing a bit of effort into understanding these. [My apologies if 
this is stating something which is obvious to you.]

HTH 

Peter Alspach

> TreatmentTimeCrmIb ...
> 
> Being the treatment and time the predictors and Cr, mIb and 
> so on the response variables. When I call
> 
> Cr.aov <- aov(Cr~Treatment, data=myMatrix)
> 
> I got this error:
> 
> Error in storage.mode(y) <- "double" :
>   invalid to change the storage mode of a factor In addition: 
> Warning message:
> In model.response(mf, "numeric") :
>   using type="numeric" with a factor response will be ignored
> 
> Can anyone help me? I take the chance to ask another 
> question. I would like to perform an ANOVA on different 
> variables. Is the call aov(Cr+mIb+...~Treatment, 
> data=myMatrix) o correct way of doing it? What I would like 
> to have is a p value for the variance for each response variable.
> Best,
> 
> Dani
> 
> --
> Daniel Valverde Saubí
> 
> Grup de Biologia Molecular de Llevats
> Facultat de Veterinària de la Universitat Autònoma de 
> Barcelona Edifici V, Campus UAB
> 08193 Cerdanyola del Vallès- SPAIN
> 
> Centro de Investigación Biomédica en Red en Bioingeniería, 
> Biomateriales y Nanomedicina (CIBER-BBN)
> 
> Grup d'Aplicacions Biomèdiques de la RMN Facultat de 
> Biociències Universitat Autònoma de Barcelona Edifici Cs, Campus UAB
> 08193 Cerdanyola del Vallès- SPAIN
> +34 93 5814126
> 
> __
> R-help@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.
> 

The contents of this e-mail are privileged and/or confidential to the named
 recipient and are not to be used by any other person and/or organisation.
 If you have received this e-mail in error, please notify the sender and delete
 all material pertaining to this e-mail.

__
R-help@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.

Re: [R] ANOVA help

[Gareth Campbell]

Thanks for that, the error I get when I use manova is this:

Error in summary.manova(fit) : residuals have rank 30 < 31

What does this mean??

When I call the formula I get 31 residuals and 31 residual standard errors.
Not sure why I get this error?

Thanks

Gareth


2008/8/11 S Ellison <[EMAIL PROTECTED]>

> test="Wilks" is available in summary.manova, but not in summary.aov;
> could that be the problem? I find in ?summary.manova the example makes a
> clear distinction between these three summaries...
>
> fit <- manova(Y ~ rate * additive)
> summary.aov(fit)   # univariate ANOVA tables
> summary(fit, test="Wilks") # ANOVA table of Wilks' lambda
> summary(fit)   # same F statistics as single-df terms
>
> Another common variable in anova using R is the way contrasts are set;
> you might check that with your colleague?.
>
> Steve E
>
> >>> "Gareth Campbell" <[EMAIL PROTECTED]> 08/10/08 11:43 PM >>>
> Hi,
>
> I'm doing anova on a matrix of multivariate data where I want to assess
> the
> effect of each column (element).
>
> My matrix is 86 rows x 31 columns.  I've created a grouping factor of
> length
> 86 containing group assignments of 6 types.
>
> Then I run:
>
> x<- aov(matrix~grouping.factor)
> summary(aov.fit.raw, test="Wilks")
>
> This is working fine enough, but I'm getting different results to
> someone
> who I'm comparing with - am I doing what I think I am doing here??
>
> What I think I'm doing is - take the first column (element) and then
> look at
> that element between the 6 groups and report F, Pvalue etc..., then move
> onto the 2nd column and repeat.
>
> Thanks
>
>
>
> --
> Gareth Campbell
> PhD Candidate
> The University of Auckland
>
> P +649 815 3670
> M +6421 256 3511
> E [EMAIL PROTECTED]
> [EMAIL PROTECTED]
>
>[[alternative HTML version deleted]]
>
> __
> R-help@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.
>
>
> ***
> This email and any attachments are confidential. Any u...{{dropped:24}}

__
R-help@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.

Re: [R] ANOVA help

[Gareth Campbell]

I've checked the function code and the error is returned form the following:

if (rss.qr$rank < ncol(resid))
stop(gettextf("residuals have rank %d < %d", rss.qr$rank, ncol(resid)),
domain=NA)

I read this as " If the Residual sums of squares $ rank ... " whatever that
is " is less than the number of columns of residuals, then give me an
error".

Thanks


2008/8/12 Gareth Campbell <[EMAIL PROTECTED]>

> Thanks for that, the error I get when I use manova is this:
>
> Error in summary.manova(fit) : residuals have rank 30 < 31
>
> What does this mean??
>
> When I call the formula I get 31 residuals and 31 residual standard
> errors.  Not sure why I get this error?
>
> Thanks
>
> Gareth
>
>
> 2008/8/11 S Ellison <[EMAIL PROTECTED]>
>
> test="Wilks" is available in summary.manova, but not in summary.aov;
>> could that be the problem? I find in ?summary.manova the example makes a
>> clear distinction between these three summaries...
>>
>> fit <- manova(Y ~ rate * additive)
>> summary.aov(fit)   # univariate ANOVA tables
>> summary(fit, test="Wilks") # ANOVA table of Wilks' lambda
>> summary(fit)   # same F statistics as single-df terms
>>
>> Another common variable in anova using R is the way contrasts are set;
>> you might check that with your colleague?.
>>
>> Steve E
>>
>> >>> "Gareth Campbell" <[EMAIL PROTECTED]> 08/10/08 11:43 PM >>>
>> Hi,
>>
>> I'm doing anova on a matrix of multivariate data where I want to assess
>> the
>> effect of each column (element).
>>
>> My matrix is 86 rows x 31 columns.  I've created a grouping factor of
>> length
>> 86 containing group assignments of 6 types.
>>
>> Then I run:
>>
>> x<- aov(matrix~grouping.factor)
>> summary(aov.fit.raw, test="Wilks")
>>
>> This is working fine enough, but I'm getting different results to
>> someone
>> who I'm comparing with - am I doing what I think I am doing here??
>>
>> What I think I'm doing is - take the first column (element) and then
>> look at
>> that element between the 6 groups and report F, Pvalue etc..., then move
>> onto the 2nd column and repeat.
>>
>> Thanks
>>
>>
>>
>> --
>> Gareth Campbell
>> PhD Candidate
>> The University of Auckland
>>
>> P +649 815 3670
>> M +6421 256 3511
>> E [EMAIL PROTECTED]
>> [EMAIL PROTECTED]
>>
>>[[alternative HTML version deleted]]
>>
>> __
>> R-help@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.
>>
>>
>> ***
>> This email and any attachments are confidential. Any use, copying or
>> disclosure other than by the intended recipient is unauthorised. If
>> you have received this message in error, please notify the sender
>> immediately via +44(0)20 8943 7000 or notify [EMAIL PROTECTED]
>> and delete this message and any copies from your computer and network.
>> LGC Limited. Registered in England 2991879.
>> Registered office: Queens Road, Teddington, Middlesex, TW11 0LY, UK
>>
>
>
>
> --
> Gareth Campbell
> PhD Candidate
> The University of Auckland
>
> P +649 815 3670
> M +6421 256 3511
> E [EMAIL PROTECTED]
> [EMAIL PROTECTED]
>



-- 
Gareth Campbell
PhD Candidate
The University of Auckland

P +649 815 3670
M +6421 256 3511
E [EMAIL PROTECTED]
[EMAIL PROTECTED]

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] ANOVA help

[Gareth Campbell]

Hmm, I have my elements (as in chemical elements) as columns and my samples
as rows.  This is the normal way round for all my other analyses.  I think I
figured it out...  I have 6 groups, all of about 10-12 samples and each with
31 elements.  So for each group there are many more elements than samples
and it makes the analysis impossible due to the problems of the singular
covariance matrix, if I've got that right.
Is that the problem?

2008/8/12 S Ellison <[EMAIL PROTECTED]>

> Rank deficiency is usually an indication of under-determination, if I've
> got it right.
>
> But something is clearly odd, because with 31-column matrix of 86 rows,
> you should have 86 residuals for each of 31 models, not 31 residuals.
> Could your matrix be the wrong way round?
>
>
> >>> "Gareth Campbell" <[EMAIL PROTECTED]> 11/08/2008 23:38 >>>
> Thanks for that, the error I get when I use manova is this:
>
> Error in summary.manova(fit) : residuals have rank 30 < 31
>
> What does this mean??
>
> When I call the formula I get 31 residuals and 31 residual standard
> errors.
> Not sure why I get this error?
>
> Thanks
>
> Gareth
>
>
> 2008/8/11 S Ellison <[EMAIL PROTECTED]>
>
> > test="Wilks" is available in summary.manova, but not in summary.aov;
> > could that be the problem? I find in ?summary.manova the example
> makes a
> > clear distinction between these three summaries...
> >
> > fit <- manova(Y ~ rate * additive)
> > summary.aov(fit)   # univariate ANOVA tables
> > summary(fit, test="Wilks") # ANOVA table of Wilks' lambda
> > summary(fit)   # same F statistics as single-df
> terms
> >
> > Another common variable in anova using R is the way contrasts are
> set;
> > you might check that with your colleague?.
> >
> > Steve E
> >
> > >>> "Gareth Campbell" <[EMAIL PROTECTED]> 08/10/08 11:43 PM >>>
> > Hi,
> >
> > I'm doing anova on a matrix of multivariate data where I want to
> assess
> > the
> > effect of each column (element).
> >
> > My matrix is 86 rows x 31 columns.  I've created a grouping factor
> of
> > length
> > 86 containing group assignments of 6 types.
> >
> > Then I run:
> >
> > x<- aov(matrix~grouping.factor)
> > summary(aov.fit.raw, test="Wilks")
> >
> > This is working fine enough, but I'm getting different results to
> > someone
> > who I'm comparing with - am I doing what I think I am doing here??
> >
> > What I think I'm doing is - take the first column (element) and then
> > look at
> > that element between the 6 groups and report F, Pvalue etc..., then
> move
> > onto the 2nd column and repeat.
> >
> > Thanks
> >
> >
> >
> > --
> > Gareth Campbell
> > PhD Candidate
> > The University of Auckland
> >
> > P +649 815 3670
> > M +6421 256 3511
> > E [EMAIL PROTECTED]
> > [EMAIL PROTECTED]
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-help@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.
> >
> >
> > ***
> > This email and any attachments are confidential. Any use, copying or
> > disclosure other than by the intended recipient is unauthorised. If
> > you have received this message in error, please notify the sender
> > immediately via +44(0)20 8943 7000 or notify [EMAIL PROTECTED]
> > and delete this message and any copies from your computer and
> network.
> > LGC Limited. Registered in England 2991879.
> > Registered office: Queens Road, Teddington, Middlesex, TW11 0LY, UK
> >
>
>
>
> --
> Gareth Campbell
> PhD Candidate
> The University of Auckland
>
> P +649 815 3670
> M +6421 256 3511
> E [EMAIL PROTECTED]
> [EMAIL PROTECTED]
>
> ***
> This email and any attachments are confidential. Any u...{{dropped:24}}

__
R-help@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.

Re: [R] ANOVA help

[Prof Brian Ripley]

On Tue, 12 Aug 2008, Gareth Campbell wrote:


Hmm, I have my elements (as in chemical elements) as columns and my samples
as rows.  This is the normal way round for all my other analyses.  I think I
figured it out...  I have 6 groups, all of about 10-12 samples and each with
31 elements.  So for each group there are many more elements than samples
and it makes the analysis impossible due to the problems of the singular
covariance matrix, if I've got that right.
Is that the problem?


If this is a MANOVA with 31 responses (elements) then the problem is that 
those responses are linearly dependent after subtracting group means (and 
quite possibly before).  With 86 samples this is not normal (MANOVA is 
using a common covariance matrix for each group, not one for each group).


I do wonder what your aim is: this seems to be using MANOVA to do LDA.


2008/8/12 S Ellison <[EMAIL PROTECTED]>


Rank deficiency is usually an indication of under-determination, if I've
got it right.

But something is clearly odd, because with 31-column matrix of 86 rows,
you should have 86 residuals for each of 31 models, not 31 residuals.
Could your matrix be the wrong way round?



"Gareth Campbell" <[EMAIL PROTECTED]> 11/08/2008 23:38 >>>

Thanks for that, the error I get when I use manova is this:

Error in summary.manova(fit) : residuals have rank 30 < 31

What does this mean??

When I call the formula I get 31 residuals and 31 residual standard
errors.
Not sure why I get this error?

Thanks

Gareth


2008/8/11 S Ellison <[EMAIL PROTECTED]>


test="Wilks" is available in summary.manova, but not in summary.aov;
could that be the problem? I find in ?summary.manova the example

makes a

clear distinction between these three summaries...

fit <- manova(Y ~ rate * additive)
summary.aov(fit)   # univariate ANOVA tables
summary(fit, test="Wilks") # ANOVA table of Wilks' lambda
summary(fit)   # same F statistics as single-df

terms


Another common variable in anova using R is the way contrasts are

set;

you might check that with your colleague?.

Steve E


"Gareth Campbell" <[EMAIL PROTECTED]> 08/10/08 11:43 PM >>>

Hi,

I'm doing anova on a matrix of multivariate data where I want to

assess

the
effect of each column (element).

My matrix is 86 rows x 31 columns.  I've created a grouping factor

of

length
86 containing group assignments of 6 types.

Then I run:

x<- aov(matrix~grouping.factor)
summary(aov.fit.raw, test="Wilks")

This is working fine enough, but I'm getting different results to
someone
who I'm comparing with - am I doing what I think I am doing here??

What I think I'm doing is - take the first column (element) and then
look at
that element between the 6 groups and report F, Pvalue etc..., then

move

onto the 2nd column and repeat.

Thanks



--
Gareth Campbell
PhD Candidate
The University of Auckland

P +649 815 3670
M +6421 256 3511
E [EMAIL PROTECTED]
[EMAIL PROTECTED]

   [[alternative HTML version deleted]]

__
R-help@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.


***
This email and any attachments are confidential. Any use, copying or
disclosure other than by the intended recipient is unauthorised. If
you have received this message in error, please notify the sender
immediately via +44(0)20 8943 7000 or notify [EMAIL PROTECTED]
and delete this message and any copies from your computer and

network.

LGC Limited. Registered in England 2991879.
Registered office: Queens Road, Teddington, Middlesex, TW11 0LY, UK





--
Gareth Campbell
PhD Candidate
The University of Auckland

P +649 815 3670
M +6421 256 3511
E [EMAIL PROTECTED]
[EMAIL PROTECTED]

***
This email and any attachments are confidential. Any u...{{dropped:24}}


__
R-help@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.



--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@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.

[R] anova() or aov()?

[jass]

Dear all,

I have a simple (I think) question that is troubling me lately:

Is there any main difference between anova() command and aov() command when 
performing an ANOVA in Experimental design 
analyses?

Thank you for your time,

Ismini

__
R-help@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.

[R] ANOVA in R

[Amit Patel]

Hi

I Have a very large dataset that I would like to conduct ANOVA tests on. Im not 
a very strong programmer so any help would be appreciated.
the format is 

Identifier             A1       A2        B1      
B2       C1   C2      Norm1         Norm2
1234                  1        1            
NA     NA      4       3        
NA               NA
4567                  2        
2              4      4         8       
8       9                    9


and so on
I have 10 runs for 3 different doses plus the normal state. Any help greatly 
appreciated

 
 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 





  
[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA in R

[Samor Gandhi]

Hi, I'm using a repeated measures ANOVA in R using lme(). The SAS code would 
be: PROC MIXED DATA=[data set below]; CLASS pid treat period time seq; 
MODEL Y = seq period treat time treat*time; REPEATED time / SUBJECT=pid 
TYPE=cs;RUN, I donot have SAS, instead I have R and I would like to try the 
following:anova(lme(response ~ seq period treat time treat*time,random= ~1|SUB, 
   correlation=corCompSymm())) Is this correct? Can I also write the model 
as Y_ijklt = m + a_l + b_k + c_j + d_t + (cd)_jt + u_ijklt Y_ijklt is the 
response variable due to pid i, treat j, period k, seq l, and time t.  Thank 
you very much in advance for your help :)



Samor




  
[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA in R

[Samor Gandhi]

Hi,I'm using a repeated measures ANOVA in R using lme(). The SAS code would be: 
 
PROC MIXED DATA=[data set below];
 CLASS pid treat period time seq;
 MODEL Y = seq period treat time treat*time;
 REPEATED time / SUBJECT=pid TYPE=cs;
RUN,  I donot have SAS, instead I have R and I would like to try the following:
anova(lme(response ~ seq period treat time treat*time,random= ~1|SUB,    
correlation=corCompSymm()))

Is this correct? Can I also write the model as

Y_ijklt = m + a_l + b_k + c_j + d_t + (cd)_jt + u_ijkltY_ijklt is the response 
variable due to pid i, treat j, period k, seq l, and time t. Thank you very 
much in advance for your help :)
Samor 






  
[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA and lmer

[Eric Imbert]

I am analyzing from a very simple experiment.
I have measured plants of two different colours (yellow and purple) in 9
different populations.

So, I have two different factors : a fixed effect (Colour with two
levels) and a random one (Population with 9 levels).

I first analyzed the data with the aov function
LargS is the variable 

aov(formula = LargS ~ Col + Error(Col/Pop))

Terms:
Col
Sum of Squares  3.440351
Deg. of Freedom1

Estimated effects are balanced

Stratum 2: Col:Pop

Terms:
Residuals
Sum  of Squares   3017.112
Deg. of Freedom16

Residual standard error: 13.73206 

Stratum 3: Within

Terms:
Residuals
Sum of Squares   3347.385
Deg. of Freedom   302

To test for the interaction Col*Pop, I used the following F-ratio =
(3017/16)/(3347/302) = 188. Highly significant !


Now, let's go to the analysis performed by lmer - First I do the linear
model without the Col*Pop interaction :
m3=lmer(LargS ~ Col + (1 | Pop)

And next with the interaction : m2=lmer(LargS ~ Col + (Col | Pop))

Comparing both models : anova(m2,m3) :

   Df AIC BIC  logLik  Chisq Chi Df Pr(>Chisq)
m3  3 1710.67 1721.97 -852.33 
m2  5 1714.59 1733.43 -852.30 0.0746  2 0.9634

=> Conclusion : the interaction Col*Pop is not significant !

I guess I am missing something.
Who can help ?

Eric


[[alternative HTML version deleted]]

__
R-help@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.

[R] anova power calculations

[Will Holcomb]

I sent a message a couple days ago about doing calculations for power of the
ANOVA. Several people got back to me very quickly which I really
appreciated.

I'm working now on a similar problem, but instead of a balanced ANOVA, I
have an unbalanced one. The first part of the question was:

You assume that the within-population standard deviations all equal 9. You
set the Type 1 error rate at á = .05. You presume that the population means
will have the following values: uA = 17.5, uB = 19, uC = 25, and uD = 20.5.
You intend to run 80 subjects in all, with equal n's across all 4 groups.
You plan on conducting a one-way ANOVA. Compute your power to reject the
null hypothesis under these conditions.

I did:

within.var = 9 ^ 2
means = c(17.5, 19, 25, 20.5)
N = 80
J = length(means)
power.anova.test(groups = J, n = N / J,
 between.var = var(means),
 within.var = within.var,
 sig.level = 0.05)

This gives me 0.6155 which agrees with SAS. The next problem though is:

You have the same Type 1 error rate and make the same assumptions about the
population standard deviation and the population means as in part a. You
still have 80 subjects in all but now you want to know how power might
change by running 10 subjects in groups A, B, and D and 50 subjects in group
C. Determine the power under this subject allocation scheme.

For this one I am doing:

# Quantile of the cutoff point in the central F
central.quant = qf(.05, J - 1, N - J, lower.tail = FALSE)
weighted.means = data.frame(Mean = means, Weight = c(10, 10, 50, 10))
# Noncentrality parameter for unbalanced ANOVA
noncentral.param = 0
for(i in 1:length(weighted.means$Mean)) {
  noncentral.param = (noncentral.param + weighted.means$Weight[i] *
  (weighted.means$Mean[i] - mean(weighted.means$Mean)) ^
2)
}
noncentral.param = noncentral.param / within.var
# Probability of central quantile in noncentral distribution
noncentral.p = pf(central.quant, J - 1, N - J, noncentral.param,
lower.tail= FALSE)
noncentral.p

The logic behind this is in my assignment at:

http://odin.himinbi.org/classes/psy304b/homework_2.xhtml#p2b

This works for a balanced ANOVA and gives the same result as
power.anova.test (and SAS). For the unbalanced ANOVA though it is giving me
a different result though than SAS, 0.8759455 versus 0.680.

So is there a straightforward way to compute the power of an unbalanced
ANOVA? If there isn't, does anyone have any idea what is wrong with my code?
The SAS I am comparing it to is:

Data Dep;
Input cue $ mean uneven_weight;
datalines;
A 17.5  1
B 191
C 255
D 20.5  1
;

proc glmpower;
class cue;
model mean = cue;
weight uneven_weight;
power
stddev = 9
alpha = 0.05
ntotal= 80
power = .;
run;

Any help would be much appreciated.

Will

[[alternative HTML version deleted]]

__
R-help@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.

[R] anova (repeated measure)

[guillaume chaumet]

Hi R people,
I'm trying to do a Anova with 3 different tests (intra subject), 6 times of
measures (intra subject) and 3 groups (inter subject). I want to obtain an
Huynh-Feldt epsilon on accuracy with >anova(lm(acc.t1~gr.dem),test="Spherical")
for test 1 # gr.dem is my group variable
>anova(lm(acc.t2~gr.dem),test="Spherical") for test 2
>anova(lm(acc.t3~gr.dem),test="Spherical") for test 3

acc.t1, acc.t2 and acc.t3 was data frames with 6 columns each

the results was something like this:

Analysis of Variance Table

Greenhouse-Geisser epsilon: 0.2746
Huynh-Feldt epsilon:0.2915

 Df F num Df den Df Pr(>F) G-G Pr H-F Pr
(Intercept)   1 7352.7017  6150 1.640e-182  6.894e-52  6.269e-55
gr.dem10.6222  61500.712300.511700.52064
Residuals25

And so what could I do with my epsilon?
Previous aov() revealed me that there no effects of my groups but is my 6
times measures factor and my test factor were always significant?

Guillaume

[[alternative HTML version deleted]]

__
R-help@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.

[R] anova repeated measure

[guillaume chaumet]

Dear R people,
I'm new user of R and I'm trying to do a Anova with 3 different tests (intra
subject), 6 times of measures (intra subject) and 3 groups (inter subject).
I want to obtain an Huynh-Feldt epsilon on accuracy with:
>anova(lm(acc.t1~gr.dem),test="Spherical") for test 1 # gr.dem is my group
variable
>anova(lm(acc.t2~gr.dem),test="Spherical") for test 2
>anova(lm(acc.t3~gr.dem),test="Spherical") for test 3

acc.t1, acc.t2 and acc.t3 was data frames with 6 columns each

the results was something like this:

Analysis of Variance Table

Greenhouse-Geisser epsilon: 0.2746
Huynh-Feldt epsilon:0.2915

 Df F num Df den Df Pr(>F) G-G Pr H-F Pr
(Intercept)   1 7352.7017  6150 1.640e-182  6.894e-52  6.269e-55
gr.dem10.6222  61500.712300.511700.52064
Residuals25

And so what could I do with my epsilon?
Previous aov() revealed me that there no effects of my groups but is my 6
times measures factor and my test factor were always significant?


Guillaume

[[alternative HTML version deleted]]

__
R-help@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.

[R] anova and contratst

[Colm G.Connolly]

Hi everybody,

I'm using the gmodels package to convert human readable contrasts  
into the format required by R and would be grateful if someone could  
confirm for me whether I've got the contrasts right in the sample  
code below.

I'm working on the assumption that the contrasts are index according  
to the way that levels reports them for a factor.
In my case levels(roi.errs$Group) reports ctrl, long, short in that  
order so I'm assuming that to compare ctrl to short the correct  
contrast is c(1, 0, -1). Am I correct?

I tried to perform 3 contrasts with my data and make.contrasts  
complained about there being too many contrasts specified. Is this  
because the number of contrasts can only be less than or equal to the  
number of degrees of freedom that go into calculation of the mean  
square for the Group term?

Is the only way to get the short vs long contrast to sacrifice one of  
the other contrasts?

rm(list=ls());
library(gmodels);

Group = c("ctrl", "ctrl", "ctrl", "short", "short", "short", "long",  
"long", "long")
Subject = c("66101_2", "66112_2", "66118_2", "66119", "66121",  
"66123", "66026_3", "66030_2", "66044_3")
Mean_1 = c(-0.399281, 6.951385, 5.986774, 3.485916, 7.081511,  
1.705302, 7.859186, 4.573201, 3.931118)
Mean_2 = c(0.350127, 5.863118, 3.826101, 0.999821, 9.099477,  
2.376836, 5.384967, 5.950012, 7.688557)
roi.errs = data.frame(Group, Subject, Mean_1, Mean_2)

# planned contrasts
cmat = rbind(
   "ctrl vs. short" = c(1,  0, -1),
   "ctrl vs. long"  = c(1, -1,  0));
groupContrasts=make.contrasts(cmat);
summarySplit=list(Group=list("ctrl vs. short" = 1, "ctrl vs. long"=2));

roi.err.aov<-aov(Mean_1 ~ Group, contrasts=list 
("Group"=groupContrasts), data=roi.errs);
print(summary(roi.err.aov, split=summarySplit))
with(roi.errs, print(pairwise.t.test(Mean_1, Group)))

Thanks in advance,

--
Dr Colm G. Connolly
School of Psychology and Institute of Neuroscience
The Lloyd Building
University of Dublin
Trinity College, Dublin 2, Éire
Tel: +353-1-896-8475
Fax: +353-1-671-3183

__
R-help@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.

Re: [R] anova help

[Mike Lawrence]

Provide some example data & code and someone may be able to help.

On Sat, Feb 14, 2009 at 2:46 PM, Joe King  wrote:
> Hi all, I am trying to run a two factor anova, but one of the factors is a
> random factor, now I am also running in SPSS and it seems its dividing by
> the wrong term to get the appropriate F term. here is my data. In SPSS the F
> scores about double the ones in R, how can I specify one of my factors as a
> random factor or change it to where it does the right model fitting? I am
> using the lm command instead of glm. I am new to R so this might seem basic.
>
>
>
> Joe King, M.A.
>   j...@joepking.com
>
> "Never give in, never give in, never; never; never; never - in nothing,
> great or small, large or petty - never give in except to convictions of
> honor and good sense" - Winston Churchill
>
>
>
> "You have enemies? Good. That means you've stood up for something, sometime
> in your life." - Winston Churchill
>
>
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@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.
>



-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University
www.thatmike.com

Looking to arrange a meeting? Check my public calendar:
http://www.thatmike.com/mikes-public-calendar

~ Certainty is folly... I think. ~

__
R-help@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.

Re: [R] anova help

[Tal Galili]

Hi Joe,
you might want to have a look at the nlme package with the lme command.
Another option is the more advanced lmer package.

Lastly, you could have a look at the ?aov command, and notice the option of
using the +Error() term (but that would only work for balanced design cases,
so I've heard, so be aware)

Tal



On Sat, Feb 14, 2009 at 8:46 PM, Joe King  wrote:

> Hi all, I am trying to run a two factor anova, but one of the factors is a
> random factor, now I am also running in SPSS and it seems its dividing by
> the wrong term to get the appropriate F term. here is my data. In SPSS the
> F
> scores about double the ones in R, how can I specify one of my factors as a
> random factor or change it to where it does the right model fitting? I am
> using the lm command instead of glm. I am new to R so this might seem
> basic.
>
>
>
> Joe King, M.A.
>   j...@joepking.com
>
> "Never give in, never give in, never; never; never; never - in nothing,
> great or small, large or petty - never give in except to convictions of
> honor and good sense" - Winston Churchill
>
>
>
> "You have enemies? Good. That means you've stood up for something, sometime
> in your life." - Winston Churchill
>
>
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@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.
>



-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] anova help

[Joe King]

Ok well heres my data, with the second variable (B variable) being the
random variable and the other fixed. Although the F values are about half of
what SPSS puts out.

 

My code is anova(lm(dependentrandom~typemusic+typemusic*musicselection))

 

This is just dummy data for a class but I am trying to use the data I am
running in SPSS to learn R. I am also in an R class but we are not going to
learn ANOVA.

 

Joe King, M.A.
 <mailto:j...@joepking.com> j...@joepking.com

"Never give in, never give in, never; never; never; never - in nothing,
great or small, large or petty - never give in except to convictions of
honor and good sense" - Winston Churchill

 

"You have enemies? Good. That means you've stood up for something, sometime
in your life." - Winston Churchill

 

From: Tal Galili [mailto:tal.gal...@gmail.com] 
Sent: Sunday, February 15, 2009 10:25 AM
To: Joe King
Cc: r-help@r-project.org
Subject: Re: [R] anova help

 

Hi Joe,

you might want to have a look at the nlme package with the lme command.

Another option is the more advanced lmer package.

 

Lastly, you could have a look at the ?aov command, and notice the option of
using the +Error() term (but that would only work for balanced design cases,
so I've heard, so be aware)

 

Tal

 

 

On Sat, Feb 14, 2009 at 8:46 PM, Joe King  wrote:

Hi all, I am trying to run a two factor anova, but one of the factors is a
random factor, now I am also running in SPSS and it seems its dividing by
the wrong term to get the appropriate F term. here is my data. In SPSS the F
scores about double the ones in R, how can I specify one of my factors as a
random factor or change it to where it does the right model fitting? I am
using the lm command instead of glm. I am new to R so this might seem basic.



Joe King, M.A.
 <mailto:j...@joepking.com> j...@joepking.com

"Never give in, never give in, never; never; never; never - in nothing,
great or small, large or petty - never give in except to convictions of
honor and good sense" - Winston Churchill



"You have enemies? Good. That means you've stood up for something, sometime
in your life." - Winston Churchill




   [[alternative HTML version deleted]]

__
R-help@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.




-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il



__
R-help@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.

Re: [R] anova help

[David Winsemius]

You probably did not forget to attach the data the first time, either.  
You, instead, probably did not remember to read the posting guide  
thoroughly regarding what constitutes an acceptable attachment type.  
Suggest that you change your mail client settings to transmit plain  
text and "follow the directions":


http://www.R-project.org/posting-guide.html

--
David Winsemius


On Feb 15, 2009, at 2:21 PM, Joe King wrote:


Ok well heres my data, with the second variable (B variable) being the
random variable and the other fixed. Although the F values are about  
half of

what SPSS puts out.



My code is anova(lm(dependentrandom~typemusic 
+typemusic*musicselection))




This is just dummy data for a class but I am trying to use the data  
I am
running in SPSS to learn R. I am also in an R class but we are not  
going to

learn ANOVA.



Joe King, M.A.
<mailto:j...@joepking.com> j...@joepking.com

"Never give in, never give in, never; never; never; never - in  
nothing,
great or small, large or petty - never give in except to convictions  
of

honor and good sense" - Winston Churchill



"You have enemies? Good. That means you've stood up for something,  
sometime

in your life." - Winston Churchill



From: Tal Galili [mailto:tal.gal...@gmail.com]
Sent: Sunday, February 15, 2009 10:25 AM
To: Joe King
Cc: r-help@r-project.org
Subject: Re: [R] anova help



Hi Joe,

you might want to have a look at the nlme package with the lme  
command.


Another option is the more advanced lmer package.



Lastly, you could have a look at the ?aov command, and notice the  
option of
using the +Error() term (but that would only work for balanced  
design cases,

so I've heard, so be aware)



Tal





On Sat, Feb 14, 2009 at 8:46 PM, Joe King  wrote:

Hi all, I am trying to run a two factor anova, but one of the  
factors is a
random factor, now I am also running in SPSS and it seems its  
dividing by
the wrong term to get the appropriate F term. here is my data. In  
SPSS the F
scores about double the ones in R, how can I specify one of my  
factors as a
random factor or change it to where it does the right model fitting?  
I am
using the lm command instead of glm. I am new to R so this might  
seem basic.




Joe King, M.A.
<mailto:j...@joepking.com> j...@joepking.com

"Never give in, never give in, never; never; never; never - in  
nothing,
great or small, large or petty - never give in except to convictions  
of

honor and good sense" - Winston Churchill



"You have enemies? Good. That means you've stood up for something,  
sometime

in your life." - Winston Churchill




  [[alternative HTML version deleted]]

__
R-help@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.




--
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il



__
R-help@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.


__
R-help@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.

Re: [R] Anova Help?

[Tal Galili]

Hi Dar.

I am not sure I got your question -
Are you asking what analysis to perform ?
Or how to perform it ?

Could you please give more details ?




On Wed, Feb 25, 2009 at 5:52 PM, Dar  wrote:

> I am conducting an experiment where students are put into 5 total
> groups (one is the control group).  They are given a task and then I'm
> measuring if there are differences (A 5 X 1 between-subject design
> will be used for the experiment).  I'm a little confused on the data
> explanation (or should I say how do I explain what is being analyzed
> versus just comparing values)  Any help?
>
> __
> R-help@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.
>



-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova Help?

[Dar]

I’m just setting up the experiment and need help explaining what the
data analysis would be.  Let me know of any questions….. what would be
compared and how it would be measured.  Is it a multi-way or 1-way
anova?
Thanks!

A 5 X 1 between-subject design will be used for the experiment. Four
tool groups (Group1, Group2, Group3, Group4) will be able to view one
category of information that their groups are allowed to in their
tools. The control group, however, will have no access to the
experimental tool.

Pre-survey has demographic and previous experience questionnaire.

The independent variable is the exposure to the information. The
independent variable has five levels: no information, appearance
information, educational information, contact information, and
personal information.

The dependent variable will be (2 dependent variables) gathered
through a questionnaire at the end of the study.
Connectedness. This variable will be operationalized as the
Connectedness score
Learning. This variable will be operationalized as the Learning score

Covariates:
Frequency of using tool
Duration of tool use
Previous experience (three types of previous experience).


On Feb 25, 11:35 am, Tal Galili  wrote:
> Hi Dar.
>
> I am not sure I got your question -
> Are you asking what analysis to perform ?
> Or how to perform it ?
>
> Could you please give more details ?
>
> On Wed, Feb 25, 2009 at 5:52 PM, Dar  wrote:
> > I am conducting an experiment where students are put into 5 total
> > groups (one is the control group).  They are given a task and then I'm
> > measuring if there are differences (A 5 X 1 between-subject design
> > will be used for the experiment).  I'm a little confused on the data
> > explanation (or should I say how do I explain what is being analyzed
> > versus just comparing values)  Any help?
>
> > __
> > r-h...@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.
>
> --
> --
>
> My contact information:
> Tal Galili
> Phone number: 972-50-3373767
> FaceBook: Tal Galili
> My Blogs:www.talgalili.comwww.biostatistics.co.il
>
>         [[alternative HTML version deleted]]
>
> __
> r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@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.

Re: [R] Anova Help?

[JLucke]

Based on your description so far and on my making a reasonable set of 
assumptions that would still need to be checked,  your analysis would 
appear to be a multivariate (connectedness, learning) analysis of 
covariance of a  5 (information exposure) by 3 (previous experience) 
design with 2 covariates (frequency, duration).  You really should consult 
one of the myriad of analysis of variance texts.

Joseph F. Lucke
Senior Statistician
Research Institute on Addictions
University at Buffalo
SUNY




Dar  
Sent by: r-help-boun...@r-project.org
02/25/2009 01:22 PM

To
r-help@r-project.org
cc

Subject
Re: [R] Anova Help?






I’m just setting up the experiment and need help explaining what the
data analysis would be.  Let me know of any questions….. what would be
compared and how it would be measured.  Is it a multi-way or 1-way
anova?
Thanks!

A 5 X 1 between-subject design will be used for the experiment. Four
tool groups (Group1, Group2, Group3, Group4) will be able to view one
category of information that their groups are allowed to in their
tools. The control group, however, will have no access to the
experimental tool.

Pre-survey has demographic and previous experience questionnaire.

The independent variable is the exposure to the information. The
independent variable has five levels: no information, appearance
information, educational information, contact information, and
personal information.

The dependent variable will be (2 dependent variables) gathered
through a questionnaire at the end of the study.
Connectedness. This variable will be operationalized as the
Connectedness score
Learning. This variable will be operationalized as the Learning score

Covariates:
Frequency of using tool
Duration of tool use
Previous experience (three types of previous experience).


On Feb 25, 11:35 am, Tal Galili  wrote:
> Hi Dar.
>
> I am not sure I got your question -
> Are you asking what analysis to perform ?
> Or how to perform it ?
>
> Could you please give more details ?
>
> On Wed, Feb 25, 2009 at 5:52 PM, Dar  wrote:
> > I am conducting an experiment where students are put into 5 total
> > groups (one is the control group).  They are given a task and then I'm
> > measuring if there are differences (A 5 X 1 between-subject design
> > will be used for the experiment).  I'm a little confused on the data
> > explanation (or should I say how do I explain what is being analyzed
> > versus just comparing values)  Any help?
>
> > __
> > r-h...@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.
>
> --
> --
>
> My contact information:
> Tal Galili
> Phone number: 972-50-3373767
> FaceBook: Tal Galili
> My Blogs:www.talgalili.comwww.biostatistics.co.il
>
> [[alternative HTML version deleted]]
>
> __
> r-h...@r-project.org mailing 
listhttps://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting 
guidehttp://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@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.



[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova Help?

[Tal Galili]

Hi Dar.

could you please specify what type of variables the dependent and
independent variables are ? (numeric/ordered/categorical , what is their
range)
could you please specify how many dependent and independent variables you
have ? (numeric/ordered/categorical , what is their range)

In either way - reading more on anova is a good way to go.


Tal






On Wed, Feb 25, 2009 at 8:22 PM, Dar  wrote:

> I’m just setting up the experiment and need help explaining what the
> data analysis would be.  Let me know of any questions….. what would be
> compared and how it would be measured.  Is it a multi-way or 1-way
> anova?
> Thanks!
>
> A 5 X 1 between-subject design will be used for the experiment. Four
> tool groups (Group1, Group2, Group3, Group4) will be able to view one
> category of information that their groups are allowed to in their
> tools. The control group, however, will have no access to the
> experimental tool.
>
> Pre-survey has demographic and previous experience questionnaire.
>
> The independent variable is the exposure to the information. The
> independent variable has five levels: no information, appearance
> information, educational information, contact information, and
> personal information.
>
> The dependent variable will be (2 dependent variables) gathered
> through a questionnaire at the end of the study.
> Connectedness. This variable will be operationalized as the
> Connectedness score
> Learning. This variable will be operationalized as the Learning score
>
> Covariates:
> Frequency of using tool
> Duration of tool use
> Previous experience (three types of previous experience).
>
>
> On Feb 25, 11:35 am, Tal Galili  wrote:
> > Hi Dar.
> >
> > I am not sure I got your question -
> > Are you asking what analysis to perform ?
> > Or how to perform it ?
> >
> > Could you please give more details ?
> >
> > On Wed, Feb 25, 2009 at 5:52 PM, Dar  wrote:
> > > I am conducting an experiment where students are put into 5 total
> > > groups (one is the control group).  They are given a task and then I'm
> > > measuring if there are differences (A 5 X 1 between-subject design
> > > will be used for the experiment).  I'm a little confused on the data
> > > explanation (or should I say how do I explain what is being analyzed
> > > versus just comparing values)  Any help?
> >
> > > __
> > > r-h...@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.
> >
> > --
> > --
> >
> > My contact information:
> > Tal Galili
> > Phone number: 972-50-3373767
> > FaceBook: Tal Galili
> > My Blogs:www.talgalili.comwww.biostatistics.co.il
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > r-h...@r-project.org mailing listhttps://
> stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guidehttp://
> www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@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.
>



-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova Help?

[Dar]

Independent Variables: 5 different levels of exposure to information:
no information, appearance
information, educational information, contact information, and
personal information.  The types Unranked Categorical Variable

Dependent Variables: 2 and are numerical

The dependent variable will be (2 dependent variables) gathered
through a questionnaire at the end of the study.
Connectedness. This variable will be operationalized as the
Connectedness score
Learning. This variable will be operationalized as the Learning score

On Feb 25, 3:24 pm, Tal Galili  wrote:
> Hi Dar.
>
> could you please specify what type of variables the dependent and
> independent variables are ? (numeric/ordered/categorical , what is their
> range)
> could you please specify how many dependent and independent variables you
> have ? (numeric/ordered/categorical , what is their range)
>
> In either way - reading more on anova is a good way to go.
>
> Tal
>
>
>
> On Wed, Feb 25, 2009 at 8:22 PM, Dar  wrote:
> > I’m just setting up the experiment and need help explaining what the
> > data analysis would be.  Let me know of any questions….. what would be
> > compared and how it would be measured.  Is it a multi-way or 1-way
> > anova?
> > Thanks!
>
> > A 5 X 1 between-subject design will be used for the experiment. Four
> > tool groups (Group1, Group2, Group3, Group4) will be able to view one
> > category of information that their groups are allowed to in their
> > tools. The control group, however, will have no access to the
> > experimental tool.
>
> > Pre-survey has demographic and previous experience questionnaire.
>
> > The independent variable is the exposure to the information. The
> > independent variable has five levels: no information, appearance
> > information, educational information, contact information, and
> > personal information.
>
> > The dependent variable will be (2 dependent variables) gathered
> > through a questionnaire at the end of the study.
> > Connectedness. This variable will be operationalized as the
> > Connectedness score
> > Learning. This variable will be operationalized as the Learning score
>
> > Covariates:
> > Frequency of using tool
> > Duration of tool use
> > Previous experience (three types of previous experience).
>
> > On Feb 25, 11:35 am, Tal Galili  wrote:
> > > Hi Dar.
>
> > > I am not sure I got your question -
> > > Are you asking what analysis to perform ?
> > > Or how to perform it ?
>
> > > Could you please give more details ?
>
> > > On Wed, Feb 25, 2009 at 5:52 PM, Dar  wrote:
> > > > I am conducting an experiment where students are put into 5 total
> > > > groups (one is the control group).  They are given a task and then I'm
> > > > measuring if there are differences (A 5 X 1 between-subject design
> > > > will be used for the experiment).  I'm a little confused on the data
> > > > explanation (or should I say how do I explain what is being analyzed
> > > > versus just comparing values)  Any help?
>
> > > > __
> > > > r-h...@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.
>
> > > --
> > > --
>
> > > My contact information:
> > > Tal Galili
> > > Phone number: 972-50-3373767
> > > FaceBook: Tal Galili
> > > My Blogs:www.talgalili.comwww.biostatistics.co.il
>
> > >         [[alternative HTML version deleted]]
>
> > > __
> > > r-h...@r-project.org mailing listhttps://
> > stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guidehttp://
> >www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
>
> > __
> > r-h...@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.
>
> --
> --
>
> My contact information:
> Tal Galili
> Phone number: 972-50-3373767
> FaceBook: Tal Galili
> My Blogs:www.talgalili.comwww.biostatistics.co.il
>
>         [[alternative HTML version deleted]]
>
> __
> r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@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.

Re: [R] Anova Help?

[Tal Galili]

Hi Dar.
It sounds like anova is a good way to go.
The command to look into in R is: ?aov
you might also go into manova (for anova for two dependent variables
together)
But make sure you understand the background of it when you work the
commands.


Cheers,
Tal






On Wed, Feb 25, 2009 at 11:06 PM, Dar  wrote:

> Independent Variables: 5 different levels of exposure to information:
> no information, appearance
> information, educational information, contact information, and
> personal information.  The types Unranked Categorical Variable
>
> Dependent Variables: 2 and are numerical
>
> The dependent variable will be (2 dependent variables) gathered
> through a questionnaire at the end of the study.
> Connectedness. This variable will be operationalized as the
> Connectedness score
> Learning. This variable will be operationalized as the Learning score
>
> On Feb 25, 3:24 pm, Tal Galili  wrote:
> > Hi Dar.
> >
> > could you please specify what type of variables the dependent and
> > independent variables are ? (numeric/ordered/categorical , what is their
> > range)
> > could you please specify how many dependent and independent variables you
> > have ? (numeric/ordered/categorical , what is their range)
> >
> > In either way - reading more on anova is a good way to go.
> >
> > Tal
> >
> >
> >
> > On Wed, Feb 25, 2009 at 8:22 PM, Dar  wrote:
> > > I’m just setting up the experiment and need help explaining what the
> > > data analysis would be.  Let me know of any questions….. what would be
> > > compared and how it would be measured.  Is it a multi-way or 1-way
> > > anova?
> > > Thanks!
> >
> > > A 5 X 1 between-subject design will be used for the experiment. Four
> > > tool groups (Group1, Group2, Group3, Group4) will be able to view one
> > > category of information that their groups are allowed to in their
> > > tools. The control group, however, will have no access to the
> > > experimental tool.
> >
> > > Pre-survey has demographic and previous experience questionnaire.
> >
> > > The independent variable is the exposure to the information. The
> > > independent variable has five levels: no information, appearance
> > > information, educational information, contact information, and
> > > personal information.
> >
> > > The dependent variable will be (2 dependent variables) gathered
> > > through a questionnaire at the end of the study.
> > > Connectedness. This variable will be operationalized as the
> > > Connectedness score
> > > Learning. This variable will be operationalized as the Learning score
> >
> > > Covariates:
> > > Frequency of using tool
> > > Duration of tool use
> > > Previous experience (three types of previous experience).
> >
> > > On Feb 25, 11:35 am, Tal Galili  wrote:
> > > > Hi Dar.
> >
> > > > I am not sure I got your question -
> > > > Are you asking what analysis to perform ?
> > > > Or how to perform it ?
> >
> > > > Could you please give more details ?
> >
> > > > On Wed, Feb 25, 2009 at 5:52 PM, Dar  wrote:
> > > > > I am conducting an experiment where students are put into 5 total
> > > > > groups (one is the control group).  They are given a task and then
> I'm
> > > > > measuring if there are differences (A 5 X 1 between-subject design
> > > > > will be used for the experiment).  I'm a little confused on the
> data
> > > > > explanation (or should I say how do I explain what is being
> analyzed
> > > > > versus just comparing values)  Any help?
> >
> > > > > __
> > > > > r-h...@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.
> >
> > > > --
> > > > --
> >
> > > > My contact information:
> > > > Tal Galili
> > > > Phone number: 972-50-3373767
> > > > FaceBook: Tal Galili
> > > > My Blogs:www.talgalili.comwww.biostatistics.co.il
> >
> > > > [[alternative HTML version deleted]]
> >
> > > > __
> > > > r-h...@r-project.org mailing listhttps://
> > > stat.ethz.ch/mailman/listinfo/r-help
> > > > PLEASE do read the posting guidehttp://
> > >www.R-project.org/posting-guide.html
> > > > and provide commented, minimal, self-contained, reproducible code.
> >
> > > __
> > > r-h...@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.
> >
> > --
> > --
> >
> > My contact information:
> > Tal Galili
> > Phone number: 972-50-3373767
> > FaceBook: Tal Galili
> > My Blogs:www.talgalili.comwww.biostatistics.co.il
> >
> > [[alternative HTML version deleted]]
> >
> > ___

[R] ANOVA/statistics question

[drmh]

(Have searched for this already)

Hi,

How do you find the strength of correlation between two variables using an
ANOVA table?  "Pr(>F)" gives the statistical significance of the
association, but not the strength of the correlation.

See data (from R) below

Readable:
  "Df" "Sum Sq""Mean Sq"   "F value"

"Pr(>F)"
"pre"  1   0.00593  0.00593936 
0.7450563   0.401636958677004
"coh" 1   0.04311  0.04311302  5.4082639
  
0.0344751749542619
"Residuals"15 0.11957  0.00797169  NA   
 
NA

Original:
"Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
"pre" 1 0.0059393604629317 0.0059393604629317 0.745056336657567
0.401636958677004
"coh" 1 0.0431130207164516 0.0431130207164516 5.40826398359156
0.0344751749542619
"Residuals" 15 0.119575396598395 0.00797169310655964 NA NA

Any help would be greatly appreciated,
Douglas Holmes
-- 
View this message in context: 
http://www.nabble.com/ANOVA-statistics-question-tp23231563p23231563.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.

[R] anova(cph(..) output

[pompon]

Hello,

I am a beginner in R and statistics, so my question may be trivial. Sorry in
advance.
I performed a Cox proportion hazard regression with 2 categorical variables
with cph{design}. Then an anova on the results.
the output is 

> anova(cph(surv(survival, censor) ~ plant + leaf.age + plant*leaf.age,
> Mpnymph)

Wald Statistics  Response: Surv(survival, censored) 

 FactorChi-Square d.f. P
 
 plant  (Factor+Higher Order Factors) 96.96 12   <.0001
  All Interactions   10.58 
6   0.1022
 leaf.age  (Factor+Higher Order Factors)  29.11  7   0.0001
  All Interactions 10.58 
6   0.1022
 plant * leaf.age  (Factor+Higher Order Factors)  10.58  6   0.1022
 TOTAL   106.63 13   <.0001

What do "All interaction" stand for?
The real df of for plant is 6 and 1 for leaf.age. Then, which chi square is
one for my main factors anf their interaction.

thank you,
Julien.
-- 
View this message in context: 
http://www.nabble.com/anova%28cph%28..%29-output-tp23563818p23563818.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.

[R] ANOVA between & within variance

[Gregor Rolshausen]

hi,
is there an option to calculate the 'within' & 'between' group variances 
for a simple ANOVA (aov) model (2 groups, 1 trait, normally distr.) ?

or do I have to calculate them from the Sum Sq ?

thanks for your time and greetings,

gregor


--
Gregor Rolshausen

PhD Student; University of Freiburg, Germany

e-mail: [EMAIL PROTECTED]   
tel.  : ++49 761 2032559

__
R-help@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.

[R] anova planned comparisons/contrasts

[Tyler Smith]

Hi,

I'm trying to figure out how anova works in R by translating the
examples in Sokal And Rohlf's (1995 3rd edition) Biometry. I've hit a
snag with planned comparisons, their box 9.4 and section 9.6. It's a
basic anova design:

treatment <- factor(rep(c("control", "glucose", "fructose", 
  "gluc+fruct", "sucrose"), each = 10))

length <- c(75, 67, 70, 75, 65, 71, 67, 67, 76, 68,
57, 58, 60, 59, 62, 60, 60, 57, 59, 61,
58, 61, 56, 58, 57, 56, 61, 60, 57, 58,
58, 59, 58, 61, 57, 56, 58, 57, 57, 59,
62, 66, 65, 63, 64, 62, 65, 65, 62, 67)

sugars <- data.frame(treatment, length)

The basic analysis is easy enough:

anova(lm(length ~ treatment, sugars))

S&R proceed to calculate planned comparisons between control and all
other groups, between gluc+fruct and the remaining sugars, and among
the three pure sugars. I can get the first two comparisons using the
contrasts() function:

contrasts(sugars$treatment) <- matrix(c(-4, 1, 1,  1,  1,
0, -1, 3, -1, -1), 5, 2) 

summary(lm(length ~ treatment, sugars))

The results appear to be the same, excepting that the book calculates
an F value, whereas R produces an equivalent t value. However, I can't
figure out how to calculate a contrast among three groups. For those
without access to Sokal and Rohlf, I've written an R function that
performs the analysis they use, copied below. Is there a better way to
do this in R?

Also, I don't know how to interpret the Estimate and Std. Error
columns from the summary of the lm with contrasts. It's clear the
intercept in this case is the grand mean, but what are the other
values? They appear to be the difference between one of the contrast
groups and the grand mean -- wouldn't an estimate of the differences
between the contrasted groups be more appropriate, or am I (likely)
misinterpreting this?

Thanks!

Tyler

MyContrast <- function(Var, Group, G1, G2, G3=NULL) {
  ## Var == data vector, Group == factor
  ## G1, G2, G3 == character vectors of factor levels to contrast
  
  nG1 = sum(Group %in% G1)
  nG2 = sum(Group %in% G2)
  GrandMean = mean(Var[Group %in% c(G1, G2, G3)]) 
  G1Mean = mean(Var[Group %in% G1])
  G2Mean = mean(Var[Group %in% G2])

  if(is.null(G3))
MScontr = nG1 * ((G1Mean - GrandMean)^2) + 
  nG2 * ((G2Mean - GrandMean)^2)
  else {
  nG3 = sum(Group %in% G3)
  G3Mean = mean(Var[Group %in% G3])
  MScontr = (nG1 * ((G1Mean - GrandMean)^2) + 
 nG2 * ((G2Mean - GrandMean)^2) + 
 nG3 * ((G3Mean - GrandMean)^2))/2
}
  
  An <- anova(lm(Var ~ Group))
  MSwithin = An[3]['Residuals',]
  DegF = An$Df[length(An$Df)]
  Fval = MScontr / MSwithin
  Pval = 1 - pf(Fval, 1, DegF)
  
  return (list(MS_contrasts = MScontr, MS_within = MSwithin, 
   F_value = Fval, P_value = Pval))
}

## The first two contrasts produce the same (+/- rounding error)
## p-values as obtained using contrasts()
MyContrast(sugars$length, sugars$treatment, 'control', 
  c("fructose", "gluc+fruct", "glucose",
  "sucrose"))
MyContrast(sugars$length, sugars$treatment, 'gluc+fruct',
  c("fructose", "glucose", "sucrose"))

## How do you do this in standard R?
MyContrast(sugars$length, sugars$treatment, "fructose", "glucose",
  "sucrose")

__
R-help@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.

[R] Anova(car) SS digits

[Tyler Smith]

Hi,

When I use Anova(car) to produce type III SS, 'Sum Sq' is reported in
integers:

> Anova(bot.lm3, type ="III")
Anova Table (Type III tests)

Response: bottemp
 Sum Sq  DfF value  Pr(>F)
(Intercept)   45295   1 29436.4440 < 2e-16
fungroup  3   2 0.8259 0.44006
numsp.fun11   2 3.4503 0.03460
block   419   468.0017 < 2e-16
fungroup:numsp.fun6   4 1.0317 0.39340
fungroup:block9   8 0.7429 0.65346
numsp.fun:block  12   8 0.9501 0.47795
fungroup:numsp.fun:block 27  16 1.0920 0.36882
Residuals   205 133   

I have tried specifying digits = 7 in the call to Anova, but that
doesn't change anything. The output from another stats program for the
same data confirms that the ouptut above is correct, but is rounded to
the nearest integer. The dataset is too large to post here, but I can
confirm that I get the expected results from the examples in ?Anova,
as well as for type = 'II' on my own data. This output, and relevant
options() and sessionInfo are pasted below. What am I doing wrong?

Thanks,

Tyler

>  mod <- lm(conformity ~ fcategory*partner.status, data=Moore, 
+contrasts=list(fcategory=contr.sum, partner.status=contr.sum))
>  Anova(mod)
Anova Table (Type II tests)

Response: conformity
 Sum Sq Df F value   Pr(>F)
fcategory 11.61  2  0.2770 0.759564
partner.status   212.21  1 10.1207 0.002874
fcategory:partner.status 175.49  2  4.1846 0.022572
Residuals817.76 39 
>  Anova(mod, type="III")
Anova Table (Type III tests)

Response: conformity
 Sum Sq Df  F valuePr(>F)
(Intercept)  5752.8  1 274.3592 < 2.2e-16
fcategory  36.0  2   0.8589  0.431492
partner.status239.6  1  11.4250  0.001657
fcategory:partner.status  175.5  2   4.1846  0.022572
Residuals 817.8 39   

> Anova(bot.lm3, type ="II")
Anova Table (Type II tests)

Response: bottemp
 Sum Sq  Df F value  Pr(>F)
fungroup   1.87   2  0.6069 0.54656
numsp.fun 11.14   2  3.6214 0.02942
block486.47   4 79.0379 < 2e-16
fungroup:numsp.fun 6.32   4  1.0275 0.39553
fungroup:block11.43   8  0.9283 0.49542
numsp.fun:block   11.29   8  0.9169 0.50473
fungroup:numsp.fun:block  26.89  16  1.0920 0.36882
Residuals204.65 133

> options('digits')
$digits
[1] 7

> options('contrasts')
$contrasts
[1] "contr.sum"  "contr.poly"

> sessionInfo()
R version 2.5.1 (2007-06-27) 
i486-pc-linux-gnu 

locale:

LC_CTYPE=en_CA.UTF-8;LC_NUMERIC=C;LC_TIME=en_CA.UTF-8;LC_COLLATE=en_CA.UTF-8;
LC_MONETARY=en_CA.UTF-8;LC_MESSAGES=en_CA.UTF-8;LC_PAPER=en_CA.UTF-8;LC_NAME=C;
LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_CA.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] "stats" "graphics"  "grDevices" "utils" "datasets"  "methods"  
[7] "base" 

other attached packages:
car 
"1.2-7"

__
R-help@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.

Re: [R] anova() or aov()?

[Chuck Cleland]

On 1/12/2009 8:57 AM, j...@in.gr wrote:
> Dear all,
> 
> I have a simple (I think) question that is troubling me lately:
> 
> Is there any main difference between anova() command and aov() command when 
> performing an ANOVA in Experimental design 
> analyses?

  The main difference is that aov() *fits* a single model and anova()
*summarizes* a single model or compares two or more models.

> Thank you for your time,
> 
> Ismini
> 
> __
> R-help@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.

-- 
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

__
R-help@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.

[R] Anova and unbalanced designs

[Skotara]

Dear R-list!

My question is related to an Anova including within and between subject 
factors and unequal group sizes.

Here is a minimal example of what I did:

library(car)
within1 <- c(1,2,3,4,5,6,4,5,3,2); within2 <- c(3,4,3,4,3,4,3,4,5,4)
values <- data.frame(w1 = within1, w2 = within2)
values <- as.matrix(values)
between <- factor(c(rep(1,4), rep(2,6)))
betweenanova <- lm(values ~ between)
with <- expand.grid(within = factor(1:2))
withinanova <- Anova(betweenanova, idata=with, idesign= 
~as.factor(within), type = "III" )


I do not know if this is the appropriate method to deal with unbalanced 
designs.


I observed, that SPSS calculates everything identically except the main 
effect of the within factor, here, the SSQ and F-value are very different
If selecting the option "show means", the means for the levels of the 
within factor in SPSS are the same as:

mean(c(mean(values$w1[1:4]),mean(values$w1[5:10]))) and
mean(c(mean(values$w2[1:4]),mean(values$w2[5:10]))).
In other words, they are calculated as if both groups would have the 
same size.


I wonder if this is a good solution and if so, how could I do the same 
thing in R?
However, I think if this is treated in SPSS as if the group sizes are 
identical,
then why not the interaction, which yields to the same result as using 
Anova()?


Many thanks in advance for your time and help!

__
R-help@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.

[R] ANOVA with subsampling question

[Wade Wall]

Hi all,

I am trying to analyze an experiment I ran, but not sure how to code in R.

I have germinated seeds in petri dishes at 3 different temperatures
(call it low, med, and high) and 2 different light levels (light and
dark).  For each seed I have recorded time to germination (not
counting those that didn't germinate because I will analyze in a
separate ANOVA).  Each temperature/light treatment has 5 petri dishes
with 10 seeds per dish, for a total of 30 dishes and 300 seeds.  The
replicate is petri dish, but I want to treat the seeds as subsampling
in the error effects.

Any help would be appreciated.  I have looked at code for mixed
models, but want to make sure that I am on the right track.

Thanks a lot,

Wade

__
R-help@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.

Re: [R] ANOVA in R

[ONKELINX, Thierry]

Dear Samor,

Note that the R-Sig-mixed-models is more suitable for that kind of questions.
 
You need to add "+", ":" or "*" between the variables in the formula

The correlation is lme() is not the same as the correlation in SAS. The 
correlation in lme() is the correlation among residuals, not among random 
effects. You need one of the pdClasses if you want to define the correlation 
among the random effects. Have a look at ?nlme::pdClasses

Hence your code should probably look like this

Model <- lme(response ~ seq + period + treat*time, random = pdCompSymm(form = 
~1|SUB))
anova(Model)

HTH,

Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and 
Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology 
and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
thierry.onkel...@inbo.be 
www.inbo.be 

To call in the statistician after the experiment is done may be no more than 
asking him to perform a post-mortem examination: he may be able to say what the 
experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure 
that a reasonable answer can be extracted from a given body of data.
~ John Tukey

-Oorspronkelijk bericht-
Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Namens 
Samor Gandhi
Verzonden: woensdag 4 februari 2009 7:43
Aan: r-help@r-project.org
Onderwerp: [R] ANOVA in R

Hi,I'm using a repeated measures ANOVA in R using lme(). The SAS code would be: 
 
PROC MIXED DATA=[data set below];
 CLASS pid treat period time seq;
 MODEL Y = seq period treat time treat*time;
 REPEATED time / SUBJECT=pid TYPE=cs;
RUN,  I donot have SAS, instead I have R and I would like to try the following:
anova(lme(response ~ seq period treat time treat*time,random= ~1|SUB,    
correlation=corCompSymm()))

Is this correct? Can I also write the model as

Y_ijklt = m + a_l + b_k + c_j + d_t + (cd)_jt + u_ijkltY_ijklt is the response 
variable due to pid i, treat j, period k, seq l, and time t. Thank you very 
much in advance for your help :)
Samor 






  
[[alternative HTML version deleted]]


Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.

__
R-help@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.

Re: [R] ANOVA in R

[Samor Gandhi]

Thank you very much.
 
Is the following model would be corresponding to the data?
 
Y_ijklt = m + a_l + b_k + c_j + d_t + (cd)_jt + u_ijklt
 
Cheers,
Samor

--- On Wed, 2/4/09, ONKELINX, Thierry  wrote:

From: ONKELINX, Thierry 
Subject: RE: [R] ANOVA in R
To: samorgan...@yahoo.com, r-help@r-project.org
Cc: r-sig-mixed-mod...@r-project.org
Date: Wednesday, February 4, 2009, 3:21 PM

Dear Samor,

Note that the R-Sig-mixed-models is more suitable for that kind of questions.
 
You need to add "+", ":" or "*" between the
variables in the formula

The correlation is lme() is not the same as the correlation in SAS. The
correlation in lme() is the correlation among residuals, not among random
effects. You need one of the pdClasses if you want to define the correlation
among the random effects. Have a look at ?nlme::pdClasses

Hence your code should probably look like this

Model <- lme(response ~ seq + period + treat*time, random = pdCompSymm(form
= ~1|SUB))
anova(Model)

HTH,

Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology
and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
thierry.onkel...@inbo.be 
www.inbo.be 

To call in the statistician after the experiment is done may be no more than
asking him to perform a post-mortem examination: he may be able to say what the
experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure
that a reasonable answer can be extracted from a given body of data.
~ John Tukey

-Oorspronkelijk bericht-
Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Namens
Samor Gandhi
Verzonden: woensdag 4 februari 2009 7:43
Aan: r-help@r-project.org
Onderwerp: [R] ANOVA in R

Hi,I'm using a repeated measures ANOVA in R using lme(). The SAS code would
be:  
PROC MIXED DATA=[data set below];
 CLASS pid treat period time seq;
 MODEL Y = seq period treat time treat*time;
 REPEATED time / SUBJECT=pid TYPE=cs;
RUN,  I donot have SAS, instead I have R and I would like to try the
following:
anova(lme(response ~ seq period treat time treat*time,random= ~1|SUB,
   correlation=corCompSymm()))

Is this correct? Can I also write the model as

Y_ijklt = m + a_l + b_k + c_j + d_t + (cd)_jt + u_ijkltY_ijklt is the response
variable due to pid i, treat j, period k, seq l, and time t. Thank you very much
in advance for your help :)
Samor 






  
[[alternative HTML version deleted]]


Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd
is
door een geldig ondertekend document. The views expressed in  this message
and any annex are purely those of the writer and may not be regarded as stating

an official position of INBO, as long as the message is not confirmed by a duly

signed document.



  
[[alternative HTML version deleted]]

__
R-help@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.

[R] anova p value extraction

[Paul Benton]

hello all,

Quick question, how do I get the p value out of the anova?

Thanks,

Paul

> pb<-aov(as.numeric(diff[5,16:33]) ~ grF)
> summary(pb)
Df Sum SqMean Sq F value  Pr(>F)
grF  3 2.7860e+10 9.2867e+09  4.2236 0.02534 *
Residuals   14 3.0783e+10 2.1988e+09
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> str(summary(pb))
List of 1
 $ :Classes 'anova' and 'data.frame':   2 obs. of  5 variables:
  ..$ Df : num [1:2] 3 14
  ..$ Sum Sq : num [1:2] 2.79e+10 3.08e+10
  ..$ Mean Sq: num [1:2] 9.29e+09 2.20e+09
  ..$ F value: num [1:2] 4.22   NA
  ..$ Pr(>F) : num [1:2] 0.0253 NA
 - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> attr(summary(pb), "Pr(>F)")
NULL

__
R-help@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.

[R] ANOVA between linear models.

[Richard Martin]

Hi All,

I'm accustomed to performing an ANOVA to aid in choosing between
linear models (for example y~x or y~x+x^2), however with different
models I can't seem to do it. I'm trying to fit an exponential model
of the form ye^(bt).

Below is a code snippet that highlights what I'm trying to do

s = rnorm(100, sd = 1, mean=10)
s = s + seq(0.1,10,0.1)
x = 1:length(s)
model.lin = lm(s ~ x)
model.poly = lm(s ~x  + I(x^2))
model.exp = lm(log(s) ~ x)

anova(model.lin, model.poly)
#gives the correct outcome

anova(model.lin, model.exp)
#doesn't.

This fails because of the transformation of the response variable. Can
anyone give any advice as to how I should proceed - is there a
different test to use in this instance, or some way of reversing the
transform after the model has been fitted?

Any help greatly appreciated!!

Richard

-- 
Contendere, Explorare, Invenire, et non Cedere

__
R-help@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.

Re: [R] ANOVA and lmer

[Douglas Bates]

On Feb 7, 2008 3:43 PM, Eric Imbert <[EMAIL PROTECTED]> wrote:

> I am analyzing from a very simple experiment.
> I have measured plants of two different colours (yellow and purple) in 9
> different populations.

> So, I have two different factors : a fixed effect (Colour with two
> levels) and a random one (Population with 9 levels).

> I first analyzed the data with the aov function
> LargS is the variable

> aov(formula = LargS ~ Col + Error(Col/Pop))

> Terms:
> Col
> Sum of Squares  3.440351
> Deg. of Freedom1

> Estimated effects are balanced

> Stratum 2: Col:Pop

> Terms:
> Residuals
> Sum  of Squares   3017.112
> Deg. of Freedom16

> Residual standard error: 13.73206

> Stratum 3: Within

> Terms:
> Residuals
> Sum of Squares   3347.385
> Deg. of Freedom   302

> To test for the interaction Col*Pop, I used the following F-ratio =
> (3017/16)/(3347/302) = 188. Highly significant !

> Now, let's go to the analysis performed by lmer - First I do the linear
> model without the Col*Pop interaction :
> m3=lmer(LargS ~ Col + (1 | Pop)

> And next with the interaction : m2=lmer(LargS ~ Col + (Col | Pop))

I don't think this model provides a comparison that is directly
comparable with the test you describe above.

There is a model that is intermediate to your m2 and m3 models in
terms of the number of parameters.  (Actually when I see these results
I realize that there is a bug in anova for lmer models in that the
count of the parameters 1 less than it should be.  I forgot to include
the residual variance as a parameter to be estimated.  However, this
will not affect comparisons here because it is only the differences in
the counts that matter.)

You model m2 is shown as having 5 parameters, which should be 6 (2
fixed effects, 3 variance components for the random effects, 1
residual variance estimate), while m2 is shown as having 3 parameters
(should be 4: 2 f.e. + 1 r.e. var + 1 resid var).  The model

LargS ~ Col + (1|Pop) + (1|Col:Pop)

which allows for random effects for each population and for each
Col:Pop combination, has 5 parameters (currently shown as 4).  The
random effects for Pop are independent of each other and with a
constant variance while the random effects for the Col:Pop
combinations are independent of each other and of the Pop random
effects and with another constant variance.  Thus there are 2 f.e. + 2
r.e. var + 1 resid. var.

I think it is important to plot the data first and decide which of
these models is indicated by the data.  It sounds like the structure
of your experiment is similar to the structure of the "barley" data
set in the lattice package.  Deepayan Sarkar's forthcoming book
"Lattice: Multivariate Data Visualization with R" (to be published by
Springer and due to ship next month) describes very effective ways of
plotting such data.  You can get a preview at his web site for the
book http://lmdvr.r-forge.r-project.org

> Comparing both models : anova(m2,m3) :
>
>Df AIC BIC  logLik  Chisq Chi Df Pr(>Chisq)
> m3  3 1710.67 1721.97 -852.33
> m2  5 1714.59 1733.43 -852.30 0.0746  2 0.9634
>
> => Conclusion : the interaction Col*Pop is not significant !
>
> I guess I am missing something.
> Who can help ?
>
> Eric
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@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.
>

__
R-help@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.

Re: [R] anova power calculations

[Will Holcomb]

To answer my own question, the notes that I was drawing from made it seem as
through the mean I was using for my sum of squares was unweighted. In
reality, I needed to weight the contribution of each contributing mean by
the sample size. Once I did that, I got a value that agreed with SAS.

So, the solution for computing power of an unbalanced ANOVA is:

within.var = 9 ^ 2
means = c(17.5, 19, 25, 20.5)
J = length(means)

# Quantile of the cutoff point in the central F
central.quart = qf(.05, J - 1, N - J, lower.tail = FALSE)

weighted.means = data.frame(Mean = means, n = c(10, 10, 50, 10))
N = sum(weighted.means$n)

weighted.mean = weighted.mean(weighted.means$Mean, weighted.means$n)

# Noncentrality parameter for unbalanced ANOVA
noncentral.param = sum(weighted.means$n * (weighted.means$Mean -
weighted.mean) ^ 2) / within.var

# Probability of central quantile in noncentral distribution
noncentral.p = pf(central.quart, J - 1, N - J, noncentral.param,
lower.tail= FALSE)

Will

2008/2/21 Will Holcomb <[EMAIL PROTECTED]>:

> I sent a message a couple days ago about doing calculations for power of
> the ANOVA. Several people got back to me very quickly which I really
> appreciated.
>
> I'm working now on a similar problem, but instead of a balanced ANOVA, I
> have an unbalanced one. The first part of the question was:
>
> You assume that the within-population standard deviations all equal 9. You
> set the Type 1 error rate at á = .05. You presume that the population means
> will have the following values: uA = 17.5, uB = 19, uC = 25, and uD = 20.5.
> You intend to run 80 subjects in all, with equal n's across all 4 groups.
> You plan on conducting a one-way ANOVA. Compute your power to reject the
> null hypothesis under these conditions.
>
> I did:
>
> within.var = 9 ^ 2
> means = c(17.5, 19, 25, 20.5)
> N = 80
> J = length(means)
> power.anova.test(groups = J, n = N / J,
>  between.var = var(means),
>  within.var = within.var,
>  sig.level = 0.05)
>
> This gives me 0.6155 which agrees with SAS. The next problem though is:
>
> You have the same Type 1 error rate and make the same assumptions about
> the population standard deviation and the population means as in part a. You
> still have 80 subjects in all but now you want to know how power might
> change by running 10 subjects in groups A, B, and D and 50 subjects in group
> C. Determine the power under this subject allocation scheme.
>
> For this one I am doing:
>
> # Quantile of the cutoff point in the central F
> central.quant = qf(.05, J - 1, N - J, lower.tail = FALSE)
> weighted.means = data.frame(Mean = means, Weight = c(10, 10, 50, 10))
> # Noncentrality parameter for unbalanced ANOVA
> noncentral.param = 0
> for(i in 1:length(weighted.means$Mean)) {
>   noncentral.param = (noncentral.param + weighted.means$Weight[i] *
>   (weighted.means$Mean[i] - mean(weighted.means$Mean))
> ^ 2)
> }
> noncentral.param = noncentral.param / within.var
> # Probability of central quantile in noncentral distribution
> noncentral.p = pf(central.quant, J - 1, N - J, noncentral.param,
> lower.tail = FALSE)
> noncentral.p
>
> The logic behind this is in my assignment at:
>
> http://odin.himinbi.org/classes/psy304b/homework_2.xhtml#p2b
>
> This works for a balanced ANOVA and gives the same result as
> power.anova.test (and SAS). For the unbalanced ANOVA though it is giving
> me a different result though than SAS, 0.8759455 versus 0.680.
>
> So is there a straightforward way to compute the power of an unbalanced
> ANOVA? If there isn't, does anyone have any idea what is wrong with my code?
> The SAS I am comparing it to is:
>
> Data Dep;
> Input cue $ mean uneven_weight;
> datalines;
> A 17.5  1
> B 191
> C 255
> D 20.5  1
> ;
>
> proc glmpower;
> class cue;
> model mean = cue;
> weight uneven_weight;
> power
> stddev = 9
> alpha = 0.05
> ntotal= 80
> power = .;
> run;
>
> Any help would be much appreciated.
>
> Will
>

[[alternative HTML version deleted]]

__
R-help@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.

[R] Anova function and glm.nb

[Nelson, Gary (FWE)]

Hi All,

I am using the glm.nb function from the MASS package (current version)
to fit and compare GLMs with negative binomial error distributions.  My
question is: what is the appropriate method to use in the anova function
to compare models? If only one fitted object like

m<-glm.nb(number<-p+sal+temp,data=data)

is specified in the anova function (anova(m)), a fixed theta is used to
generate the analysis of deviance:

Analysis of Deviance Table
Model: Negative Binomial(0.2345), link: log
Response: number
Terms added sequentially (first to last)

  Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL117122.707  
p  1   11.327   116111.380 0.001
sal12.286   115109.094 0.131
tem11.979   114107.115 0.159
ph 12.567   113104.549 0.109
Warning message:
In anova.negbin(m) : tests made without re-estimating 'theta'


If multiple fitted objects like

m1<-glm.nb(number~1,data=data)
m2<-glm.nb(number~p,data=data)
m3<-glm.nb(number~p+sal,data=data)
m4<-glm.nb(number~p+sal+temp,data=data)


is specified (anova(m1,m2,m3,4)), the theta is assumed re-estimated in
each case to calculate the likelihood ratio tests:

Likelihood ratio tests of Negative Binomial Models
Response: number
   Model theta Resid. df2 x log-lik.   Testdf LR
stat. Pr(Chi)
1  1 0.1892056   117   -527.7463

2  p 0.2153105   116   -517.9349 1 vs 2 1
9.811382 0.001734351
3p + sal 0.2214626   115   -515.7942 2 vs 3 1
2.140706 0.143435894
4  p + sal + tem 0.2261900   114   -513.8846 3 vs 4 1
1.909643 0.167002884
5 p + sal + tem + ph 0.2344827   113   -511.3633 4 vs 5 1
2.521237 0.112322429


The conclusions are the same, but I'd like to know if one method is
favored over the other. 

Thanks,

Gary Nelson.

__
R-help@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.

[R] anova with ordered groups

[Michele Pinelli]

Hi,
I have to do a comparison among three groups of genetic transcription
levels.

I have a situation like this:
group 0: baseline
group 1: first treatment
group 2: second treatment

In the first group, I have only 2 samples, in the second one 4 samples and
in the last group I have 10 samples.

I would check if the trnascription of a gene increases from the baseline
situation to the fist treatment AND from the first treatment to the second
treatment.

I did a table like this:
GENE-LEVEL, GROUP
10,0
12,0
22,1
23,1
21,1
24,1
31,2
32,2
33,2
32,2
31,2
32,2
33,2
31,2
32,2
33,2

and I did this analysis:
lm(GENE-LEVEL~GROUP)
using the GROUP as integer and not as level.

Is it correct? and if the relationship is not strictly linear, ie the second
treatment level is not the double of the first treatment, is this analysis
still good?
Mainly I would only check whether group0group1>group2)

Thank you
Michele Pinelli


Michele Pinelli, MD, PhD
Complex Disease Genetics Unit
Department of Cellular and Molecular Biology and Pathology "L. Califano"
University "Federico II", Naples
Italy
[EMAIL PROTECTED]
Skype: dr.sson (Michele Pinelli/Naples, IT)


"People still respect science and technology. … There is a problem however.
They have no idea what science is."
by Alan Leshner,
CEO of the American Association for the Advancement of Science

[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA and TukeyHSD disagrees?

[Fredrik Karlsson]

Dear list,

Sorry for posting a borderline statistical question on the list, but hte
SPSS people around me just stares at me blankly when refering to tests with
any term other than ANOVA and post-hoc. I would appreciate any insight on
how this all is possible:

I have a model fitted by aov() stored in "ppdur", which gives this result
when using ANOVA:

> anova(ppdur)
Analysis of Variance Table

Response: PAPositionPercentOfVoweldur
  Df  Sum Sq Mean Sq F value  Pr(>F)
UtteranceType   4   247316183  2.7642 0.02696 *
SyllLable   1   14584   14584  6.5202 0.01094 *
Cycle   1 798 798  0.3566 0.55067
Speaker 299754987  2.2297 0.10855
Label   120082008  0.8979 0.34377
UtteranceType:SyllLable 4   152103803  1.7001 0.14854
UtteranceType:Cycle 4   131923298  1.4745 0.20855
SyllLable:Cycle 1   11306   11306  5.0545 0.02497 *
UtteranceType:Speaker   7   137211960  0.8764 0.52488
SyllLable:Speaker   21291 645  0.2885 0.74951
Cycle:Speaker   2   107535377  2.4038 0.09135 .
UtteranceType:Label 43579 895  0.4000 0.80871
SyllLable:Label 144994499  2.0114 0.15670
Cycle:Label 1 229 229  0.1022 0.74929
Speaker:Label   21241 620  0.2774 0.75788
UtteranceType:SyllLable:Cycle   3 473 158  0.0705 0.97571
UtteranceType:SyllLable:Speaker 6   139192320  1.0372 0.40006
UtteranceType:Cycle:Speaker 31221 407  0.1820 0.90865
SyllLable:Cycle:Speaker 21457 729  0.3258 0.72210
UtteranceType:SyllLable:Label   238231911  0.8545 0.42607
UtteranceType:Cycle:Label   385662855  1.2766 0.28160
SyllLable:Cycle:Label   135753575  1.5983 0.20669
UtteranceType:Speaker:Label 42658 664  0.2970 0.87990
SyllLable:Speaker:Label 2 139  70  0.0311 0.96938
Cycle:Speaker:Label 2   135996800  3.0400 0.04866 *
UtteranceType:SyllLable:Cycle:Speaker   220151008  0.4505 0.63757
UtteranceType:SyllLable:Cycle:Label 1  11  11  0.0051 0.94328
UtteranceType:SyllLable:Speaker:Label   1 603 603  0.2695 0.60386
Residuals 539 12056052237
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ok now, when I want to know where the differences are, I get this result
from TukeyHSD:


> TukeyHSD(ppdur,c("Cycle:Speaker:Label" ),ordered=TRUE)
  Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered

Fit: aov(formula = PAPositionPercentOfVoweldur ~ UtteranceType * SyllLable *
Cycle * Speaker * Label, data =PABTSub)

$`Cycle:Speaker:Label`
difflwr   upr p adj
3:Andrea:!H*L-1:Lavinia:!H*L   2.6069140 -37.499300  42.71313 1.000
1:Vito:!H*L-1:Lavinia:!H*L 8.7764075 -85.090794 102.64361 1.000
3:Andrea:H*L-1:Lavinia:!H*L   12.3411960 -18.883688  43.56608 0.9792814
1:Vito:H*L-1:Lavinia:!H*L 15.0416018 -32.746962  62.83017 0.9968844
1:Andrea:H*L-1:Lavinia:!H*L   15.2934976 -17.987036  48.57403 0.9379977
1:Lavinia:H*L-1:Lavinia:!H*L  16.9297832 -14.670124  48.52969 0.8394874
3:Lavinia:H*L-1:Lavinia:!H*L  18.3218965 -13.445765  50.08956 0.7631935
3:Lavinia:!H*L-1:Lavinia:!H*L 20.9338365 -19.932636  61.80031 0.8761167
3:Vito:!H*L-1:Lavinia:!H*L24.3874104 -18.890036  67.66486 0.7894161
3:Vito:H*L-1:Lavinia:!H*L 27.8865684  -6.758302  62.53144 0.2589397
1:Andrea:!H*L-1:Lavinia:!H*L  28.8093072 -18.979256  76.59787 0.7077134
1:Vito:!H*L-3:Andrea:!H*L  6.1694934 -87.982875 100.32186 1.000
3:Andrea:H*L-3:Andrea:!H*L 9.7342820 -22.337674  41.80624 0.9977466
1:Vito:H*L-3:Andrea:!H*L  12.4346877 -35.911602  60.78098 0.9995138
1:Andrea:H*L-3:Andrea:!H*L12.6865836 -21.389960  46.76313 0.9870844
1:Lavinia:H*L-3:Andrea:!H*L   14.3228692 -18.114318  46.76006 0.9529437
3:Lavinia:H*L-3:Andrea:!H*L   15.7149825 -16.885651  48.31562 0.9149770
3:Lavinia:!H*L-3:Andrea:!H*L  18.3269225 -23.190369  59.84421 0.9530401
3:Vito:!H*L-3:Andrea:!H*L 21.7804964 -22.112036  65.67303 0.8978922
3:Vito:H*L-3:Andrea:!H*L  25.2796544 -10.130570  60.68988 0.4470040
1:Andrea:!H*L-3:Andrea:!H*L   26.2023932 -22.143897  74.54868 0.821
3:Andrea:H*L-1:Vito:!H*L   3.5647885 -87.160916  94.29049 1.000
1:Vito:H*L-1:Vito:!H*L 6.2651943 -91.407252 103.93764 1.000
1:Andrea:H*L-1:Vito:!H*L   6.5170902 -84.936471  97.97065 1.000
1:Lavinia:H*L-1:Vito:!H*L  8.1533757 -82.702082  99.00883 1.000
3:Lavinia:

[R] Anova interaction not tested

[Gabriel Murray]

I've noticed with certain datasets that when I try to do an anova and test
for main effects and interaction for two explanatory variables, sometimes
the main effect results are given but not the interaction results. For
example,

ex1 = aov(Score ~ var1*var2, data=myData)
summary(ex1)

gives me only the main effects for var1 and var2, but not the interaction. I
also tried

ex1 = aov(Score ~ var1+var2+var1:var2, data=myData)
summary(ex1)

but it still only gives the main effects. The only way I can get the
interaction results is if I test for ONLY that, e.g.

ex1 = aov(Score ~ var1:var2, data=myData)
summary(ex1)

Why would the interaction not be tested in the first two examples?

Thanks,
Gabe

[[alternative HTML version deleted]]

__
R-help@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.

[R] anova table with lmer

[Julien Beguin]

Hi everyone,

I am wandering if someone was able to get an anova table using the function lmer
in lme4 package with an other distribution than gaussian (ex: poisson). In fact,
I am using a GLMM with a poisson distribution where I have 2 independant fixed
variables (coded as factors with respectively two and three levels) and I am
interested with simple and interaction effects (I also have a couple of random
variables). However, rather than to have an estimate of each combination (plus
intercept), I would like first to know the effects of my 3 sources of variation
(not the combination of each level). This works well with a gaussian
distribution (using anova()) but not with another one (ex: poisson). I have
seen exemple where it seems to work pretty well even with factor variables so I
guess there is a way to do it... If someone know how to do, I would greatly
appreciate your help.

Best regards,

Julien

--
Julien Beguin
Etudiant au doctorat
Faculté de Foresterie et de Géomatique
Université Laval, local 2113
2405 rue de la Terrasse, G1V 0A6 Québec (Qc)
Tel: (418) 656-2131 poste 2620
http://www.cen.ulaval.ca/anticosti/jbeguin.html

__
R-help@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.

Re: [R] ANOVA/statistics question

[Tal Galili]

Hi Douglas
I would go for a different command then aov.
something like:
?cor
or
?cor.test
To also get the p value of the correlation.

Cheers,
Tal



On Sat, Apr 25, 2009 at 8:27 AM, drmh wrote:

>
> (Have searched for this already)
>
> Hi,
>
> How do you find the strength of correlation between two variables using an
> ANOVA table?  "Pr(>F)" gives the statistical significance of the
> association, but not the strength of the correlation.
>
> See data (from R) below
>
> Readable:
>  "Df" "Sum Sq""Mean Sq"   "F value"
> "Pr(>F)"
> "pre"  1   0.00593  0.00593936
> 0.7450563   0.401636958677004
> "coh" 1   0.04311  0.04311302
>  5.4082639
> 0.0344751749542619
> "Residuals"15 0.11957  0.00797169  NA
> NA
>
> Original:
> "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
> "pre" 1 0.0059393604629317 0.0059393604629317 0.745056336657567
> 0.401636958677004
> "coh" 1 0.0431130207164516 0.0431130207164516 5.40826398359156
> 0.0344751749542619
> "Residuals" 15 0.119575396598395 0.00797169310655964 NA NA
>
> Any help would be greatly appreciated,
> Douglas Holmes
> --
> View this message in context:
> http://www.nabble.com/ANOVA-statistics-question-tp23231563p23231563.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@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.
>



-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
http://www.r-statistics.com/
http://www.talgalili.com
http://www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] ANOVA/statistics question

[drmh]

Hi, thanks for your prompt reply

In my situation, the dependent variable is "post-test" and the independent
variables are "pre" and "coh".
Howw would I find the correlation between coh and post with the effect of
"pre" regressed using your commands?



Tal Galili wrote:
> 
> Hi Douglas
> I would go for a different command then aov.
> something like:
> ?cor
> or
> ?cor.test
> To also get the p value of the correlation.
> 
> Cheers,
> Tal
> 
> 
> 
> On Sat, Apr 25, 2009 at 8:27 AM, drmh
> wrote:
> 
>>
>> (Have searched for this already)
>>
>> Hi,
>>
>> How do you find the strength of correlation between two variables using
>> an
>> ANOVA table?  "Pr(>F)" gives the statistical significance of the
>> association, but not the strength of the correlation.
>>
>> See data (from R) below
>>
>> Readable:
>>  "Df" "Sum Sq""Mean Sq"   "F
>> value"
>> "Pr(>F)"
>> "pre"  1   0.00593  0.00593936
>> 0.7450563   0.401636958677004
>> "coh" 1   0.04311  0.04311302
>>  5.4082639
>> 0.0344751749542619
>> "Residuals"15 0.11957  0.00797169  NA
>> NA
>>
>> Original:
>> "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
>> "pre" 1 0.0059393604629317 0.0059393604629317 0.745056336657567
>> 0.401636958677004
>> "coh" 1 0.0431130207164516 0.0431130207164516 5.40826398359156
>> 0.0344751749542619
>> "Residuals" 15 0.119575396598395 0.00797169310655964 NA NA
>>
>> Any help would be greatly appreciated,
>> Douglas Holmes
>> --
>> View this message in context:
>> http://www.nabble.com/ANOVA-statistics-question-tp23231563p23231563.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@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.
>>
> 
> 
> 
> -- 
> --
> 
> 
> My contact information:
> Tal Galili
> Phone number: 972-50-3373767
> FaceBook: Tal Galili
> My Blogs:
> http://www.r-statistics.com/
> http://www.talgalili.com
> http://www.biostatistics.co.il
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/ANOVA-statistics-question-tp23231563p23234421.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.

Re: [R] ANOVA/statistics question

[Tal Galili]

Hi Douglas.
So you want to check for correlation or regression ?

how many levels does "pre" have ?
you could subset the variables you want to check correlation on, by the pre
levels.

for example:
Let's say pre has two levels: 1 and 2. then you can do:
cor(y[pre == 1], x[pre == 1])
cor(y[pre == 2], x[pre == 2])


Also, If you want to go for regression, you can go to something like:
summary(lm(y ~ x + pre))
But I suggest getting much more understanding of what linear regression is
before using it's methods...

Good luck,
Tal



On Sat, Apr 25, 2009 at 1:26 PM, drmh wrote:

>
> Hi, thanks for your prompt reply
>
> In my situation, the dependent variable is "post-test" and the independent
> variables are "pre" and "coh".
> Howw would I find the correlation between coh and post with the effect of
> "pre" regressed using your commands?
>
>
>
> Tal Galili wrote:
> >
> > Hi Douglas
> > I would go for a different command then aov.
> > something like:
> > ?cor
> > or
> > ?cor.test
> > To also get the p value of the correlation.
> >
> > Cheers,
> > Tal
> >
> >
> >
> > On Sat, Apr 25, 2009 at 8:27 AM, drmh
> > wrote:
> >
> >>
> >> (Have searched for this already)
> >>
> >> Hi,
> >>
> >> How do you find the strength of correlation between two variables using
> >> an
> >> ANOVA table?  "Pr(>F)" gives the statistical significance of the
> >> association, but not the strength of the correlation.
> >>
> >> See data (from R) below
> >>
> >> Readable:
> >>  "Df" "Sum Sq""Mean Sq"   "F
> >> value"
> >> "Pr(>F)"
> >> "pre"  1   0.00593  0.00593936
> >> 0.7450563   0.401636958677004
> >> "coh" 1   0.04311  0.04311302
> >>  5.4082639
> >> 0.0344751749542619
> >> "Residuals"15 0.11957  0.00797169  NA
> >> NA
> >>
> >> Original:
> >> "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
> >> "pre" 1 0.0059393604629317 0.0059393604629317 0.745056336657567
> >> 0.401636958677004
> >> "coh" 1 0.0431130207164516 0.0431130207164516 5.40826398359156
> >> 0.0344751749542619
> >> "Residuals" 15 0.119575396598395 0.00797169310655964 NA NA
> >>
> >> Any help would be greatly appreciated,
> >> Douglas Holmes
> >> --
> >> View this message in context:
> >>
> http://www.nabble.com/ANOVA-statistics-question-tp23231563p23231563.html
> >> Sent from the R help mailing list archive at Nabble.com.
> >>
> >> __
> >> R-help@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.
> >>
> >
> >
> >
> > --
> > --
> >
> >
> > My contact information:
> > Tal Galili
> > Phone number: 972-50-3373767
> > FaceBook: Tal Galili
> > My Blogs:
> > http://www.r-statistics.com/
> > http://www.talgalili.com
> > http://www.biostatistics.co.il
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@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.
> >
> >
>
> --
> View this message in context:
> http://www.nabble.com/ANOVA-statistics-question-tp23231563p23234421.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@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.
>



-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
http://www.r-statistics.com/
http://www.talgalili.com
http://www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] ANOVA/statistics question

[drmh]

Hello again,
In my situation, I have three variables: pretest, posttest, and cohesion. 

I want to work out the correlation between postest and cohesion. 

I looked at multiple sets of data and created ANOVA tables of them. However,
as pretest and postest are sometimes correlated (with a statistical
significance < 0.05), it is necessary to discount the effect of pretest to
work out the real correlation of posttest and coherence.. I need a system
for working out the strength of the correlation between posttest and
coherence, when does actually occur. 

According to my understanding level refers the amount or magnitude of
experimental units. Pretest, posttest are scores - range from any value from
0 to 1. Cohesion can be any value.

What exactly would
cor(y[pre == 1], x[pre == 1])
cor(y[pre == 2], x[pre == 2])
give me?

I understand that my lack of understanding may be exasperating, but it is
not for want of effort - I have put in hours trying to understand this
stuff.

Thanks,
Douglas Holmes



Tal Galili wrote:
> 
> Hi Douglas.
> So you want to check for correlation or regression ?
> 
> how many levels does "pre" have ?
> you could subset the variables you want to check correlation on, by the
> pre
> levels.
> 
> for example:
> Let's say pre has two levels: 1 and 2. then you can do:
> cor(y[pre == 1], x[pre == 1])
> cor(y[pre == 2], x[pre == 2])
> 
> 
> Also, If you want to go for regression, you can go to something like:
> summary(lm(y ~ x + pre))
> But I suggest getting much more understanding of what linear regression is
> before using it's methods...
> 
> Good luck,
> Tal
> 
> 
> 
> On Sat, Apr 25, 2009 at 1:26 PM, drmh
> wrote:
> 
>>
>> Hi, thanks for your prompt reply
>>
>> In my situation, the dependent variable is "post-test" and the
>> independent
>> variables are "pre" and "coh".
>> Howw would I find the correlation between coh and post with the effect of
>> "pre" regressed using your commands?
>>
>>
>>
>> Tal Galili wrote:
>> >
>> > Hi Douglas
>> > I would go for a different command then aov.
>> > something like:
>> > ?cor
>> > or
>> > ?cor.test
>> > To also get the p value of the correlation.
>> >
>> > Cheers,
>> > Tal
>> >
>> >
>> >
>> > On Sat, Apr 25, 2009 at 8:27 AM, drmh
>> > wrote:
>> >
>> >>
>> >> (Have searched for this already)
>> >>
>> >> Hi,
>> >>
>> >> How do you find the strength of correlation between two variables
>> using
>> >> an
>> >> ANOVA table?  "Pr(>F)" gives the statistical significance of the
>> >> association, but not the strength of the correlation.
>> >>
>> >> See data (from R) below
>> >>
>> >> Readable:
>> >>  "Df" "Sum Sq""Mean Sq"   "F
>> >> value"
>> >> "Pr(>F)"
>> >> "pre"  1   0.00593  0.00593936
>> >> 0.7450563   0.401636958677004
>> >> "coh" 1   0.04311  0.04311302
>> >>  5.4082639
>> >> 0.0344751749542619
>> >> "Residuals"15 0.11957  0.00797169  NA
>> >> NA
>> >>
>> >> Original:
>> >> "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
>> >> "pre" 1 0.0059393604629317 0.0059393604629317 0.745056336657567
>> >> 0.401636958677004
>> >> "coh" 1 0.0431130207164516 0.0431130207164516 5.40826398359156
>> >> 0.0344751749542619
>> >> "Residuals" 15 0.119575396598395 0.00797169310655964 NA NA
>> >>
>> >> Any help would be greatly appreciated,
>> >> Douglas Holmes
>> >> --
>> >> View this message in context:
>> >>
>> http://www.nabble.com/ANOVA-statistics-question-tp23231563p23231563.html
>> >> Sent from the R help mailing list archive at Nabble.com.
>> >>
>> >> __
>> >> R-help@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.
>> >>
>> >
>> >
>> >
>> > --
>> > --
>> >
>> >
>> > My contact information:
>> > Tal Galili
>> > Phone number: 972-50-3373767
>> > FaceBook: Tal Galili
>> > My Blogs:
>> > http://www.r-statistics.com/
>> > http://www.talgalili.com
>> > http://www.biostatistics.co.il
>> >
>> >   [[alternative HTML version deleted]]
>> >
>> > __
>> > R-help@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.
>> >
>> >
>>
>> --
>> View this message in context:
>> http://www.nabble.com/ANOVA-statistics-question-tp23231563p23234421.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@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, se

Re: [R] ANOVA/statistics question

[Peter Flom]

drmh  wrote
>
>Hello again,
>In my situation, I have three variables: pretest, posttest, and cohesion. 
>
>I want to work out the correlation between postest and cohesion. 
>

cor(cohesion, posttest) gives you this.

>I looked at multiple sets of data and created ANOVA tables of them. However,
>as pretest and postest are sometimes correlated (with a statistical
>significance < 0.05), it is necessary to discount the effect of pretest to
>work out the real correlation of posttest and coherence.. I need a system
>for working out the strength of the correlation between posttest and
>coherence, when does actually occur. 

Whether pretest and posttest are correlated, and whether that correlation is 
statistically significant, is irrelevant to your question as posed.  Correlation
is defined between two variables, not among three.  

You might want some sort of regression such as 

lm(cohesion~pretest+posttest)

but you might not

>
>According to my understanding level refers the amount or magnitude of
>experimental units.


What is level?  You mention pretest, posttest and cohesion - now you mention 
level.
What are these experimental units?

 Pretest, posttest are scores - range from any value from
>0 to 1. Cohesion can be any value.
>
>What exactly would
>cor(y[pre == 1], x[pre == 1])
>cor(y[pre == 2], x[pre == 2])
>give me?
>

well, you said above that pretest and posttest can range from 0 to 1; if this 
is the case, pre would rarely be 1 and never be 2, so the first line above 
wouldn't give you much, and the second wouldn't give you anything.  Also, you 
are now using y and x instead of (presumably) cohesion and posttest, and pre 
instead of, presumably, pretest.


Peter

Peter L. Flom, PhD
Statistical Consultant
www DOT peterflomconsulting DOT com

__
R-help@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.

Re: [R] anova(cph(..) output

[Frank E Harrell Jr]

pompon wrote:

Hello,

I am a beginner in R and statistics, so my question may be trivial. Sorry in
advance.
I performed a Cox proportion hazard regression with 2 categorical variables
with cph{design}. Then an anova on the results.
the output is 


anova(cph(surv(survival, censor) ~ plant + leaf.age + plant*leaf.age,
Mpnymph)


Wald Statistics  Response: Surv(survival, censored) 

 FactorChi-Square d.f. P 
 plant  (Factor+Higher Order Factors) 96.96 12   <.0001
  All Interactions   10.58 
6   0.1022

 leaf.age  (Factor+Higher Order Factors)  29.11  7   0.0001
  All Interactions 10.58 
6   0.1022

 plant * leaf.age  (Factor+Higher Order Factors)  10.58  6   0.1022
 TOTAL   106.63 13   <.0001

What do "All interaction" stand for?
The real df of for plant is 6 and 1 for leaf.age. Then, which chi square is
one for my main factors anf their interaction.

thank you,
Julien.


Julien,

I know what you mean when you say 'real df' but that's not the whole 
story as plant has 6 more df by interacting with a single df variable. 
There is no such thing as 'the' main effect test for plant.  The 12 df 
test is unique and tests whether plant is associated with Y for any 
level of leaf.age.


You can see exactly what is being tested by using various print options 
for anova.Design, as described in the help file.  The "dots" option is 
easy on the eyes.


Frank
--
Frank E Harrell Jr   Professor and Chair   School of Medicine
 Department of Biostatistics   Vanderbilt University

__
R-help@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.

Re: [R] anova(cph(..) output

[pompon]

Hi,

Thank you very much for the answer.

However, I have still some misunderstandings.
from the output, can we say that plant and leaf age are significant but not
their interaction?
And the last question I promise, what would you advise me to write in the
paper to explain the different method and ackonwledge for the df?

Thank you again,
julien.

 

Frank E Harrell Jr wrote:
> 
> pompon wrote:
>> Hello,
>> 
>> I am a beginner in R and statistics, so my question may be trivial. Sorry
>> in
>> advance.
>> I performed a Cox proportion hazard regression with 2 categorical
>> variables
>> with cph{design}. Then an anova on the results.
>> the output is 
>> 
>>> anova(cph(surv(survival, censor) ~ plant + leaf.age + plant*leaf.age,
>>> Mpnymph)
>> 
>> Wald Statistics  Response: Surv(survival,
>> censored) 
>> 
>>  FactorChi-Square
>> d.f. P 
>>  plant  (Factor+Higher Order Factors) 96.96 12   <.0001
>>   All Interactions   10.58 
>> 6   0.1022
>>  leaf.age  (Factor+Higher Order Factors)  29.11  7   0.0001
>>   All Interactions 10.58 
>> 6   0.1022
>>  plant * leaf.age  (Factor+Higher Order Factors)  10.58  6   0.1022
>>  TOTAL   106.63 13   <.0001
>> 
>> What do "All interaction" stand for?
>> The real df of for plant is 6 and 1 for leaf.age. Then, which chi square
>> is
>> one for my main factors anf their interaction.
>> 
>> thank you,
>> Julien.
> 
> Julien,
> 
> I know what you mean when you say 'real df' but that's not the whole 
> story as plant has 6 more df by interacting with a single df variable. 
> There is no such thing as 'the' main effect test for plant.  The 12 df 
> test is unique and tests whether plant is associated with Y for any 
> level of leaf.age.
> 
> You can see exactly what is being tested by using various print options 
> for anova.Design, as described in the help file.  The "dots" option is 
> easy on the eyes.
> 
> Frank
> -- 
> Frank E Harrell Jr   Professor and Chair   School of Medicine
>   Department of Biostatistics   Vanderbilt University
> 
> __
> R-help@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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/anova%28cph%28..%29-output-tp23563818p23617483.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.

Re: [R] anova(cph(..) output

[Frank E Harrell Jr]

pompon wrote:

Hi,

Thank you very much for the answer.

However, I have still some misunderstandings.
from the output, can we say that plant and leaf age are significant but not
their interaction?
And the last question I promise, what would you advise me to write in the
paper to explain the different method and ackonwledge for the df?

Thank you again,
julien.


I would say there is moderate evidence for an interaction (P=0.10) and 
strong evidence for both a plant effect (at least at some level of leaf) 
and a leaf effect (at least at some level of plant).


Frank



 


Frank E Harrell Jr wrote:

pompon wrote:

Hello,

I am a beginner in R and statistics, so my question may be trivial. Sorry
in
advance.
I performed a Cox proportion hazard regression with 2 categorical
variables
with cph{design}. Then an anova on the results.
the output is 


anova(cph(surv(survival, censor) ~ plant + leaf.age + plant*leaf.age,
Mpnymph)

Wald Statistics  Response: Surv(survival,
censored) 


 FactorChi-Square
d.f. P 
 plant  (Factor+Higher Order Factors) 96.96 12   <.0001
  All Interactions   10.58 
6   0.1022

 leaf.age  (Factor+Higher Order Factors)  29.11  7   0.0001
  All Interactions 10.58 
6   0.1022

 plant * leaf.age  (Factor+Higher Order Factors)  10.58  6   0.1022
 TOTAL   106.63 13   <.0001

What do "All interaction" stand for?
The real df of for plant is 6 and 1 for leaf.age. Then, which chi square
is
one for my main factors anf their interaction.

thank you,
Julien.

Julien,

I know what you mean when you say 'real df' but that's not the whole 
story as plant has 6 more df by interacting with a single df variable. 
There is no such thing as 'the' main effect test for plant.  The 12 df 
test is unique and tests whether plant is associated with Y for any 
level of leaf.age.


You can see exactly what is being tested by using various print options 
for anova.Design, as described in the help file.  The "dots" option is 
easy on the eyes.


Frank
--
Frank E Harrell Jr   Professor and Chair   School of Medicine
  Department of Biostatistics   Vanderbilt University

__
R-help@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.







--
Frank E Harrell Jr   Professor and Chair   School of Medicine
 Department of Biostatistics   Vanderbilt University

__
R-help@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.

[R] ANOVA with replicate measures

[João Fadista]

Hi,

I have the following dataset (shown only part of it):

WFPSsampling dayslurry  depth   replicates  N2O (mg/m2/h)
1   1   1   3   1   0.007
1   1   1   3   2   0.022
1   1   1   3   3   0.015
1   1   2   3   1   0.033
1   1   2   3   2   0.010
1   1   2   3   3   0.022
1   1   2   6   1   0.038
1   1   2   6   2   0.048
1   1   2   6   3   0.043
1   1   3   0   1   0.024
1   2   1   3   1   0.017
1   2   1   3   2   0.016
1   2   1   3   3   0.015
1   2   2   3   1   0.047
1   2   2   3   2   0.083
1   2   2   3   3   0.042
1   2   2   6   1   0.036
1   2   2   6   2   0.048
1   2   2   6   3   0.030
1   2   3   0   1   0.020
1   3   1   3   1   0.006
1   3   1   3   2   0.007
1   3   1   3   3   0.019
1   3   2   3   1   0.055
1   3   2   3   2   0.004
1   3   2   3   3   0.025
1   3   2   6   1   0.030
1   3   2   6   2   0.024
1   3   2   6   3   0.042
1   3   3   0   1   0.017


What it is: daily N2O emission rate after the injection of 2 different types of 
slurry into the soil moistened at 2 different levels.
So I want to test all the factors on the N2O emissions:
1) injection depth
2) slurry type
3) water content into the soil (=WFPS)
Bu I want to take into account the replicates. How can I write the R syntax in 
this case? Something like:
 model2<-aov(N2O.mg.m2.h. ~ WFPS*sampling_day*slurry*depth + 
Error(replicates/(WFPS + sampling_day + slurry + depth)))



Thanks in advance,
Joao



[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA tables - storing F values

[Gareth Campbell]

When I run a summary(anova) I get output for all of the elements (columns)
as these are multiple - single anova results.  Can I store the F values?  I
can't find the attribute of the fitted model attributes(fit) that stores
these F values, and for that matter, P values.

Thanks

-- 
Gareth Campbell
PhD Candidate
The University of Auckland

P +649 815 3670
M +6421 256 3511
E [EMAIL PROTECTED]
[EMAIL PROTECTED]

[[alternative HTML version deleted]]

__
R-help@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.

[R] ANOVA contrast matrix vs. TukeyHSD?

[Sam Yeaman]

Dear Help List,

Thanks in advance for reading...I hope my questions are not too ignorant.

I have an experiment looking at evolution of wing size [centroid] in 
fruitflies and the effect of 6 different experimental treatments 
[treatment]. I have five replicate populations [replic] in each 
treatment and have reared the flies in two different temperatures [cond] 
to assay the wing size, making measurements on males and females 
[gender]. My design can be summarized as follows:


This is my model (I think it's right, ignoring interaction terms for 
simplicity):


> lm1 ~ aov (centroid ~ gender + cond + treatment/replic, data = parents)

The treatments are:

> levels (parents$treatment)
[1] "c"  "h"  "mc" "mh" "s"  "t"

I only care about a few of the pairwise comparisons between the levels 
of "treatment", as only certain contrasts are scientifically interesting:


c vs. h
mh vs. mc
(c + h) vs. s   [I would like to compare the mean of c and h (my 
controls) to s and t)

(c + h) vs. t
s vs. t
h vs. mh
c vs. mc

These are two more than I can specify using "contrasts()" and they are 
not orthogonal. I can use the TukeyHSD and only look at the comparisons 
I care about, but I think this should give me much less power than 
specifying a few a priori contrasts (?). Also, I don't know how to 
combine my controls (c + h) into a single comparison using TukeyHSD.



My first problem is that when I specify the matrix shown below (the 
first 5 comparisons from above), I get a much higher p-value on some of 
the planned contrasts than I do on the TukeyHSD:


contrasts (parents$treatment) <- cbind 
(c(-1,1,0,0,0,0),c(-1,-1,0,0,2,0),c(-1,-1,0,0,0,2),c(0,0,-1,1,0,0),c(0,0,0,0,1,-1)) 



> contrasts(parents$treatment)
 [,1] [,2] [,3] [,4] [,5]
c-1   -1   -100
h 1   -1   -100
mc000   -10
mh00010
s 02001
t 0020   -1


 THE OUTPUT (truncated) 

Call:
lm(formula = centroid ~ gender + cond + treatment/replic, data = parents)

Residuals:
Min1QMedian3Q   Max
-81.58846  -4.53540   0.00803   4.76568  39.84177

Coefficients: (1 not defined because of singularities)
   Estimate Std. Error  t value Pr(>|t|)   
(Intercept) 328.730960.26303 1249.770  < 2e-16 ***

genders -37.390690.19661 -190.179  < 2e-16 ***
condu   -37.477400.19693 -190.308  < 2e-16 ***
treatment10.510260.400841.273 0.203079   
treatment2   -0.173330.23175   -0.748 0.454541   
treatment30.077610.225350.344 0.730566   
treatment4   -1.960200.38524   -5.088 3.73e-07 ***
treatment5 NA NA   NA   NA  


## The TukeyHSD output (truncated) #

Tukey multiple comparisons of means
  95% family-wise confidence level

Fit: aov(formula = centroid ~ gender + cond + treatment/replic, data = 
parents)

*
*
*
$treatment
   diff  lwrupr p adj
h-c   -1.38085941 -2.382732615 -0.3789862 0.0012123
mc-c  -2.22026936 -3.198423972 -1.2421147 0.000
mh-c  -2.27157901 -3.268013478 -1.2751445 0.000
s-c   -1.19540471 -2.170272952 -0.2205365 0.0063382
t-c   -0.39899955 -1.374954044  0.5769549 0.8533107
mc-h  -0.83940995 -1.813993060  0.1351732 0.1378366
mh-h  -0.89071960 -1.883648319  0.1022091 0.1081954
s-h0.18545470 -0.785829956  1.1567394 0.9943136
t-h0.98185986  0.009484949  1.9542348 0.0462121
mh-mc -0.05130965 -1.020300865  0.9176816 0.892
s-mc   1.02486465  0.078064558  1.9716647 0.0249356
t-mc   1.82126980  0.873351301  2.7691883 0.007
s-mh   1.07617430  0.110500644  2.0418480 0.0187007
t-mh   1.87257946  0.905809220  2.8393497 0.005
t-s0.79640515 -0.148121782  1.7409321 0.1550137


When I specify the c vs. h comparison, I am getting a p-value of 
0.203079, but the TukeyHSD gives the same contrast a p-value of 
0.0012123. Also, the fifth comparison gives "NA"; I assume this is due 
to it being non-orthogonal? I feel like I am either misunderstanding the 
point of contrasts() completely or I have done something wrong, so I 
would really appreciate any help.


My other question is related...just wondering why I need to limit myself 
to only orthogonal comparisons using contrasts()? This eliminates 
comparisons of scientific interest, for example if c vs. h, mc vs. c, 
and mh vs. h are all different, I have no way of knowing if mc vs. mh is 
significantly different by examining the other contrasts.


Sorry if these questions are ignorant...I have spent a long time trying 
to figure it out and haven't found the answer in either the available 
books or the help list.


Many thanks,
Sam Yeaman

__
R-help@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, reprod

Re: [R] ANOVA between & within variance

[Gregor Rolshausen]

dear Dr. Kubovy,

I am sorry. but the Variance table is not exactly what I want. I want  
the partitioned VARIANCE for between and within the groups. the anova 
()-table just gives me the SumSq and the mean Sq... I know how to run  
t.test and ANOVA!
in the nlme-package there is the VarCorr function, which extracts the  
between and within variances, but only for nested ANOVAs. so my  
question was, if there is a function like that for not-nested ANOVAS ?

sorry. maybe I should reformulate the question.

cheers ,
gregor






Am Sep 27, 2008 um 7:19 AM schrieb Michael Kubovy:

> Than all you need is to run a t-test, no?  More generally (from ?lm):
>
> ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
> trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
> group <- gl(2,10,20, labels=c("Ctl","Trt"))
> weight <- c(ctl, trt)
> anova(lm.D9 <- lm(weight ~ group))
> This gives you what you need:
> Analysis of Variance Table
>
> Response: weight
>   Df Sum Sq Mean Sq F value Pr(>F)
> group  1   0.690.691.42   0.25
> Residuals 18   8.730.48
> I am concerned that you have not spent enough time either studying  
> stats or reading up on R. There are many good introductions to  
> stats using R.
> _
> Professor Michael Kubovy
> University of Virginia
> Department of Psychology
> USPS: P.O.Box 400400Charlottesville, VA 22904-4400
> Parcels:Room 102Gilmer Hall
> McCormick RoadCharlottesville, VA 22903
> Office:B011+1-434-982-4729
> Lab:B019+1-434-982-4751
> Fax:+1-434-982-4766
> WWW:http://www.people.virginia.edu/~mk9y/
>

Gregor Rolshausen
PhD Student
Department of Evolutionary Ecology
University of Freiburg im Breisgau; Hauptstrasse 1, 79108 Freiburg
phone - +49.761.2559
email - [EMAIL PROTECTED]




[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] ANOVA between & within variance

[Kingsford Jones]

So, rather than the estimated between and within group variances from
a standard fixed-effects ANOVA (i.e. the Mean Sqs in the anova table),
you are looking for the estimated variance components from a model
with crossed random effects?  If so, you may find the lmer function
found in package lme4 to be useful (formulas for crossed random
effects are much more straightforward than in nlme).  If not, this
page provides some helpful tips on eliciting useful responses from
this list: http://www.r-project.org/posting-guide.html

HTH,

Kingsford Jones

On Sat, Sep 27, 2008 at 2:21 AM, Gregor Rolshausen
<[EMAIL PROTECTED]> wrote:
> dear Dr. Kubovy,
>
> I am sorry. but the Variance table is not exactly what I want. I want
> the partitioned VARIANCE for between and within the groups. the anova
> ()-table just gives me the SumSq and the mean Sq... I know how to run
> t.test and ANOVA!
> in the nlme-package there is the VarCorr function, which extracts the
> between and within variances, but only for nested ANOVAs. so my
> question was, if there is a function like that for not-nested ANOVAS ?
>
> sorry. maybe I should reformulate the question.
>
> cheers ,
> gregor
>
>
>
>
>
>
> Am Sep 27, 2008 um 7:19 AM schrieb Michael Kubovy:
>
>> Than all you need is to run a t-test, no?  More generally (from ?lm):
>>
>> ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
>> trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
>> group <- gl(2,10,20, labels=c("Ctl","Trt"))
>> weight <- c(ctl, trt)
>> anova(lm.D9 <- lm(weight ~ group))
>> This gives you what you need:
>> Analysis of Variance Table
>>
>> Response: weight
>>   Df Sum Sq Mean Sq F value Pr(>F)
>> group  1   0.690.691.42   0.25
>> Residuals 18   8.730.48
>> I am concerned that you have not spent enough time either studying
>> stats or reading up on R. There are many good introductions to
>> stats using R.
>> _
>> Professor Michael Kubovy
>> University of Virginia
>> Department of Psychology
>> USPS: P.O.Box 400400Charlottesville, VA 22904-4400
>> Parcels:Room 102Gilmer Hall
>> McCormick RoadCharlottesville, VA 22903
>> Office:B011+1-434-982-4729
>> Lab:B019+1-434-982-4751
>> Fax:+1-434-982-4766
>> WWW:http://www.people.virginia.edu/~mk9y/
>>
>
> Gregor Rolshausen
> PhD Student
> Department of Evolutionary Ecology
> University of Freiburg im Breisgau; Hauptstrasse 1, 79108 Freiburg
> phone - +49.761.2559
> email - [EMAIL PROTECTED]
>
>
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@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.
>

__
R-help@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.

Re: [R] anova planned comparisons/contrasts

[Peter Alspach]

Tyler 

For balanced data like this you might find aov() gives an output which
is more comparable to Sokal and Rohlf (which I don't have):

> trtCont <- C(sugars$treatment, matrix(c(-4,1,1,1,1, 0,-1,3,-1,-1), 5,
2))
> sugarsAov <- aov(length ~ trtCont, sugars)
> summary(sugarsAov, split=list(trtCont=list('control vs rest'=1, 'gf vs
others'=2)))
   Df  Sum Sq Mean Sq  F valuePr(>F)
trtCont 4 1077.32  269.33  49.3680 6.737e-16 ***
  trtCont: control vs rest  1  832.32  832.32 152.5637 4.680e-16 ***
  trtCont: gf vs others 1   48.13   48.13   8.8228  0.004759 ** 
Residuals  45  245.505.46   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> model.tables(sugarsAov, type='mean', se=T)
Tables of means
Grand mean
  
61.94 

 trtCont 
trtCont
   control   fructose gluc+fructglucosesucrose 
  70.1   58.2   58.0   59.3   64.1 

Standard errors for differences of means
trtCont
  1.045
replic.  10

Since the two specified contrasts are orthogonal, the difference among
the remaining three sugars can be obtained by substraction:

> sugarsSum <- summary(sugarsAov,
   split=list(trtCont=list('control vs rest'=1, 'gf
vs others'=2)))
# Sum Sq
> sugarsSum[[1]][1,2]-sum(sugarsSum[[1]][2:3,2])
 
196.8667 
# Mean Sq
> (sugarsSum[[1]][1,2]-sum(sugarsSum[[1]][2:3,2]))/2
 
98.4 
# F value
> ((sugarsSum[[1]][1,2]-sum(sugarsSum[[1]][2:3,2]))/2) /
sugarsSum[[1]][4,3]
 
18.04277 
# Pr(>F)
> 1-pf(((sugarsSum[[1]][1,2]-sum(sugarsSum[[1]][2:3,2]))/2) /
sugarsSum[[1]][4,3], 2, 45)
 
1.762198e-06 

HTH 

Peter Alspach

> [mailto:[EMAIL PROTECTED] On Behalf Of Tyler Smith
> Subject: [R] anova planned comparisons/contrasts
> 
> Hi,
> 
> I'm trying to figure out how anova works in R by translating 
> the examples in Sokal And Rohlf's (1995 3rd edition) 
> Biometry. I've hit a snag with planned comparisons, their box 
> 9.4 and section 9.6. It's a basic anova design:
> 
> treatment <- factor(rep(c("control", "glucose", "fructose", 
> "gluc+fruct", "sucrose"), each = 10))
> 
> length <- c(75, 67, 70, 75, 65, 71, 67, 67, 76, 68,
> 57, 58, 60, 59, 62, 60, 60, 57, 59, 61,
> 58, 61, 56, 58, 57, 56, 61, 60, 57, 58,
> 58, 59, 58, 61, 57, 56, 58, 57, 57, 59,
> 62, 66, 65, 63, 64, 62, 65, 65, 62, 67)
> 
> sugars <- data.frame(treatment, length)
> 
> The basic analysis is easy enough:
> 
> anova(lm(length ~ treatment, sugars))
> 
> S&R proceed to calculate planned comparisons between control 
> and all other groups, between gluc+fruct and the remaining 
> sugars, and among the three pure sugars. I can get the first 
> two comparisons using the
> contrasts() function:
> 
> contrasts(sugars$treatment) <- matrix(c(-4, 1, 1,  1,  1,
> 0, -1, 3, -1, -1), 5, 2) 
> 
> summary(lm(length ~ treatment, sugars))
> 
> The results appear to be the same, excepting that the book 
> calculates an F value, whereas R produces an equivalent t 
> value. However, I can't figure out how to calculate a 
> contrast among three groups. For those without access to 
> Sokal and Rohlf, I've written an R function that performs the 
> analysis they use, copied below. Is there a better way to do 
> this in R?
> 
> Also, I don't know how to interpret the Estimate and Std. 
> Error columns from the summary of the lm with contrasts. It's 
> clear the intercept in this case is the grand mean, but what 
> are the other values? They appear to be the difference 
> between one of the contrast groups and the grand mean -- 
> wouldn't an estimate of the differences between the 
> contrasted groups be more appropriate, or am I (likely) 
> misinterpreting this?
> 
> Thanks!
> 
> Tyler
[Code deleted]

__
R-help@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.

Re: [R] anova planned comparisons/contrasts

[Tyler Smith]

On 2007-11-22, Peter Alspach <[EMAIL PROTECTED]> wrote:
> Tyler 
>
> For balanced data like this you might find aov() gives an output which
> is more comparable to Sokal and Rohlf (which I don't have):
>
>> trtCont <- C(sugars$treatment, matrix(c(-4,1,1,1,1, 0,-1,3,-1,-1), 5,
> 2))
>> sugarsAov <- aov(length ~ trtCont, sugars)
>> summary(sugarsAov, split=list(trtCont=list('control vs rest'=1, 'gf vs
> others'=2)))

>> model.tables(sugarsAov, type='mean', se=T)

Thank you Peter, that's a big help! To confirm that I understand
correctly, aov is identical to lm, but provides better summary
information for balanced anova designs. As such, it is preferred to lm
for balanced anova designs, but should be avoided otherwise. Is that
correct?

Also, it appears that C and contrasts serve pretty much the same
purpose. Is there a context in which one is preferable to the other? 

Cheers,

Tyler

__
R-help@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.

Re: [R] anova planned comparisons/contrasts

[Dieter Menne]

Tyler Smith  mail.mcgill.ca> writes:

> I'm trying to figure out how anova works in R by translating the
> examples in Sokal And Rohlf's (1995 3rd edition) Biometry. I've hit a
> snag with planned comparisons, their box 9.4 and section 9.6. It's a
> basic anova design:
> 

 how to do contrast

It's easier to use function estimable in package gmodels, or glht in package
multcomp.

Dieter

__
R-help@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.

Re: [R] anova planned comparisons/contrasts

[Tyler Smith]

On 2007-11-22, Dieter Menne <[EMAIL PROTECTED]> wrote:
> Tyler Smith  mail.mcgill.ca> writes:
>
>> I'm trying to figure out how anova works in R by translating the
>> examples in Sokal And Rohlf's (1995 3rd edition) Biometry. I've hit a
>> snag with planned comparisons, their box 9.4 and section 9.6. It's a
>> basic anova design:
>> 
>
>  how to do contrast
>
> It's easier to use function estimable in package gmodels, or glht in package
> multcomp.
>

Isn't glht calculating unplanned comparisons, as opposed to the
planned comparisons with contrasts() and summary(..., split= ...)? Can
you do planned comparisons with glht, or unplanned comparisons with
summary()? 

contrasts(warpbreaks$tension) <- matrix(c(-1,1,0,1,0,-1,0,-1,1), 3, 3)
amod <- aov(formula = breaks ~ tension, data = warpbreaks)

## this isn't right:
summary(amod, split = list(tension = list('L vs M' =1, 'L vs H'=2, 
  'M vs H' = 3)))

## posthoc contrasts (three ways to do the same test, I think):
summary(glht(amod, linfct = mcp(tension = 'Tukey')))
summary(glht(amod, linfct = mcp(tension = 
matrix(c(-1,1,0,1,0,-1,0,-1,1), 3, 3
TukeyHSD(amod)

Thanks,

Tyler

__
R-help@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.

Re: [R] anova planned comparisons/contrasts

[Greg Snow]

Others have shown some approaches that work well for after you fit the
model.  Here is another approach starting with the model fit itself:

tmp <- c("control", "glucose", "fructose", 
  "gluc+fruct", "sucrose")

treatment <- factor(rep( tmp, each=10 ), levels=tmp)

length <- c(75, 67, 70, 75, 65, 71, 67, 67, 76, 68,
57, 58, 60, 59, 62, 60, 60, 57, 59, 61,
58, 61, 56, 58, 57, 56, 61, 60, 57, 58,
58, 59, 58, 61, 57, 56, 58, 57, 57, 59,
62, 66, 65, 63, 64, 62, 65, 65, 62, 67)



cm1 <- rbind( int=1/5, cntrl = c(4, -1,-1,-1,-1)/4,
plus=c( 0, -1, -1, 3, -1 )/3,
sug1=c( 0, -1, 1, 0, 0),
sug2=c( 0, -1, -1, 0, 2)/2 )

cm2 <- zapsmall(solve(cm1))

contrasts(treatment) <- cm2[,-1]

treatment

sugars <- data.frame( length=length, treatment=treatment )

fit <- aov( length ~ treatment, data=sugars )
summary.lm(fit)
summary(fit)
summary(fit, split=list(treatment=list(
  'Control vs. Sugars'=1, 'Gluc+Fruc vs. sugars'=2, 'Sugars'=3:4)))


Is this more what you were looking for?

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111
 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Tyler Smith
> Sent: Thursday, November 22, 2007 12:38 PM
> To: [EMAIL PROTECTED]
> Subject: [R] anova planned comparisons/contrasts
> 
> Hi,
> 
> I'm trying to figure out how anova works in R by translating 
> the examples in Sokal And Rohlf's (1995 3rd edition) 
> Biometry. I've hit a snag with planned comparisons, their box 
> 9.4 and section 9.6. It's a basic anova design:
> 
> treatment <- factor(rep(c("control", "glucose", "fructose", 
> "gluc+fruct", "sucrose"), each = 10))
> 
> length <- c(75, 67, 70, 75, 65, 71, 67, 67, 76, 68,
> 57, 58, 60, 59, 62, 60, 60, 57, 59, 61,
> 58, 61, 56, 58, 57, 56, 61, 60, 57, 58,
> 58, 59, 58, 61, 57, 56, 58, 57, 57, 59,
> 62, 66, 65, 63, 64, 62, 65, 65, 62, 67)
> 
> sugars <- data.frame(treatment, length)
> 
> The basic analysis is easy enough:
> 
> anova(lm(length ~ treatment, sugars))
> 
> S&R proceed to calculate planned comparisons between control 
> and all other groups, between gluc+fruct and the remaining 
> sugars, and among the three pure sugars. I can get the first 
> two comparisons using the
> contrasts() function:
> 
> contrasts(sugars$treatment) <- matrix(c(-4, 1, 1,  1,  1,
> 0, -1, 3, -1, -1), 5, 2) 
> 
> summary(lm(length ~ treatment, sugars))
> 
> The results appear to be the same, excepting that the book 
> calculates an F value, whereas R produces an equivalent t 
> value. However, I can't figure out how to calculate a 
> contrast among three groups. For those without access to 
> Sokal and Rohlf, I've written an R function that performs the 
> analysis they use, copied below. Is there a better way to do 
> this in R?
> 
> Also, I don't know how to interpret the Estimate and Std. 
> Error columns from the summary of the lm with contrasts. It's 
> clear the intercept in this case is the grand mean, but what 
> are the other values? They appear to be the difference 
> between one of the contrast groups and the grand mean -- 
> wouldn't an estimate of the differences between the 
> contrasted groups be more appropriate, or am I (likely) 
> misinterpreting this?
> 
> Thanks!
> 
> Tyler
> 
> MyContrast <- function(Var, Group, G1, G2, G3=NULL) {
>   ## Var == data vector, Group == factor
>   ## G1, G2, G3 == character vectors of factor levels to contrast
>   
>   nG1 = sum(Group %in% G1)
>   nG2 = sum(Group %in% G2)
>   GrandMean = mean(Var[Group %in% c(G1, G2, G3)])
>   G1Mean = mean(Var[Group %in% G1])
>   G2Mean = mean(Var[Group %in% G2])
> 
>   if(is.null(G3))
> MScontr = nG1 * ((G1Mean - GrandMean)^2) + 
> nG2 * ((G2Mean - GrandMean)^2)
>   else {
>   nG3 = sum(Group %in% G3)
>   G3Mean = mean(Var[Group %in% G3])
>   MScontr = (nG1 * ((G1Mean - GrandMean)^2) + 
>nG2 * ((G2Mean - GrandMean)^2) + 
>nG3 * ((G3Mean - GrandMean)^2))/2
> }
>   
>   An <- anova(lm(Var ~ Group))
>   MSwithin = An[3]['Residuals',]
>   DegF = An$Df[length(An$Df)]
>   Fval = MScontr / MSwithin
>   Pval = 1 - pf(Fval, 1, DegF)
>   
>   return (list(MS_contrasts = MScontr, MS_within = MSwithin, 
>  F_value = Fval, P_va

Re: [R] Anova(car) SS digits

[Tyler Smith]

On 2007-11-29, John Fox <[EMAIL PROTECTED]> wrote:
>
> Anova() creates an object of class "anova" which then gets printed by the
> print method for "anova" objects. The reason that the sums of squares for
> the type II tests are rounded like this is because of the large value for
> the intercept. Although Anova() doesn't take a digits argument,
> print.anova() does, and so, e.g., you could use the command
>
> print(Anova(bot.lm3, type ="III"), digits=7)
>
> I hope this helps,
>  John

Very much. Thanks!

Tyler

__
R-help@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.

Re: [R] Anova(car) SS digits

[John Fox]

Dear Tyler,

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Tyler Smith
> Sent: Thursday, November 29, 2007 2:07 PM
> To: [EMAIL PROTECTED]
> Subject: [R] Anova(car) SS digits
> 
> Hi,
> 
> When I use Anova(car) to produce type III SS, 'Sum Sq' is reported in
> integers:
> 
> > Anova(bot.lm3, type ="III")
> Anova Table (Type III tests)
> 
> Response: bottemp
>  Sum Sq  DfF value  Pr(>F)
> (Intercept)   45295   1 29436.4440 < 2e-16
> fungroup  3   2 0.8259 0.44006
> numsp.fun11   2 3.4503 0.03460
> block   419   468.0017 < 2e-16
> fungroup:numsp.fun6   4 1.0317 0.39340
> fungroup:block9   8 0.7429 0.65346
> numsp.fun:block  12   8 0.9501 0.47795
> fungroup:numsp.fun:block 27  16 1.0920 0.36882
> Residuals   205 133   
> 
> I have tried specifying digits = 7 in the call to Anova, but 
> that doesn't change anything. The output from another stats 
> program for the same data confirms that the ouptut above is 
> correct, but is rounded to the nearest integer. The dataset 
> is too large to post here, but I can confirm that I get the 
> expected results from the examples in ?Anova, as well as for 
> type = 'II' on my own data. This output, and relevant
> options() and sessionInfo are pasted below. What am I doing wrong?

Anova() creates an object of class "anova" which then gets printed by the
print method for "anova" objects. The reason that the sums of squares for
the type II tests are rounded like this is because of the large value for
the intercept. Although Anova() doesn't take a digits argument,
print.anova() does, and so, e.g., you could use the command

print(Anova(bot.lm3, type ="III"), digits=7)

I hope this helps,
 John

> 
> Thanks,
> 
> Tyler
> 
> >  mod <- lm(conformity ~ fcategory*partner.status, data=Moore,
> +contrasts=list(fcategory=contr.sum, 
> partner.status=contr.sum))
> >  Anova(mod)
> Anova Table (Type II tests)
> 
> Response: conformity
>  Sum Sq Df F value   Pr(>F)
> fcategory 11.61  2  0.2770 0.759564
> partner.status   212.21  1 10.1207 0.002874
> fcategory:partner.status 175.49  2  4.1846 0.022572
> Residuals817.76 39 
> >  Anova(mod, type="III")
> Anova Table (Type III tests)
> 
> Response: conformity
>  Sum Sq Df  F valuePr(>F)
> (Intercept)  5752.8  1 274.3592 < 2.2e-16
> fcategory  36.0  2   0.8589  0.431492
> partner.status239.6  1  11.4250  0.001657
> fcategory:partner.status  175.5  2   4.1846  0.022572
> Residuals 817.8 39   
> 
> > Anova(bot.lm3, type ="II")
> Anova Table (Type II tests)
> 
> Response: bottemp
>  Sum Sq  Df F value  Pr(>F)
> fungroup   1.87   2  0.6069 0.54656
> numsp.fun 11.14   2  3.6214 0.02942
> block486.47   4 79.0379 < 2e-16
> fungroup:numsp.fun 6.32   4  1.0275 0.39553
> fungroup:block11.43   8  0.9283 0.49542
> numsp.fun:block   11.29   8  0.9169 0.50473
> fungroup:numsp.fun:block  26.89  16  1.0920 0.36882
> Residuals204.65 133
> 
> > options('digits')
> $digits
> [1] 7
> 
> > options('contrasts')
> $contrasts
> [1] "contr.sum"  "contr.poly"
> 
> > sessionInfo()
> R version 2.5.1 (2007-06-27)
> i486-pc-linux-gnu 
> 
> locale:
> 
> LC_CTYPE=en_CA.UTF-8;LC_NUMERIC=C;LC_TIME=en_CA.UTF-8;LC_COLLA
> TE=en_CA.UTF-8;
> LC_MONETARY=en_CA.UTF-8;LC_MESSAGES=en_CA.UTF-8;LC_PAPER=en_CA
> .UTF-8;LC_NAME=C;
> LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_CA.UTF-8;LC_IDEN
> TIFICATION=C
> 
> attached base packages:
> [1] "stats" "graphics"  "grDevices" "utils" 
> "datasets"  "methods"  
> [7] "base" 
> 
> other attached packages:
> car
> "1.2-7"
> 
> __
> R-help@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.
>

__
R-help@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.

[R] Anova for stratified Cox regression

[Matthias Gondan]

Dear List,

I have tried a stratified Cox Regression, it is working fine, except for
the "Anova"-Tests:

Here the commands (should work out of the box):

library(survival)
d = colon[colon$etype==2, ]
m = coxph(Surv(time, status) ~ strata(sex) + rx, data=d)
summary(m)
# Printout ok
anova(m, test='Chisq')

This is the output of the anova command:

> Analysis of Deviance Table
>  Cox model: response is Surv(time, status)
> Terms added sequentially (first to last)
>
>  Df Deviance Resid. Df Resid. Dev P(>|Chi|)
> NULL   929 5233.5  
> strata(sex)   0929 
> rx2927 5221.2  
> Warning message:
> In is.na(coef(fit)) :
>   is.na() auf nicht-(Liste oder Vektor) des Typs 'NULL' angewendet

It should be possible to do Chi-Square-Tests in a stratified analysis, 
right?

Best wishes,

Matthias

__
R-help@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.

Re: [R] Anova and unbalanced designs

[John Fox]

Dear Nils,

This is a pretty simple design, and I wouldn't have thought that there was
much room for getting different results. More generally, but not here (since
there's only one between-subject factor), one shouldn't use
contr.treatment() with "type-III" tests, as you did. Is it possible that you
got "type-II" tests from SPSS:

-- snip --

> summary(Anova(betweenanova, idata=with, idesign= ~within, type = "II" ))

Type II Repeated Measures MANOVA Tests:

--
 
Term: between 

 Response transformation matrix:
   (Intercept)
w1   1
w2   1

Sum of squares and products for the hypothesis:
(Intercept)
(Intercept) 9.6

Sum of squares and products for error:
(Intercept)
(Intercept)  18

Multivariate Tests: between
 Df test stat approx F num Df den Df   Pr(>F)  
Pillai1  0.347826 4.27  1  8 0.072726 .
Wilks 1  0.652174 4.27  1  8 0.072726 .
Hotelling-Lawley  1  0.53 4.27  1  8 0.072726 .
Roy   1  0.53 4.27  1  8 0.072726 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

--
 
Term: within 

 Response transformation matrix:
   within1
w1   1
w2  -1

Sum of squares and products for the hypothesis:
within1
within1 0.4

Sum of squares and products for error:
 within1
within1 21.3

Multivariate Tests: within
 Df test stat  approx F num Df den Df  Pr(>F)
Pillai1 0.0184049 0.150  1  8 0.70864
Wilks 1 0.9815951 0.150  1  8 0.70864
Hotelling-Lawley  1 0.0187500 0.150  1  8 0.70864
Roy   1 0.0187500 0.150  1  8 0.70864

--
 
Term: between:within 

 Response transformation matrix:
   within1
w1   1
w2  -1

Sum of squares and products for the hypothesis:
 within1
within1 4.27

Sum of squares and products for error:
 within1
within1 21.3

Multivariate Tests: between:within
 Df test stat  approx F num Df den Df  Pr(>F)
Pillai1 0.167 1.600  1  8 0.24150
Wilks 1 0.833 1.600  1  8 0.24150
Hotelling-Lawley  1 0.200 1.600  1  8 0.24150
Roy   1 0.200 1.600  1  8 0.24150

Univariate Type II Repeated-Measures ANOVA Assuming Sphericity

SS num Df Error SS den Df  F  Pr(>F)  
between 4.8000  1   9.  8 4.2667 0.07273 .
within  0.2000  1  10.6667  8 0.1500 0.70864  
between:within  2.1333  1  10.6667  8 1.6000 0.24150  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

-- snip --

I hope this helps,
 John

--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
> Behalf Of Skotara
> Sent: January-23-09 12:16 PM
> To: r-help@r-project.org
> Subject: [R] Anova and unbalanced designs
> 
> Dear R-list!
> 
> My question is related to an Anova including within and between subject
> factors and unequal group sizes.
> Here is a minimal example of what I did:
> 
> library(car)
> within1 <- c(1,2,3,4,5,6,4,5,3,2); within2 <- c(3,4,3,4,3,4,3,4,5,4)
> values <- data.frame(w1 = within1, w2 = within2)
> values <- as.matrix(values)
> between <- factor(c(rep(1,4), rep(2,6)))
> betweenanova <- lm(values ~ between)
> with <- expand.grid(within = factor(1:2))
> withinanova <- Anova(betweenanova, idata=with, idesign=
> ~as.factor(within), type = "III" )
> 
> I do not know if this is the appropriate method to deal with unbalanced
> designs.
> 
> I observed, that SPSS calculates everything identically except the main
> effect of the within factor, here, the SSQ and F-value are very different
> If selecting the option "show means", the means for the levels of the
> within factor in SPSS are the same as:
> mean(c(mean(values$w1[1:4]),mean(values$w1[5:10]))) and
> mean(c(mean(values$w2[1:4]),mean(values$w2[5:10]))).
> In other words, they are calculated as if both groups would have the
> same size.
> 
> I wonder if this is a good solution and if so, how could I do the same
> thing in R?
> However, I think if this is treated in SPSS as if the group sizes are
> identical,
> then why not the interaction, which yields to the same result as using
> Anova()?
> 

Re: [R] Anova and unbalanced designs

[Skotara]

Dear John,

thank you for your answer. You are right, I also would not have expected 
a divergent result.

I have double-checked it again. No, I got type-III tests.
When I use type II, I get the same results in SPSS as in 'Anova' (using 
also type-II tests).
My guess was that the somehow weighted means SPSS shows could be 
responsible for this difference.
Or that using 'Anova' would not be correct for unequal group n's, which 
was not the case I think.

Do you have any further ideas?

Thank you!
Nils

John Fox schrieb:

Dear Nils,

This is a pretty simple design, and I wouldn't have thought that there was
much room for getting different results. More generally, but not here (since
there's only one between-subject factor), one shouldn't use
contr.treatment() with "type-III" tests, as you did. Is it possible that you
got "type-II" tests from SPSS:

-- snip --

  

summary(Anova(betweenanova, idata=with, idesign= ~within, type = "II" ))



Type II Repeated Measures MANOVA Tests:

--
 
Term: between 


 Response transformation matrix:
   (Intercept)
w1   1
w2   1

Sum of squares and products for the hypothesis:
(Intercept)
(Intercept) 9.6

Sum of squares and products for error:
(Intercept)
(Intercept)  18

Multivariate Tests: between
 Df test stat approx F num Df den Df   Pr(>F)  
Pillai1  0.347826 4.27  1  8 0.072726 .

Wilks 1  0.652174 4.27  1  8 0.072726 .
Hotelling-Lawley  1  0.53 4.27  1  8 0.072726 .
Roy   1  0.53 4.27  1  8 0.072726 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 


--
 
Term: within 


 Response transformation matrix:
   within1
w1   1
w2  -1

Sum of squares and products for the hypothesis:
within1
within1 0.4

Sum of squares and products for error:
 within1
within1 21.3

Multivariate Tests: within
 Df test stat  approx F num Df den Df  Pr(>F)
Pillai1 0.0184049 0.150  1  8 0.70864
Wilks 1 0.9815951 0.150  1  8 0.70864
Hotelling-Lawley  1 0.0187500 0.150  1  8 0.70864
Roy   1 0.0187500 0.150  1  8 0.70864

--
 
Term: between:within 


 Response transformation matrix:
   within1
w1   1
w2  -1

Sum of squares and products for the hypothesis:
 within1
within1 4.27

Sum of squares and products for error:
 within1
within1 21.3

Multivariate Tests: between:within
 Df test stat  approx F num Df den Df  Pr(>F)
Pillai1 0.167 1.600  1  8 0.24150
Wilks 1 0.833 1.600  1  8 0.24150
Hotelling-Lawley  1 0.200 1.600  1  8 0.24150
Roy   1 0.200 1.600  1  8 0.24150

Univariate Type II Repeated-Measures ANOVA Assuming Sphericity

SS num Df Error SS den Df  F  Pr(>F)  
between 4.8000  1   9.  8 4.2667 0.07273 .
within  0.2000  1  10.6667  8 0.1500 0.70864  
between:within  2.1333  1  10.6667  8 1.6000 0.24150  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 


-- snip --

I hope this helps,
 John

--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


  

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]


On
  

Behalf Of Skotara
Sent: January-23-09 12:16 PM
To: r-help@r-project.org
Subject: [R] Anova and unbalanced designs

Dear R-list!

My question is related to an Anova including within and between subject
factors and unequal group sizes.
Here is a minimal example of what I did:

library(car)
within1 <- c(1,2,3,4,5,6,4,5,3,2); within2 <- c(3,4,3,4,3,4,3,4,5,4)
values <- data.frame(w1 = within1, w2 = within2)
values <- as.matrix(values)
between <- factor(c(rep(1,4), rep(2,6)))
betweenanova <- lm(values ~ between)
with <- expand.grid(within = factor(1:2))
withinanova <- Anova(betweenanova, idata=with, idesign=
~as.factor(within), type = "III" )

I do not know if this is the appropriate method to deal with unbalanced
designs.

I observed, that SPSS calculates everything identically except the main
effect of the within factor, here, the SSQ and F-value are very different
If selecting the option "show means", the means for the levels of the
within factor in SPSS are the same as:
mean(c(mean(values$w1[1:4]),mean(values$w1[5:10]))) and
mean(c(mea

Re: [R] Anova and unbalanced designs

[John Fox]

imple, that I find it hard to understand
where there's room for error, but I wanted to check against SAS to test my
sanity (a procedure that will likely get a rise out of some list members).

Maybe you should send a message to the SPSS help list.

Regards,
 John

--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
> Behalf Of Skotara
> Sent: January-24-09 6:30 AM
> To: John Fox
> Cc: r-help@r-project.org
> Subject: Re: [R] Anova and unbalanced designs
> 
> Dear John,
> 
> thank you for your answer. You are right, I also would not have expected
> a divergent result.
> I have double-checked it again. No, I got type-III tests.
> When I use type II, I get the same results in SPSS as in 'Anova' (using
> also type-II tests).
> My guess was that the somehow weighted means SPSS shows could be
> responsible for this difference.
> Or that using 'Anova' would not be correct for unequal group n's, which
> was not the case I think.
> Do you have any further ideas?
> 
> Thank you!
> Nils
> 
> John Fox schrieb:
> > Dear Nils,
> >
> > This is a pretty simple design, and I wouldn't have thought that there
was
> > much room for getting different results. More generally, but not here
> (since
> > there's only one between-subject factor), one shouldn't use
> > contr.treatment() with "type-III" tests, as you did. Is it possible that
> you
> > got "type-II" tests from SPSS:
> >
> > -- snip --
> >
> >
> >> summary(Anova(betweenanova, idata=with, idesign= ~within, type = "II"
))
> >>
> >
> > Type II Repeated Measures MANOVA Tests:
> >
> > --
> >
> > Term: between
> >
> >  Response transformation matrix:
> >(Intercept)
> > w1   1
> > w2   1
> >
> > Sum of squares and products for the hypothesis:
> > (Intercept)
> > (Intercept) 9.6
> >
> > Sum of squares and products for error:
> > (Intercept)
> > (Intercept)  18
> >
> > Multivariate Tests: between
> >  Df test stat approx F num Df den Df   Pr(>F)
> > Pillai1  0.347826 4.27  1  8 0.072726 .
> > Wilks 1  0.652174 4.27  1  8 0.072726 .
> > Hotelling-Lawley  1  0.53 4.27  1  8 0.072726 .
> > Roy   1  0.53 4.27  1  8 0.072726 .
> > ---
> > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >
> > --
> >
> > Term: within
> >
> >  Response transformation matrix:
> >within1
> > w1   1
> > w2  -1
> >
> > Sum of squares and products for the hypothesis:
> > within1
> > within1 0.4
> >
> > Sum of squares and products for error:
> >  within1
> > within1 21.3
> >
> > Multivariate Tests: within
> >  Df test stat  approx F num Df den Df  Pr(>F)
> > Pillai1 0.0184049 0.150  1  8 0.70864
> > Wilks 1 0.9815951 0.150  1  8 0.70864
> > Hotelling-Lawley  1 0.0187500 0.150  1  8 0.70864
> > Roy   1 0.0187500 0.150  1  8 0.70864
> >
> > --
> >
> > Term: between:within
> >
> >  Response transformation matrix:
> >within1
> > w1   1
> > w2  -1
> >
> > Sum of squares and products for the hypothesis:
> >  within1
> > within1 4.27
> >
> > Sum of squares and products for error:
> >  within1
> > within1 21.3
> >
> > Multivariate Tests: between:within
> >  Df test stat  approx F num Df den Df  Pr(>F)
> > Pillai1 0.167 1.600  1  8 0.24150
> > Wilks 1 0.833 1.600  1  8 0.24150
> > Hotelling-Lawley  1 0.200 1.600  1  8 0.24150
> > Roy   1 0.200 1.600  1      8 0.24150
> >
> > Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
> >
> > SS num Df Error SS den Df  F  Pr(>F)
> > between 4.8000  1   9.  8 4

Re: [R] Anova and unbalanced designs

[Nils Skotara]

riance
>Univariate Tests of Hypotheses for Within Subject Effects
> 
>  Source DFType III SSMean Square   F Value   Pr
> > F
> 
>  within  1 0.5333 0.5333  0.40
> 0.5447
>  within*between  1 2.1333 2.1333  1.60
> 0.2415
>  Error(within)   810.6667 1.
> 
> --- snip 
> 
> As you can see, these agree with Anova():
> 
> --- snip 
> 
> Type III Repeated Measures MANOVA Tests: Pillai test statistic
>Df test stat approx F num Df den DfPr(>F)
> (Intercept) 1 0.963  209.067  1  8 5.121e-07 ***
> between 1 0.3484.267  1  8   0.07273 .
> within  1 0.0480.400  1  8   0.54474
> between:within  1 0.1671.600  1  8   0.24150
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> 
> Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
> 
> SS num Df Error SS den DfFPr(>F)
> (Intercept)235.200  19.000  8 209.0667 5.121e-07 ***
> between  4.800  19.000  8   4.2667   0.07273 .
> within   0.533  1   10.667  8   0.4000   0.54474
> between:within   2.133  1   10.667  8   1.6000   0.24150
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> --- snip 
> 
> So, unless Anova() and SAS are making the same error, I guess SPSS is doing
> something strange (or perhaps you didn't do what you intended in SPSS). As I
> said before, this problem is so simple, that I find it hard to understand
> where there's room for error, but I wanted to check against SAS to test my
> sanity (a procedure that will likely get a rise out of some list members).
> 
> Maybe you should send a message to the SPSS help list.
> 
> Regards,
>  John
> 
> --
> John Fox, Professor
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> web: socserv.mcmaster.ca/jfox
> 
> 
> > -Original Message-
> > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On
> > Behalf Of Skotara
> > Sent: January-24-09 6:30 AM
> > To: John Fox
> > Cc: r-help@r-project.org
> > Subject: Re: [R] Anova and unbalanced designs
> >
> > Dear John,
> >
> > thank you for your answer. You are right, I also would not have expected
> > a divergent result.
> > I have double-checked it again. No, I got type-III tests.
> > When I use type II, I get the same results in SPSS as in 'Anova' (using
> > also type-II tests).
> > My guess was that the somehow weighted means SPSS shows could be
> > responsible for this difference.
> > Or that using 'Anova' would not be correct for unequal group n's, which
> > was not the case I think.
> > Do you have any further ideas?
> >
> > Thank you!
> > Nils
> >
> > John Fox schrieb:
> > > Dear Nils,
> > >
> > > This is a pretty simple design, and I wouldn't have thought that there
> was
> > > much room for getting different results. More generally, but not here
> > (since
> > > there's only one between-subject factor), one shouldn't use
> > > contr.treatment() with "type-III" tests, as you did. Is it possible that
> > you
> > > got "type-II" tests from SPSS:
> > >
> > > -- snip --
> > >
> > >
> > >> summary(Anova(betweenanova, idata=with, idesign= ~within, type = "II"
> ))
> > >>
> > >
> > > Type II Repeated Measures MANOVA Tests:
> > >
> > > --
> > >
> > > Term: between
> > >
> > >  Response transformation matrix:
> > >(Intercept)
> > > w1   1
> > > w2   1
> > >
> > > Sum of squares and products for the hypothesis:
> > > (Intercept)
> > > (Intercept) 9.6
> > >
> > > Sum of squares and products for error:
> > > (Intercept)
> > > (Intercept)  18
> > >
> > > Multivariate Tests: between
> > >  Df test stat approx F num Df den Df   Pr(>F)
> > > Pillai1  0.347826 4.27  

Re: [R] Anova and unbalanced designs

[Peter Dalgaard]

Nils Skotara wrote:
Dear John, 


thank you again! You replicated the type III result I got in SPSS! When I
calculate Anova() type II:

Univariate Type II Repeated-Measures ANOVA Assuming Sphericity

SS num Df Error SS den Df  F  Pr(>F)  
between 4.8000  1   9.  8 4.2667 0.07273 .
within  0.2000  1  10.6667  8 0.1500 0.70864  
between:within  2.1333  1  10.6667  8 1.6000 0.24150  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

I see the exact same values as you had written. 
However, and now I am really lost, type III (I did not change anything else)
leads to the following: 


Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

  SS num Df Error SS den Df   FPr(>F)
(Intercept)   72.000  19.000  8 64. 4.367e-05 ***
between4.800  19.000  8  4.2667   0.07273 .  
as.factor(within)  2.000  1   10.667  8  1.5000   0.25551
between:as.factor(within)  2.133  1   10.667  8  1.6000   0.24150
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

How is this possible? 


This looks like a contrast parametrization issue: If we look at the 
per-group mean within-differences and their SE, we get


> summary(lm(within1-within2~between - 1))
..
Coefficients:
 Estimate Std. Error t value Pr(>|t|)
between1  -1. 0.8165  -1.2250.256
between2   0. 0.6667   0.5000.631
..
> table(between)
between
1 2
4 6

Now, the type II F test is based on weighting the two means as you would 
after testing for no interaction


> (4*-1+6*.)^2/(4^2*0.8165^2+6^2*0.6667^2)
[1] 0.1500205

and type III is to weight them as if there had been equal counts

> (5*-1+5*.)^2/(5^2*0.8165^2+5^2*0.6667^2)
[1] 0.400022

However, the result above corresponds to looking at group1 only

> (-1)^2/(0.8165^2)
[1] 1.499987

It helps if you choose orhtogonal contrast parametrizations:

> options(contrasts=c("contr.sum","contr.helmert"))
> betweenanova <- lm(values ~ between)> Anova(betweenanova, idata=with, 
idesign= ~as.factor(within), type = "III" )


Type III Repeated Measures MANOVA Tests: Pillai test statistic
  Df test stat approx F num Df den DfPr(>F)
(Intercept)1 0.963  209.067  1  8 5.121e-07 ***
between1 0.3484.267  1  8   0.07273 .
as.factor(within)  1 0.0480.400  1  8   0.54474
between:as.factor(within)  1 0.1671.600  1  8   0.24150
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1




--
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

__
R-help@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.

Re: [R] Anova and unbalanced designs

[John Fox]

Dear Peter and Nils,

In my initial message, I stated misleadingly that the contrast coding didn't
matter for the "type-III" tests here since there is just one
between-subjects factor, but that's not right: The between type-III SS is
correct using contr.treatment(), but the within SS is not. As is generally
the case, to get reasonable type-III tests (i.e., tests of reasonable
hypotheses), it's necessary to have contrasts that are orthogonal in the
row-basis of the design, such as contr.sum(),  contr.helmert(), or
contr.poly(). The "type-II" tests, however, are insensitive to the contrast
parametrization. Anova() always uses an orthogonal parametrization for the
within-subjects design.

The general advice in ?Anova is, "Be very careful in formulating the model
for type-III tests, or the hypotheses tested will not make sense."

Thanks, Peter, for pointing this out.

John

--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


> -Original Message-
> From: Peter Dalgaard [mailto:p.dalga...@biostat.ku.dk]
> Sent: January-24-09 6:31 PM
> To: Nils Skotara
> Cc: John Fox; r-help@r-project.org; 'Michael Friendly'
> Subject: Re: [R] Anova and unbalanced designs
> 
> Nils Skotara wrote:
> > Dear John,
> >
> > thank you again! You replicated the type III result I got in SPSS! When
I
> > calculate Anova() type II:
> >
> > Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
> >
> > SS num Df Error SS den Df  F  Pr(>F)
> > between 4.8000  1   9.  8 4.2667 0.07273 .
> > within  0.2000  1  10.6667  8 0.1500 0.70864
> > between:within  2.1333  1  10.6667  8 1.6000 0.24150
> > ---
> > Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> > I see the exact same values as you had written.
> > However, and now I am really lost, type III (I did not change anything
> else)
> > leads to the following:
> >
> > Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
> >
> >   SS num Df Error SS den Df   F
Pr(>F)
> > (Intercept)   72.000  19.000  8 64.
4.367e-05
> ***
> > between4.800  19.000  8  4.2667
0.07273 .
> > as.factor(within)  2.000  1   10.667  8  1.5000
0.25551
> > between:as.factor(within)  2.133  1   10.667  8  1.6000
0.24150
> > ---
> > Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> > How is this possible?
> 
> This looks like a contrast parametrization issue: If we look at the
> per-group mean within-differences and their SE, we get
> 
>  > summary(lm(within1-within2~between - 1))
> ..
> Coefficients:
>   Estimate Std. Error t value Pr(>|t|)
> between1  -1. 0.8165  -1.2250.256
> between2   0. 0.6667   0.5000.631
> ..
>  > table(between)
> between
> 1 2
> 4 6
> 
> Now, the type II F test is based on weighting the two means as you would
> after testing for no interaction
> 
>  > (4*-1+6*.)^2/(4^2*0.8165^2+6^2*0.6667^2)
> [1] 0.1500205
> 
> and type III is to weight them as if there had been equal counts
> 
>  > (5*-1+5*.)^2/(5^2*0.8165^2+5^2*0.6667^2)
> [1] 0.400022
> 
> However, the result above corresponds to looking at group1 only
> 
>  > (-1)^2/(0.8165^2)
> [1] 1.499987
> 
> It helps if you choose orhtogonal contrast parametrizations:
> 
>  > options(contrasts=c("contr.sum","contr.helmert"))
>  > betweenanova <- lm(values ~ between)> Anova(betweenanova, idata=with,
> idesign= ~as.factor(within), type = "III" )
> 
> Type III Repeated Measures MANOVA Tests: Pillai test statistic
>Df test stat approx F num Df den DfPr(>F)
> (Intercept)1 0.963  209.067  1  8 5.121e-07
***
> between1 0.3484.267  1  8   0.07273 .
> as.factor(within)  1 0.0480.400  1  8   0.54474
> between:as.factor(within)  1 0.1671.600  1  8   0.24150
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> 
> 
> 
> --
> O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
>c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
>   (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
> ~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

__
R-help@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.

Re: [R] anova p value extraction

[H. Paul Benton]

Yea the anova object seems to be odd. It's not S4 so that's why I tried 
originally the attr() funtion but

summary(pb)$Pr(>F)

Error: unexpected '>' in "summary(pb)$Pr(>"
> summary(pb)$Pr

NULL
> summary(pb)@Pr(>F)

Error: unexpected '>' in "summary(pb)@Pr(>"
> summary(pb)@Pr(F)

Error: no slot of name "Pr" for this object of class "summary.aov"
In addition: Warning message:
trying to get slot "Pr" from an object (class "summary.aov") that is
not an S4 object
>

Research Programmer & Technician
The Scripps Research Institute
Mass Spectrometry Core Facility
  o The
 /
o Scripps
 \
  o Research
 /
o Institute



[EMAIL PROTECTED] wrote:
>  Hi: it's probably the Pr(> F) element so just access it by
>
> sum<-summary(whatever).
>
> then sum$Pr(>F) will probably work. But make sure that's it because 
> usually the name is pval or pvalue etc so I'm
> surprised about the weird name.
>
>
>
>
> On Wed, May 7, 2008 at  8:47 PM, Paul Benton wrote:
>
>> hello all,
>>
>> Quick question, how do I get the p value out of the anova?
>>
>> Thanks,
>>
>> Paul
>>
>>> pb<-aov(as.numeric(diff[5,16:33]) ~ grF)
>>> summary(pb)
>> Df Sum SqMean Sq F value  Pr(>F)
>> grF  3 2.7860e+10 9.2867e+09  4.2236 0.02534 *
>> Residuals   14 3.0783e+10 2.1988e+09
>> ---
>> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>> str(summary(pb))
>> List of 1
>>  $ :Classes 'anova' and 'data.frame':   2 obs. of  5 variables:
>>   ..$ Df : num [1:2] 3 14
>>   ..$ Sum Sq : num [1:2] 2.79e+10 3.08e+10
>>   ..$ Mean Sq: num [1:2] 9.29e+09 2.20e+09
>>   ..$ F value: num [1:2] 4.22   NA
>>   ..$ Pr(>F) : num [1:2] 0.0253 NA
>>  - attr(*, "class")= chr [1:2] "summary.aov" "listof"
>>> attr(summary(pb), "Pr(>F)")
>> NULL
>>
>> __
>> R-help@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.
>
>

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] anova p value extraction

[Jorge Ivan Velez]

Hi Paul,

Perhaps

# Data set
set.seed(123)
x=rnorm(100)
y=x+2*rnorm(100)

# ANOVA
AOV=anova(lm (y ~ x))
p1=AOV$"Pr(>F)"[1]
p1

or

# ANOVA 2
AOV1=aov(y ~ x)
p2=unlist(summary(AOV1))["Pr(>F)1"]
p2


HTH,

Jorge



On Wed, May 7, 2008 at 8:47 PM, Paul Benton <[EMAIL PROTECTED]> wrote:

> hello all,
>
> Quick question, how do I get the p value out of the anova?
>
> Thanks,
>
> Paul
>
> > pb<-aov(as.numeric(diff[5,16:33]) ~ grF)
> > summary(pb)
>Df Sum SqMean Sq F value  Pr(>F)
> grF  3 2.7860e+10 9.2867e+09  4.2236 0.02534 *
> Residuals   14 3.0783e+10 2.1988e+09
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > str(summary(pb))
> List of 1
>  $ :Classes 'anova' and 'data.frame':   2 obs. of  5 variables:
>  ..$ Df : num [1:2] 3 14
>  ..$ Sum Sq : num [1:2] 2.79e+10 3.08e+10
>  ..$ Mean Sq: num [1:2] 9.29e+09 2.20e+09
>  ..$ F value: num [1:2] 4.22   NA
>  ..$ Pr(>F) : num [1:2] 0.0253 NA
>  - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> > attr(summary(pb), "Pr(>F)")
> NULL
>
> __
> R-help@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.

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] anova p value extraction

[Abhijit Dasgupta]

What works, amazingly, is
summary(pb)$"Pr(>F)" [note the quotes]

Abhijit

On Wed, 07 May 2008 18:17:09 -0700
"H. Paul Benton" <[EMAIL PROTECTED]> wrote:

> Yea the anova object seems to be odd. It's not S4 so that's why I tried 
> originally the attr() funtion but
> 
> summary(pb)$Pr(>F)
> 
> Error: unexpected '>' in "summary(pb)$Pr(>"
> > summary(pb)$Pr
> 
> NULL
> > summary(pb)@Pr(>F)
> 
> Error: unexpected '>' in "summary(pb)@Pr(>"
> > summary(pb)@Pr(F)
> 
> Error: no slot of name "Pr" for this object of class "summary.aov"
> In addition: Warning message:
> trying to get slot "Pr" from an object (class "summary.aov") that is
> not an S4 object
> >
> 
> Research Programmer & Technician
> The Scripps Research Institute
> Mass Spectrometry Core Facility
>   o The
>  /
> o Scripps
>  \
>   o Research
>  /
> o Institute
> 
> 
> 
> [EMAIL PROTECTED] wrote:
> >  Hi: it's probably the Pr(> F) element so just access it by
> >
> > sum<-summary(whatever).
> >
> > then sum$Pr(>F) will probably work. But make sure that's it because 
> > usually the name is pval or pvalue etc so I'm
> > surprised about the weird name.
> >
> >
> >
> >
> > On Wed, May 7, 2008 at  8:47 PM, Paul Benton wrote:
> >
> >> hello all,
> >>
> >> Quick question, how do I get the p value out of the anova?
> >>
> >> Thanks,
> >>
> >> Paul
> >>
> >>> pb<-aov(as.numeric(diff[5,16:33]) ~ grF)
> >>> summary(pb)
> >> Df Sum SqMean Sq F value  Pr(>F)
> >> grF  3 2.7860e+10 9.2867e+09  4.2236 0.02534 *
> >> Residuals   14 3.0783e+10 2.1988e+09
> >> ---
> >> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >>> str(summary(pb))
> >> List of 1
> >>  $ :Classes 'anova' and 'data.frame':   2 obs. of  5 variables:
> >>   ..$ Df : num [1:2] 3 14
> >>   ..$ Sum Sq : num [1:2] 2.79e+10 3.08e+10
> >>   ..$ Mean Sq: num [1:2] 9.29e+09 2.20e+09
> >>   ..$ F value: num [1:2] 4.22   NA
> >>   ..$ Pr(>F) : num [1:2] 0.0253 NA
> >>  - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> >>> attr(summary(pb), "Pr(>F)")
> >> NULL
> >>
> >> __
> >> R-help@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.
> >
> >
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@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.


-- 
Abhijit Dasgupta

__
R-help@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.

Re: [R] anova p value extraction

[H. Paul Benton]

yep That worked Jorge. Thanks!

summary(pb)$"Pr(>F)1"
NULL
> summ<-summary(pb)
> class(unlist(summ))

[1] "numeric"
> unlist(summ)["Pr(>F)1"]

   Pr(>F)1
0.02533637
> names(summ)
NULL 



Jorge Ivan Velez wrote:
>
> Hi Paul,
>
> Perhaps
>
> # Data set
> set.seed(123)
> x=rnorm(100)
> y=x+2*rnorm(100)
>
> # ANOVA
> AOV=anova(lm (y ~ x))
> p1=AOV$"Pr(>F)"[1]
> p1
>  
> or
>
> # ANOVA 2
> AOV1=aov(y ~ x)
> p2=unlist(summary(AOV1))["Pr(>F)1"]
> p2
>
>
> HTH,
>
> Jorge
>
>
>
> On Wed, May 7, 2008 at 8:47 PM, Paul Benton <[EMAIL PROTECTED] 
> > wrote:
>
> hello all,
>
> Quick question, how do I get the p value out of the anova?
>
> Thanks,
>
> Paul
>
> > pb<-aov(as.numeric(diff[5,16:33]) ~ grF)
> > summary(pb)
>Df Sum SqMean Sq F value  Pr(>F)
> grF  3 2.7860e+10 9.2867e+09  4.2236 0.02534 *
> Residuals   14 3.0783e+10 2.1988e+09
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > str(summary(pb))
> List of 1
>  $ :Classes 'anova' and 'data.frame':   2 obs. of  5 variables:
>  ..$ Df : num [1:2] 3 14
>  ..$ Sum Sq : num [1:2] 2.79e+10 3.08e+10
>  ..$ Mean Sq: num [1:2] 9.29e+09 2.20e+09
>  ..$ F value: num [1:2] 4.22   NA
>  ..$ Pr(>F) : num [1:2] 0.0253 NA
>  - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> > attr(summary(pb), "Pr(>F)")
> NULL
>
> __
> R-help@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.
>

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] anova p value extraction

[Gavin Simpson]

quote the variable name or index the results from anova as a dataframe
using [. Someone (Prof. Ripley IIRC, apologies if I got this wrong) once
told me that backticks ` are the preferred, portable way of doing this,
but in this case " quotes work as well.

> example(anova.lm) ## produces fit
> tmp <- anova(fit)
> str(tmp) ## notice the Classes line below
Classes ‘anova’ and 'data.frame':   5 obs. of  5 variables:
 $ Df : int  1 1 1 1 45
 $ Sum Sq : num  204.1  53.3  12.4  63.1 650.7
 $ Mean Sq: num  204.1  53.3  12.4  63.1  14.5
 $ F value: num  14.116  3.689  0.858  4.360 NA
 $ Pr(>F) : num  0.000492 0.061125 0.359355 0.042471   NA
 - attr(*, "heading")= chr  "Analysis of Variance Table\n" "Response: sr"
> tmp$`Pr(>F)`
[1] 0.0004921955 0.0611254598 0.3593550848 0.0424711387   NA
> tmp$"Pr(>F)"
[1] 0.0004921955 0.0611254598 0.3593550848 0.0424711387   NA
> tmp[,5]
[1] 0.0004921955 0.0611254598 0.3593550848 0.0424711387   NA

##drop the last entry corresponding to the row for Residuals:
> tmp[1:4,5]
[1] 0.0004921955 0.0611254598 0.3593550848 0.0424711387

This has nothing to do with S3 or S4, just knowing how to reference
non-standard (probably not the correct terminology) names.

HTH

G

On Wed, 2008-05-07 at 18:17 -0700, H. Paul Benton wrote:
> Yea the anova object seems to be odd. It's not S4 so that's why I
> tried originally the attr() funtion but
> 
> summary(pb)$Pr(>F)
> 
> Error: unexpected '>' in "summary(pb)$Pr(>"
> > summary(pb)$Pr
> 
> NULL
> > summary(pb)@Pr(>F)
> 
> Error: unexpected '>' in "summary(pb)@Pr(>"
> > summary(pb)@Pr(F)
> 
> Error: no slot of name "Pr" for this object of class "summary.aov"
> In addition: Warning message:
> trying to get slot "Pr" from an object (class "summary.aov") that is
> not an S4 object
> >
> 
> Research Programmer & Technician
> The Scripps Research Institute
> Mass Spectrometry Core Facility
>   o The
>  /
> o Scripps
>  \
>   o Research
>  /
> o Institute
> 
> 
> 
> [EMAIL PROTECTED] wrote:
> >  Hi: it's probably the Pr(> F) element so just access it by
> >
> > sum<-summary(whatever).
> >
> > then sum$Pr(>F) will probably work. But make sure that's it because 
> > usually the name is pval or pvalue etc so I'm
> > surprised about the weird name.
> >
> >
> >
> >
> > On Wed, May 7, 2008 at  8:47 PM, Paul Benton wrote:
> >
> >> hello all,
> >>
> >> Quick question, how do I get the p value out of the anova?
> >>
> >> Thanks,
> >>
> >> Paul
> >>
> >>> pb<-aov(as.numeric(diff[5,16:33]) ~ grF)
> >>> summary(pb)
> >> Df Sum SqMean Sq F value  Pr(>F)
> >> grF  3 2.7860e+10 9.2867e+09  4.2236 0.02534 *
> >> Residuals   14 3.0783e+10 2.1988e+09
> >> ---
> >> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >>> str(summary(pb))
> >> List of 1
> >>  $ :Classes 'anova' and 'data.frame':   2 obs. of  5 variables:
> >>   ..$ Df : num [1:2] 3 14
> >>   ..$ Sum Sq : num [1:2] 2.79e+10 3.08e+10
> >>   ..$ Mean Sq: num [1:2] 9.29e+09 2.20e+09
> >>   ..$ F value: num [1:2] 4.22   NA
> >>   ..$ Pr(>F) : num [1:2] 0.0253 NA
> >>  - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> >>> attr(summary(pb), "Pr(>F)")
> >> NULL
> >>
> >> __
> >> R-help@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.
> >
> >
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@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.
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

__
R-help@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.

Re: [R] ANOVA between linear models.

[Rune Haubo]

Hi Richard

You are trying to compare two models, that are not nested. This means
that all usual asymptotics of the test statistics break down, hence
the (second) test you are attempting is not meaningful. Usually one
decides on the form of the response on other grounds such as residual
analysis or based on a box-cox analysis. You could take a look at
boxcox() in MASS (book as well as package).

That said it *is* possible to compare completely different models
using the likelihood or AIC in case of difference in number of
parameters although tests are seldom (if ever) available. In this case
you need to be sure to include the additive constants of the log
likelihood - see ?extractAIC.

In case the log-transform makes sense for your problem, you might want
to contemplate a gamma-GLM with a log link.

Hope this helps

/Rune

2008/5/15 Richard Martin <[EMAIL PROTECTED]>:
> Hi All,
>
> I'm accustomed to performing an ANOVA to aid in choosing between
> linear models (for example y~x or y~x+x^2), however with different
> models I can't seem to do it. I'm trying to fit an exponential model
> of the form ye^(bt).
>
> Below is a code snippet that highlights what I'm trying to do
>
> s = rnorm(100, sd = 1, mean=10)
> s = s + seq(0.1,10,0.1)
> x = 1:length(s)
> model.lin = lm(s ~ x)
> model.poly = lm(s ~x  + I(x^2))
> model.exp = lm(log(s) ~ x)
>
> anova(model.lin, model.poly)
> #gives the correct outcome
>
> anova(model.lin, model.exp)
> #doesn't.
>
> This fails because of the transformation of the response variable. Can
> anyone give any advice as to how I should proceed - is there a
> different test to use in this instance, or some way of reversing the
> transform after the model has been fitted?
>
> Any help greatly appreciated!!
>
> Richard
>
> --
> Contendere, Explorare, Invenire, et non Cedere
>
> __
> R-help@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.
>

__
R-help@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.

[R] Anova of a nls object

[Eric Ferreira]

Dear useRs,

How can I perform the analysis of variance of just one nls object?

Regards


-- 
Eric B Ferreira
Exact Sciences Department
Federal University of Lavras
Minas Gerais - Brazil

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova function and glm.nb

[Prof Brian Ripley]

The re-estimated version is the principled one -- the other is a 
computational shortcut.

However, what you model might choose depends on why you are doing model 
selection, and also if there is a subject-specific reason to consider the 
terms in this order.

And if you look in MASS (the book) you will see a similar example and 
discussion of how to treat it.  It uses dropterm(), which is likely to be 
preferable to either of these approaches.

On Mon, 7 Apr 2008, Nelson, Gary (FWE) wrote:

> Hi All,
>
> I am using the glm.nb function from the MASS package (current version)
> to fit and compare GLMs with negative binomial error distributions.  My
> question is: what is the appropriate method to use in the anova function
> to compare models? If only one fitted object like
>
> m<-glm.nb(number<-p+sal+temp,data=data)
>
> is specified in the anova function (anova(m)), a fixed theta is used to
> generate the analysis of deviance:
>
> Analysis of Deviance Table
> Model: Negative Binomial(0.2345), link: log
> Response: number
> Terms added sequentially (first to last)
>
>  Df Deviance Resid. Df Resid. Dev P(>|Chi|)
> NULL117122.707
> p  1   11.327   116111.380 0.001
> sal12.286   115109.094 0.131
> tem11.979   114107.115 0.159
> ph 12.567   113104.549 0.109
> Warning message:
> In anova.negbin(m) : tests made without re-estimating 'theta'
>
>
> If multiple fitted objects like
>
> m1<-glm.nb(number~1,data=data)
> m2<-glm.nb(number~p,data=data)
> m3<-glm.nb(number~p+sal,data=data)
> m4<-glm.nb(number~p+sal+temp,data=data)
>
>
> is specified (anova(m1,m2,m3,4)), the theta is assumed re-estimated in
> each case to calculate the likelihood ratio tests:
>
> Likelihood ratio tests of Negative Binomial Models
> Response: number
>   Model theta Resid. df2 x log-lik.   Testdf LR
> stat. Pr(Chi)
> 1  1 0.1892056   117   -527.7463
>
> 2  p 0.2153105   116   -517.9349 1 vs 2 1
> 9.811382 0.001734351
> 3p + sal 0.2214626   115   -515.7942 2 vs 3 1
> 2.140706 0.143435894
> 4  p + sal + tem 0.2261900   114   -513.8846 3 vs 4 1
> 1.909643 0.167002884
> 5 p + sal + tem + ph 0.2344827   113   -511.3633 4 vs 5 1
> 2.521237 0.112322429
>
>
> The conclusions are the same, but I'd like to know if one method is
> favored over the other.
>
> Thanks,
>
> Gary Nelson.
>
> __
> R-help@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.
>

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@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.

[R] ANOVA within a GAMM framework

[Maik Eisenbeiß]

Hi R user, 

I am using the gamm() function of the mgcv-package. Now I would like to
decide on the random effects to include in the model. Within a GAMM
framework, is it allowed to compare the following two models 

inv_1<-gamm(y~te(sat,inv),data=daten_final, random=list(proband=~1))

inv_2<-gamm(y~te(sat,inv),data=daten_final, random=list(proband=~sat))

with a likelihood ratio test for a traditional GLMM, like this:

anova(inv_1$lme, inv_2$lme)

The output is as follows: 

  Model df  AIC  BIClogLik   Test  L.Ratio p-value
inv_2$lme 1 10 21495.90 21557.59 -10737.95
inv_1$lme 2  8 23211.12 23260.46 -11597.56 1 vs 2 1719.214  <.0001


Or is this not in tune with the automatic smoothing parameter selection
(i.e. it is not exactly the same for both model alternatives)?


If not, how can I decide on the selection of random effects?
   

Thanks in advance for your help.

Best regards 
Maik


Dipl.-Kfm. Maik Eisenbeiß 
Marketing Centrum Münster 
Institut für Anlagen und Systemtechnologien 
Westfälische Wilhelms-Universität Münster
Am Stadtgraben 13–15
48143 Münster 

Telefon: +49 251 83-29920
Telefax: +49 251 83-22903

E-Mail: [EMAIL PROTECTED]
Web: http://www.marketing-centrum.de/ias 

__
R-help@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.

Re: [R] Anova and unbalanced designs

[Tal Galili]

 type-III SS is
> correct using contr.treatment(), but the within SS is not. As is generally
> the case, to get reasonable type-III tests (i.e., tests of reasonable
> hypotheses), it's necessary to have contrasts that are orthogonal in the
> row-basis of the design, such as contr.sum(),  contr.helmert(), or
> contr.poly(). The "type-II" tests, however, are insensitive to the contrast
> parametrization. Anova() always uses an orthogonal parametrization for the
> within-subjects design.
>
> The general advice in ?Anova is, "Be very careful in formulating the model
> for type-III tests, or the hypotheses tested will not make sense."
>
> Thanks, Peter, for pointing this out.
>
> John
>
> ------
> John Fox, Professor
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> web: socserv.mcmaster.ca/jfox
>
>
> > -Original Message-
> > From: Peter Dalgaard [mailto:p.dalga...@biostat.ku.dk]
> > Sent: January-24-09 6:31 PM
> > To: Nils Skotara
> > Cc: John Fox; r-help@r-project.org; 'Michael Friendly'
> > Subject: Re: [R] Anova and unbalanced designs
> >
> > Nils Skotara wrote:
> > > Dear John,
> > >
> > > thank you again! You replicated the type III result I got in SPSS! When
> I
> > > calculate Anova() type II:
> > >
> > > Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
> > >
> > > SS num Df Error SS den Df  F  Pr(>F)
> > > between 4.8000  1   9.  8 4.2667 0.07273 .
> > > within  0.2000  1  10.6667  8 0.1500 0.70864
> > > between:within  2.1333  1  10.6667  8 1.6000 0.24150
> > > ---
> > > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > >
> > > I see the exact same values as you had written.
> > > However, and now I am really lost, type III (I did not change anything
> > else)
> > > leads to the following:
> > >
> > > Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
> > >
> > >   SS num Df Error SS den Df   F
> Pr(>F)
> > > (Intercept)   72.000  19.000  8 64.
> 4.367e-05
> > ***
> > > between4.800  19.000  8  4.2667
> 0.07273 .
> > > as.factor(within)  2.000  1   10.667  8  1.5000
> 0.25551
> > > between:as.factor(within)  2.133  1   10.667  8  1.6000
> 0.24150
> > > ---
> > > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > >
> > > How is this possible?
> >
> > This looks like a contrast parametrization issue: If we look at the
> > per-group mean within-differences and their SE, we get
> >
> >  > summary(lm(within1-within2~between - 1))
> > ..
> > Coefficients:
> >   Estimate Std. Error t value Pr(>|t|)
> > between1  -1. 0.8165  -1.2250.256
> > between2   0. 0.6667   0.5000.631
> > ..
> >  > table(between)
> > between
> > 1 2
> > 4 6
> >
> > Now, the type II F test is based on weighting the two means as you would
> > after testing for no interaction
> >
> >  > (4*-1+6*.)^2/(4^2*0.8165^2+6^2*0.6667^2)
> > [1] 0.1500205
> >
> > and type III is to weight them as if there had been equal counts
> >
> >  > (5*-1+5*.)^2/(5^2*0.8165^2+5^2*0.6667^2)
> > [1] 0.400022
> >
> > However, the result above corresponds to looking at group1 only
> >
> >  > (-1)^2/(0.8165^2)
> > [1] 1.499987
> >
> > It helps if you choose orhtogonal contrast parametrizations:
> >
> >  > options(contrasts=c("contr.sum","contr.helmert"))
> >  > betweenanova <- lm(values ~ between)> Anova(betweenanova, idata=with,
> > idesign= ~as.factor(within), type = "III" )
> >
> > Type III Repeated Measures MANOVA Tests: Pillai test statistic
> >Df test stat approx F num Df den DfPr(>F)
> > (Intercept)1 0.963  209.067  1  8 5.121e-07
> ***
> > between1 0.3484.267  1  8   0.07273 .
> > as.factor(within)  1 0.0480.400  1  8   0.54474
> > between:as.factor(within)  1 0.1671.600  1  8   0.24150
> > ---
> > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >
> >
> >
> >
> > --
> > O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
> >c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
> >   (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
> > ~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907
>
> __
> R-help@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.
>



-- 
--


My contact information:
Tal Galili
Phone number: 972-50-3373767
FaceBook: Tal Galili
My Blogs:
www.talgalili.com
www.biostatistics.co.il

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Anova and unbalanced designs

[John Fox]

Dear Tal,

> -Original Message-
> From: Tal Galili [mailto:tal.gal...@gmail.com]
> Sent: February-14-09 10:23 AM
> To: John Fox
> Cc: Peter Dalgaard; Nils Skotara; r-help@r-project.org; Michael Friendly
> Subject: Re: [R] Anova and unbalanced designs
> 
> Hello John and other R mailing list members.
> 
> I've been following your discussions regarding the Anova command for the
SS
> type 2/3 repeated measures Anova, and I have a question:
> 
> I found that when I go from using type II to using type III, the summary
> model is suddenly added with an "intercept" term (example in the end of
the
> e-mail). So my question is
> 1) why is this "intercept" term added (in SS type "III" vs the type "II")?

The computational approach taken in Anova() makes it simpler to include the
intercept in the "type-III" tests and not to include it in the "type-II"
tests.

> 2) Can/should this "intercept" term be removed ? (or how should it be
> interpreted ?)

The test for the intercept is rarely of interest. A "type-II" test for the
intercept would test that the unconditional mean of the response is 0; a
"type-III" test for the intercept would test that the constant term in the
full model fit to the data is 0. The latter depends upon the parametrization
of the model (in the case of an ANOVA model, what kind of "contrasts" are
used). You state that the example that you give is taken from ?Anova but
there's a crucial detail that's omitted: The help file only gives the
"type-II" tests; the "type-III" tests are also reasonable here, but they
depend upon having used "contr.sum" (or another set of contrasts that's
orthogonal in the row basis of the model matrix) for the between-subject
factors, treatment and gender. This detail is in the data set:

> OBrienKaiser$gender
 [1] M M M F F M M F F M M M F F F F
attr(,"contrasts")
[1] contr.sum
Levels: F M

> OBrienKaiser$treatment
 [1] control control control control control A   A   A   A
B   B  
[12] B   B   B   B   B  
attr(,"contrasts")
[,1] [,2]
control   -20
A  1   -1
B  11
Levels: control A B

With proper contrast coding, the "type-III" test for the intercept tests
that the mean of the cell means (the "grand mean") is 0.

Had the default dummy-coded contrasts (from contr.treatment) been used, the
tests would not have tested reasonable hypotheses. My advice, from the help
file: "Be very careful in formulating the model for type-III tests, or the
hypotheses tested will not make sense."

I hope this helps,
 John

> 
> My purpose is to be able to use the Anova for analyzing an experiment with
a
> 2 between and 3 within factors, where the between factors are not
balanced,
> and the within factors are (that is why I can't use the aov command).
> 
> 
> #---code start
> 
> #---code start
> 
> #---code start
> 
> # (taken from the ?Anova help file)
> 
> phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
> levels=c("pretest", "posttest", "followup"))
> hour <- ordered(rep(1:5, 3))
> idata <- data.frame(phase, hour)
> idata
> mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5,
>  post.1, post.2, post.3, post.4, post.5,
>  fup.1, fup.2, fup.3, fup.4, fup.5) ~
treatment*gender,
> data=OBrienKaiser)
> 
> # now we have two options
> # option one is to use type II:
> 
> (av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = "II"))
> 
> 
> #output:
> Type II Repeated Measures MANOVA Tests: Pillai test statistic
> Df test stat approx F num Df den DfPr(>F)
> treatment20.4809   4.6323  2 10 0.0376868
*
> gender   10.2036   2.5558  1 10 0.1409735
> treatment:gender 20.3635   2.8555  2 10 0.1044692
> phase10.8505  25.6053  2  9 0.0001930
***
> treatment:phase  20.6852   2.6056  4 20 0.0667354
.
> gender:phase 10.0431   0.2029  2  9 0.8199968
> treatment:gender:phase   20.3106   0.9193  4 20 0.4721498
> hour 10.9347  25.0401  4  7 0.0003043
***
> treatment:hour   20.3014   0.3549  8 16 0.9295212
> gender:hour  10.2927   0.7243  4  7 0.6023742
> treatment:gender:hour20.5702   0.7976  8 16 0.6131884
> phase:hour