Re: [R] contrast matrix for aov
On Wed, 9 Mar 2005, Darren Weber wrote: How do we specify a contrast interaction matrix for an ANOVA model? We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) Cue Direction (2) x Brain Hemisphere(2) Each of these has 2 levels, 'left' and 'right', so it's a simple 2x2 design matrix. We have 8 subjects in each cell (a balanced design) and we want to specify the interaction contrast so that: CueLeftCueRght for the Right Hemisphere CueRghtCueLeft for the Left Hemisphere. Here is a copy of the relevant commands for R: lh_cueL - rowMeans( LHroi.cueL[,t1:t2] ) lh_cueR - rowMeans( LHroi.cueR[,t1:t2] ) rh_cueL - rowMeans( RHroi.cueL[,t1:t2] ) rh_cueR - rowMeans( RHroi.cueR[,t1:t2] ) roiValues - c( lh_cueL, lh_cueR, rh_cueL, rh_cueR ) cuelabels - c(CueLeft, CueRight) hemlabels - c(LH, RH) roiDataframe - data.frame( roi=roiValues, Subject=gl(8,1,32,subjectlabels), Hemisphere=gl(2,16,32,hemlabels), Cue=gl(2,8,32,cuelabels) ) roi.aov - aov(roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), data=roiDataframe) I think the error model should be Error(Subject). In what sense are `Cue' and `Cue:Hemisphere' random effects nested inside `Subject'? Let me fake some `data': set.seed(1); roiValues - rnorm(32) subjectlabels - paste(V1:8, sep = ) options(contrasts = c(contr.helmert, contr.poly)) roi.aov - aov(roi ~ Cue*Hemisphere + Error(Subject), data=roiDataframe) roi.aov Call: aov(formula = roi ~ Cue * Hemisphere + Error(Subject), data = roiDataframe) Grand Mean: 0.1165512 Stratum 1: Subject Terms: Residuals Sum of Squares 4.200946 Deg. of Freedom 7 Residual standard error: 0.7746839 Stratum 2: Within Terms: Cue Hemisphere Cue:Hemisphere Residuals Sum of Squares 0.216453 0.019712 0.057860 21.896872 Deg. of Freedom 1 1 121 Residual standard error: 1.021131 Estimated effects are balanced Note that all the action is in one stratum, and the SSQs are the same as aov(roi ~ Subject + Cue * Hemisphere, data = roiDataframe) (and also the same as for your fit). print(summary(roi.aov)) It auto-prints, so you don't need print(). I've tried to create a contrast matrix like this: cm - contrasts(roiDataframe$Cue) which gives the main effect contrasts for the Cue factor. I really want to specify the interaction contrasts, so I tried this: # c( lh_cueL, lh_cueR, rh_cueL, rh_cueR ) # CueRightCueLeft for the Left Hemisphere. # CueLeftCueRight for the Right Hemisphere cm - c(-1, 1, 1, -1) dim(cm) - c(2,2) (That is up to sign what Helmert contrasts give you.) roi.aov - aov( roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), contrasts=cm, data=roiDataframe) print(summary(roi.aov)) but the results of these two aov commands are identical. Is it the case that the 2x2 design matrix is always going to give the same F values for the interaction regardless of the contrast direction? Yes, as however you code the design (via `contrasts') you are fitting the same subspaces. Not sure what you mean by `contrast direction', though. However, you have not specified `contrasts' correctly: contrasts: A list of contrasts to be used for some of the factors in the formula. and cm is not a list, and an interaction is not a factor. OR, is there some way to get a summary output for the contrasts that is not available from the print method? For more than two levels, yes: see `split' under ?summary.aov. Also, see se.contrasts which allows you to find the standard error for any contrast. For the fixed-effects model you can use summary.lm: fit - aov(roi ~ Subject + Cue * Hemisphere, data = roiDataframe) summary(fit) Df Sum Sq Mean Sq F value Pr(F) Subject 7 4.2009 0.6001 0.5756 0.7677 Cue 1 0.2165 0.2165 0.2076 0.6533 Hemisphere 1 0.0197 0.0197 0.0189 0.8920 Cue:Hemisphere 1 0.0579 0.0579 0.0555 0.8161 Residuals 21 21.8969 1.0427 summary.lm(fit) Call: aov(formula = roi ~ Subject + Cue * Hemisphere, data = roiDataframe) Residuals: Min 1Q Median 3Q Max -1.7893 -0.4197 0.1723 0.5868 1.3033 Coefficients: Estimate Std. Error t value Pr(|t|) [...] Cue1 -0.082240.18051 -0.4560.653 Hemisphere1 0.024820.18051 0.1370.892 Cue1:Hemisphere1 -0.042520.18051 -0.2360.816 where the F values are the squares of the t values. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44
Re: [R] contrast matrix for aov
Prof Brian Ripley [EMAIL PROTECTED] writes: On Wed, 9 Mar 2005, Darren Weber wrote: How do we specify a contrast interaction matrix for an ANOVA model? We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) Actually, that's not usual in SAS (and not SPSS either, I believe) in things like proc glm; model y1-y4= ; repeated row 2 col 2; [Not that SAS/SPSS is the Gospel, but they do tend to set the terminology in these matters.] There you'd get the analysis split up as analyses of three contrasts corresponding to the main effects and interaction, c(-1,-1,1,1), c(-1,1,-1,1), and c(-1,1,1,-1) in the 2x2 case (up to scale and sign). In the 2x2 case, this corresponds exactly to the 4-stratum model row*col + Error(block/(row*col)). (It is interesting to note that it is still not the optimal analysis for arbitrary covariance patterns because dependence between contrasts is not utilized - it is basically assumed to be absent.) With more than two levels, you get the additional complications of sphericity testing and GG/HF epsilons and all that. -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] contrast matrix for aov
Prof Brian Ripley wrote: On Wed, 9 Mar 2005, Darren Weber wrote: We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) roi.aov - aov(roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), data=roiDataframe) I think the error model should be Error(Subject). In what sense are `Cue' and `Cue:Hemisphere' random effects nested inside `Subject'? I do not understand this, and I think I am probably not the only one. That is why I would be grateful if you could give a bit more information. My understanding is that the fixed factors Cue and Hemisphere are crossed with the random factor Subject (in other words, Cue and Hemisphere are within-subjects factors, and this is probably why Darren called it a repeated measure design). In this case, it seems to me from the various textbooks I read on Anova, that the appropriate MS to test the interaction Cue:Hemisphere is Subject:Cue:Hemisphere (with 7 degress of freedom, as there are 8 independent subjects). If you input Error(Subject/(Cue*Hemisphere)) in the aov formula, then the test for the interaction indeed uses the Subject:Cue:Hemisphere source of variation in demoninator. This fits with the ouput of other softwares. If you include only 'Subjet', then the test for the interaction has 21 degrees of Freedom, and I do not understand what this tests. I apologize in if my terminology is not accurate. But I hope you can clarify what is wrong with the Error(Subject/(Cue*Hemisphere)) term, or maybe just point us to the relevant textbooks. Thanks in advance, Christophe Pallier Let me fake some `data': set.seed(1); roiValues - rnorm(32) subjectlabels - paste(V1:8, sep = ) options(contrasts = c(contr.helmert, contr.poly)) roi.aov - aov(roi ~ Cue*Hemisphere + Error(Subject), data=roiDataframe) roi.aov Call: aov(formula = roi ~ Cue * Hemisphere + Error(Subject), data = roiDataframe) Grand Mean: 0.1165512 Stratum 1: Subject Terms: Residuals Sum of Squares 4.200946 Deg. of Freedom 7 Residual standard error: 0.7746839 Stratum 2: Within Terms: Cue Hemisphere Cue:Hemisphere Residuals Sum of Squares 0.216453 0.019712 0.057860 21.896872 Deg. of Freedom 1 1 121 Residual standard error: 1.021131 Estimated effects are balanced Note that all the action is in one stratum, and the SSQs are the same as aov(roi ~ Subject + Cue * Hemisphere, data = roiDataframe) (and also the same as for your fit). print(summary(roi.aov)) It auto-prints, so you don't need print(). I've tried to create a contrast matrix like this: cm - contrasts(roiDataframe$Cue) which gives the main effect contrasts for the Cue factor. I really want to specify the interaction contrasts, so I tried this: # c( lh_cueL, lh_cueR, rh_cueL, rh_cueR ) # CueRightCueLeft for the Left Hemisphere. # CueLeftCueRight for the Right Hemisphere cm - c(-1, 1, 1, -1) dim(cm) - c(2,2) (That is up to sign what Helmert contrasts give you.) roi.aov - aov( roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), contrasts=cm, data=roiDataframe) print(summary(roi.aov)) but the results of these two aov commands are identical. Is it the case that the 2x2 design matrix is always going to give the same F values for the interaction regardless of the contrast direction? Yes, as however you code the design (via `contrasts') you are fitting the same subspaces. Not sure what you mean by `contrast direction', though. However, you have not specified `contrasts' correctly: contrasts: A list of contrasts to be used for some of the factors in the formula. and cm is not a list, and an interaction is not a factor. OR, is there some way to get a summary output for the contrasts that is not available from the print method? For more than two levels, yes: see `split' under ?summary.aov. Also, see se.contrasts which allows you to find the standard error for any contrast. For the fixed-effects model you can use summary.lm: fit - aov(roi ~ Subject + Cue * Hemisphere, data = roiDataframe) summary(fit) Df Sum Sq Mean Sq F value Pr(F) Subject 7 4.2009 0.6001 0.5756 0.7677 Cue 1 0.2165 0.2165 0.2076 0.6533 Hemisphere 1 0.0197 0.0197 0.0189 0.8920 Cue:Hemisphere 1 0.0579 0.0579 0.0555 0.8161 Residuals 21 21.8969 1.0427 summary.lm(fit) Call: aov(formula = roi ~ Subject + Cue * Hemisphere, data = roiDataframe) Residuals: Min 1Q Median 3Q Max -1.7893 -0.4197 0.1723 0.5868 1.3033 Coefficients: Estimate Std. Error t value Pr(|t|) [...] Cue1
Re: [R] contrast matrix for aov
On Thu, 10 Mar 2005, Christophe Pallier wrote: Prof Brian Ripley wrote: On Wed, 9 Mar 2005, Darren Weber wrote: We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) roi.aov - aov(roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), data=roiDataframe) I think the error model should be Error(Subject). In what sense are `Cue' and `Cue:Hemisphere' random effects nested inside `Subject'? I do not understand this, and I think I am probably not the only one. That is why I would be grateful if you could give a bit more information. My understanding is that the fixed factors Cue and Hemisphere are crossed with the random factor Subject (in other words, Cue and Hemisphere are within-subjects factors, and this is probably why Darren called it a repeated measure design). The issue is whether the variance of the error really depends on the treatment combination, which is what the Error(Subject/(Cue*Hemisphere)) assumes. With that model Error: Subject:Cue Df Sum Sq Mean Sq F value Pr(F) Cue1 0.2165 0.2165 0.1967 0.6708 Residuals 7 7.7041 1.1006 Error: Subject:Hemisphere Df Sum Sq Mean Sq F value Pr(F) Hemisphere 1 0.0197 0.0197 0.0154 0.9047 Residuals 7 8.9561 1.2794 Error: Subject:Cue:Hemisphere Df Sum Sq Mean Sq F value Pr(F) Cue:Hemisphere 1 0.0579 0.0579 0.0773 0.789 Residuals 7 5.2366 0.7481 you are assuming different variances for three contrasts. In this case, it seems to me from the various textbooks I read on Anova, that the appropriate MS to test the interaction Cue:Hemisphere is Subject:Cue:Hemisphere (with 7 degress of freedom, as there are 8 independent subjects). If you input Error(Subject/(Cue*Hemisphere)) in the aov formula, then the test for the interaction indeed uses the Subject:Cue:Hemisphere source of variation in demoninator. This fits with the ouput of other softwares. If you include only 'Subjet', then the test for the interaction has 21 degrees of Freedom, and I do not understand what this tests. It uses a common variance for all treatment combinations. I apologize in if my terminology is not accurate. But I hope you can clarify what is wrong with the Error(Subject/(Cue*Hemisphere)) term, or maybe just point us to the relevant textbooks. Nothing is `wrong' with it, it just seems discordant with the description of the experiment. -- 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] contrast matrix for aov
On Thu, 10 Mar 2005, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: On Wed, 9 Mar 2005, Darren Weber wrote: How do we specify a contrast interaction matrix for an ANOVA model? We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) Actually, that's not usual in SAS (and not SPSS either, I believe) in things like proc glm; model y1-y4= ; repeated row 2 col 2; [Not that SAS/SPSS is the Gospel, but they do tend to set the terminology in these matters.] That seems to be appropriate only if the four treatments are done in a particular order (`repeated') and one expects correlations in the responses. However, here the measurements seem to have been averages of replications. It may be usual to (mis?)specify experiments in SAS that way: I don't know what end users do, but it is not the only way possible in SAS. There you'd get the analysis split up as analyses of three contrasts corresponding to the main effects and interaction, c(-1,-1,1,1), c(-1,1,-1,1), and c(-1,1,1,-1) in the 2x2 case (up to scale and sign). In the 2x2 case, this corresponds exactly to the 4-stratum model row*col + Error(block/(row*col)). (It is interesting to note that it is still not the optimal analysis for arbitrary covariance patterns because dependence between contrasts is not utilized - it is basically assumed to be absent.) It also assumes that there is a difference between variances of the contrasts, that is there is either correlation between results or a difference in variances under different treatments. Nothing in the description led me to expect either of those, but I was asking why it was specified that way. (If the variance does differ with the mean then there are probably more appropriate analyses.) -- 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] contrast matrix for aov
Prof Brian Ripley wrote: On Wed, 9 Mar 2005, Darren Weber wrote: How do we specify a contrast interaction matrix for an ANOVA model? We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) Each subject experiences multiple instances of a visual stimulus that cues them to orient their attention to the left or right visual field. For each instance, we acquire brain activity measures from the left and right hemisphere (to put it simply). For each condition in this 2x2 design, we actually have about 400 instances for each cell for each subject. These are averaged in several ways, so we obtain a single scalar value for each subject in the final analysis. So, at this point, it does not appear to be a repeated measures analysis. Perhaps another way to put it - each subject experiences every cell of the design, so we have within-subjects factors rather than between-subjects factors. That is, each subject experiences all experimental conditions, we do not separate and randomly allocate subjects to only one cell or another of the experiment. I hope this clarifies the first question. I will try to digest further points in this thread, if they are not too confounded by this first point of contention already. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] contrast matrix for aov
As an R newbie (formerly SPSS), I was pleased to find some helpful notes on ANOVA here: http://personality-project.org/r/r.anova.html In my case, I believe the relevant section is: Example 4. Two-Way Within-Subjects ANOVA This is where I noted and copied the error notation. Sorry for any confusion about terms - I did mean within-subjects factors, rather than repeated measures (although, as noted earlier, we do have both in this experiment). Prof Brian Ripley wrote: On Thu, 10 Mar 2005, Christophe Pallier wrote: Prof Brian Ripley wrote: On Wed, 9 Mar 2005, Darren Weber wrote: We have a two-factor, repeated measures design, with Where does `repeated measures' come into this? You appear to have repeated a 2x2 experiment in each of 8 blocks (subjects). Such a design is usually analysed with fixed effects. (Perhaps you averaged over repeats in the first few lines of your code?) roi.aov - aov(roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), data=roiDataframe) I think the error model should be Error(Subject). In what sense are `Cue' and `Cue:Hemisphere' random effects nested inside `Subject'? I do not understand this, and I think I am probably not the only one. That is why I would be grateful if you could give a bit more information. My understanding is that the fixed factors Cue and Hemisphere are crossed with the random factor Subject (in other words, Cue and Hemisphere are within-subjects factors, and this is probably why Darren called it a repeated measure design). The issue is whether the variance of the error really depends on the treatment combination, which is what the Error(Subject/(Cue*Hemisphere)) assumes. With that model Error: Subject:Cue Df Sum Sq Mean Sq F value Pr(F) Cue1 0.2165 0.2165 0.1967 0.6708 Residuals 7 7.7041 1.1006 Error: Subject:Hemisphere Df Sum Sq Mean Sq F value Pr(F) Hemisphere 1 0.0197 0.0197 0.0154 0.9047 Residuals 7 8.9561 1.2794 Error: Subject:Cue:Hemisphere Df Sum Sq Mean Sq F value Pr(F) Cue:Hemisphere 1 0.0579 0.0579 0.0773 0.789 Residuals 7 5.2366 0.7481 you are assuming different variances for three contrasts. In this case, it seems to me from the various textbooks I read on Anova, that the appropriate MS to test the interaction Cue:Hemisphere is Subject:Cue:Hemisphere (with 7 degress of freedom, as there are 8 independent subjects). If you input Error(Subject/(Cue*Hemisphere)) in the aov formula, then the test for the interaction indeed uses the Subject:Cue:Hemisphere source of variation in demoninator. This fits with the ouput of other softwares. If you include only 'Subjet', then the test for the interaction has 21 degrees of Freedom, and I do not understand what this tests. It uses a common variance for all treatment combinations. I apologize in if my terminology is not accurate. But I hope you can clarify what is wrong with the Error(Subject/(Cue*Hemisphere)) term, or maybe just point us to the relevant textbooks. Nothing is `wrong' with it, it just seems discordant with the description of the experiment. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] contrast reference level
array chip [EMAIL PROTECTED] writes: suppose I have a factor with 4 levels: 'a','b','c','d'. I would like to do analysis of variance using aov() with the factor as independent variable. How can I specify the level b as the reference level just like the level a would be the reference level if using contr.treatment as the contrast? relevel() should do just that. -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Contrast matrices for nested factors
You should be able to get the behavior you want using the fit.glh() (short
for fit general linear hypothesis) function from the gregmisc/gmodels
package.
-G
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Behalf Of
Fernando Henrique
Ferraz P. da Rosa
Sent: Monday, September 06, 2004 11:02 PM
To: r-help
Subject: [R] Contrast matrices for nested factors
Hi, I'd like to know if it's possible to specify different
contrast matrices in lm() for a factor that is nested within
another one. This
is useful when we have a model where the nested factor has a different
number of levels, depending on the main factor.
Let me illustrate with an example to make it clearer. Consider
the following data set:
set.seed(1)
y - rnorm(14)
a - factor(c(rep(1,7),rep(2,3),rep(3,4)))
b - factor(c(1,1,1,2,2,3,3,1,1,2,1,1,2,2))
k - factor(c(1,2,3,1,2,1,2,1,2,1,1,2,1,2))
internal - data.frame(y,a,b,k)
Where y is an arbitrary response, a is a main factor, b is a
factor nested within a, and k is the replicate number. It is
easy to see
that depending on the level of a, b has different numbers of
levels. For
instance, when a = 1, we have that b might assume values 1, 2 or 3,
while a = 2 or 3, b might assume only 1 or 2.
I'd like then to use contrasts summing to 0, so I issue:
z - lm(y ~ a + a/b,data=internal,contrasts=list(a=contr.sum,
b=contr.sum))
The problem is, the design matrix is not quite what I
expected.
What happens is, instead of using a different contrast matrix for each
level of a where b is nested, it's using the same contrast matrix for
every b, namely:
contr.sum(3)
[,1] [,2]
110
201
3 -1 -1
So, when a=1, the columns of the design matrix are as
expected.
It sums to 0, because there are levels of b 1, 2 and 3, when a=1. But,
when a=2 or a=3, the same contrast matrix is being used, and then, the
factor effects do not sum to 0. That's obviously because there are no
values for b equal 3, when a != 1, and then the coding that
gets done is
'0' or '1'.
The design matrix lm() is creating is:
model.matrix(z)
(Intercept) a1 a2 a1:b1 a2:b1 a3:b1 a1:b2 a2:b2 a3:b2
11 1 0 1 0 0 0 0 0
21 1 0 1 0 0 0 0 0
31 1 0 1 0 0 0 0 0
41 1 0 0 0 0 1 0 0
51 1 0 0 0 0 1 0 0
61 1 0-1 0 0-1 0 0
71 1 0-1 0 0-1 0 0
81 0 1 0 1 0 0 0 0
91 0 1 0 1 0 0 0 0
10 1 0 1 0 0 0 0 1 0
11 1 -1 -1 0 0 1 0 0 0
12 1 -1 -1 0 0 1 0 0 0
13 1 -1 -1 0 0 0 0 0 1
14 1 -1 -1 0 0 0 0 0 1
What I would like to use is:
(Intercept) a1 a2 a1:b1 a2:b1 a3:b1 a1:b2
11 1 0 1 0 0 0 0 0
21 1 0 1 0 0 0 0 0
31 1 0 1 0 0 0 0 0
41 1 0 0 0 0 1 0 0
51 1 0 0 0 0 1 0 0
61 1 0-1 0 0-1 0 0
71 1 0-1 0 0-1 0 0
81 0 1 0 1 0 0 0 0
91 0 1 0 1 0 0 0 0
10 1 0 1 0-1 0 0 0 0
11 1 -1 -1 0 0 1 0 0 0
12 1 -1 -1 0 0 1 0 0 0
13 1 -1 -1 0 0-1 0 0 0
14 1 -1 -1 0 0-1 0 0 0
(notice that in the second matrix all collumns sum to
0, in the
first they don't).
Thank you,
--
Fernando Henrique Ferraz P. da Rosa
http://www.ime.usp.br/~feferraz
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] contrast lme and glmmPQL and getting additional results...
On 20 Mar 2004 at 11:06, Paul Johnson wrote: From what you say, it seems like you have a linear normal (mixed) model. This is what lme is made for, and there is no reason to use glmmPQL (which calls lme iteratively). Kjetil Halvorsen I have a longitudinal data analysis project. There are 10 observations on each of 15 units, and I'm estimating this with randomly varying intercepts along with an AR1 correction for the error terms within units. There is no correlation across units. Blundering around in R for a long time, I found that for linear/gaussian models, I can use either the MASS method glmmPQL (thanks to Venables and Ripley) or the lme from nlme (thanks to Pinheiro and Bates). (I also find that the package lme4 has GLMM, but I can't get the correlation structure to work with that, so I gave up on that one.) The glmmPQL and lme results are quite similar, but not identical. Here are my questions. 1. I believe that both of these offer consistent estimates. Does one have preferrable small sample properties? Is the lme the preferred method in this case because it is more narrowly designed to this gaussian family model with an identity link? If there's an argument in favor of PQL, I'd like to know it, because a couple of the Hypothesis tests based on t-statistics are affected. 2. Is there a pre-made method for calculation of the robust standard errors? I notice that model.matrix() command does not work for either lme or glmmPQL results, and so I start to wonder how people calculate sandwich estimators of the standard errors. 3. Are the AIC (or BIC) statistics comparable across models? Can one argue in favor of the glmmPQL results (with, say, a log link) if the AIC is more favorable than the AIC from an lme fit? In JK Lindsey's Models for Repeated Measurements, the AIC is sometimes used to make model selections, but I don't know where the limits of this application might lie. -- Paul E. Johnson email: [EMAIL PROTECTED] Dept. of Political Sciencehttp://lark.cc.ku.edu/~pauljohn 1541 Lilac Lane, Rm 504 University of Kansas Office: (785) 864-9086 Lawrence, Kansas 66044-3177 FAX: (785) 864-5700 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Contrast
Or ?contrast in the Design library.
-G
-Original Message-
From: Simon Blomberg [mailto:[EMAIL PROTECTED]
Sent: Tuesday, November 04, 2003 8:32 PM
To: Igor Roytberg; [EMAIL PROTECTED]
Subject: RE: [R] Contrast
see ?fit.contrast in library gregmisc.
Cheers,
Simon.
Simon Blomberg, PhD
Depression Anxiety Consumer Research Unit
Centre for Mental Health Research
Australian National University
http://www.anu.edu.au/cmhr/
[EMAIL PROTECTED] +61 (2) 6125 3379
-Original Message-
From: Igor Roytberg [mailto:[EMAIL PROTECTED]
Sent: Wednesday, 5 November 2003 12:04 PM
To: [EMAIL PROTECTED]
Subject: [R] Contrast
Could anyone please explain how to set up contrasts between
means in R. I want to know if before I conduct an experiment
and believe the mean for 1 and 2 will be different from means
3 and 4, Is this true? That is what I have to prove or
disprove, I thought that contrasts would be the way to go.
Thanks for the help.
Igor
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Contrast
see ?fit.contrast in library gregmisc. Cheers, Simon. Simon Blomberg, PhD Depression Anxiety Consumer Research Unit Centre for Mental Health Research Australian National University http://www.anu.edu.au/cmhr/ [EMAIL PROTECTED] +61 (2) 6125 3379 -Original Message- From: Igor Roytberg [mailto:[EMAIL PROTECTED] Sent: Wednesday, 5 November 2003 12:04 PM To: [EMAIL PROTECTED] Subject: [R] Contrast Could anyone please explain how to set up contrasts between means in R. I want to know if before I conduct an experiment and believe the mean for 1 and 2 will be different from means 3 and 4, Is this true? That is what I have to prove or disprove, I thought that contrasts would be the way to go. Thanks for the help. Igor [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Contrast specified with C() - R vs S-Plus problem
On Wed, 8 Oct 2003, Pascal A. Niklaus wrote: Hi, For a n-level factor, I'd like to specify the first contrast and have the remaining n-2 constructed automatically so that the set is orthogonal. I then test the contrasts with summary.lm(anova-object). In S-Plus, the following works: y.anova - aov( y ~ C(CO2,c(1,0,-1)) ) summary.lm(y.anova) I can't reproduce that in S-PLUS 6.1, and it is not as documented: contr what contrasts to use. May be one of four standard names ( helmert, poly, treatment, or sum), a function, or a matrix with as many rows as there are levels to the factor. You haven't given a matrix, and your vector gets converted to a 3x1 matrix when there are four levels. I suspect you have not got the same data in S-PLUS and R, but you haven't given us anything to check that. In R, it fails with the following error: levels(CO2) [1] A C E y.anova - aov(y + C(CO2,c(1,0,-1)) ) Error in contrasts-(*tmp*, value = contr) : wrong number of contrast matrix rows What is the way to do this in R? There are four levels, so CO2 - factor(c(, A, C, E)) attributes(C(CO2, as.matrix(1:4))) $levels [1] A C E $class [1] factor $contrasts [,1] [,2] [,3] 1 0.0236068 0.5472136 A2 -0.4393447 -0.7120227 C3 0.8078689 -0.2175955 E4 -0.3921311 0.3824045 does as you ask. -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Contrast specified with C() - R vs S-Plus problem
Thanks for the reply. Below is the solution and the S-Plus and R code that does the same (for documentation). I can't reproduce that in S-PLUS 6.1, and it is not as documented: In S-Plus 2000, C() complements the contrast matrix with orthogonal contrasts if only the first is given. CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) summary.aov(aov(m ~ C(CO2,c(1,0,-1 Df Sum of Sq Mean Sq F Value Pr(F) C(CO2, c(1, 0, -1)) 2 16.0 8.00 7.259259 0.004013532 Residuals 21 23.14286 1.102041 summary.lm(aov(m ~ C(CO2,c(1,0,-1 Coefficients: Value Std. Error t value Pr(|t|) (Intercept) 2.5000 0.214311.6667 0. C(CO2, c(1, 0, -1))1 -1. 0.2624-3.8103 0.0010 C(CO2, c(1, 0, -1))2 0. 0.3712 0. 1. Residual standard error: 1.05 on 21 degrees of freedom Multiple R-Squared: 0.4088 F-statistic: 7.259 on 2 and 21 degrees of freedom, the p-value is 0.004014 Correlation of Coefficients: (Intercept) C(CO2, c(1, 0, -1))1 C(CO2, c(1, 0, -1))1 0 C(CO2, c(1, 0, -1))2 0 0 In the meanwhile, I found that the same can be done in R like this: library(gregmisc) CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) contrastsCO2 - rbind(A vs E=c(1,0,-1)) a - aov(m ~ CO2, contrasts=list( CO2=make.contrasts(contrastsCO2))) summary(a) Df Sum Sq Mean Sq F value Pr(F) CO2 2 16.000 8.000 7.2593 0.004014 ** Residuals 21 23.143 1.102 summary.lm(a) Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 2.500e+00 2.143e-0111.67 1.22e-10 *** CO2A vs E -2.000e+00 5.249e-01-3.81 0.00102 ** CO2C29.774e-17 3.712e-01 2.63e-16 1.0 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 1.05 on 21 degrees of freedom Multiple R-Squared: 0.4088, Adjusted R-squared: 0.3525 F-statistic: 7.259 on 2 and 21 DF, p-value: 0.004014 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Contrast specified with C() - R vs S-Plus problem
On Wed, 8 Oct 2003, Pascal A. Niklaus wrote: Thanks for the reply. Below is the solution and the S-Plus and R code that does the same (for documentation). I can't reproduce that in S-PLUS 6.1, and it is not as documented: In S-Plus 2000, C() complements the contrast matrix with orthogonal contrasts if only the first is given. And so does R. You just did not specify it correctly in R. The code you quote works in R as well, but your original example had four levels, not three. CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) summary.aov(aov(m ~ C(CO2,c(1,0,-1 Df Sum Sq Mean Sq F value Pr(F) C(CO2, c(1, 0, -1)) 2 16.000 8.000 7.2593 0.004014 Residuals 21 23.143 1.102 Please do not make out your user errors to be S-R differences. Did you actually try your example in R? CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) summary.aov(aov(m ~ C(CO2,c(1,0,-1 Df Sum of Sq Mean Sq F Value Pr(F) C(CO2, c(1, 0, -1)) 2 16.0 8.00 7.259259 0.004013532 Residuals 21 23.14286 1.102041 summary.lm(aov(m ~ C(CO2,c(1,0,-1 Coefficients: Value Std. Error t value Pr(|t|) (Intercept) 2.5000 0.214311.6667 0. C(CO2, c(1, 0, -1))1 -1. 0.2624-3.8103 0.0010 C(CO2, c(1, 0, -1))2 0. 0.3712 0. 1. Residual standard error: 1.05 on 21 degrees of freedom Multiple R-Squared: 0.4088 F-statistic: 7.259 on 2 and 21 degrees of freedom, the p-value is 0.004014 Correlation of Coefficients: (Intercept) C(CO2, c(1, 0, -1))1 C(CO2, c(1, 0, -1))1 0 C(CO2, c(1, 0, -1))2 0 0 -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
