[R] path-diagram wirh Rgraphvis
Dear List,
how can I draw the following path diagram
A B C D E F
\ | / \ | /
G -- H
/ \/ \
I JK L
the problem I've got is that G and H need to be horizontally alingned
but the best I've done is diagonally or vertically alingned,
kind regards,
Rene
## Begin: R-code
#source(http://bioconductor.org/biocLite.R;)
#biocLite(Rgraphviz)
library(Rgraphviz)
library(graph)
# observed x
ox=c('RParAsp', 'RIQ', 'RSES','FSES', 'FIQ', 'FParAsp')
# observed y
oy=c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp')
# latend
l=c('FGenAsp','RGenAsp')
# pure nodes in a line
l1 - new(graphNEL, nodes = c(ox,oy,l), edgemode = directed)
plot(l1, dot)
l1 - addEdge(RParAsp, RGenAsp, l1)
l1 - addEdge(RIQ, RGenAsp, l1)
l1 - addEdge(RSES, RGenAsp, l1)
l1 - addEdge(FParAsp, FGenAsp, l1)
l1 - addEdge(FIQ, FGenAsp, l1)
l1 - addEdge(FSES, FGenAsp, l1)
l1 - addEdge(FGenAsp, FOccAsp, l1)
l1 - addEdge(FGenAsp, FEdAsp, l1)
l1 - addEdge(RGenAsp, ROccAsp, l1)
l1 - addEdge(RGenAsp, REdAsp, l1)
l1 - addEdge(RGenAsp, FGenAsp, l1)
l1 - addEdge(FGenAsp, RGenAsp, l1)
plot(l1, recipEdges=distinct)
sub1 - subGraph(ox, l1)
sub2 - subGraph(oy, l1)
sub3 - subGraph(l, l1)
sublist -
list( list(graph=sub1, cluster = T) ,
list(graph=sub2, cluster = T),
list(graph=sub3, cluster = T)
)
plot(l1,
subGList = sublist)
globalattrs$graph - list(rankdir = LR)
plot(l1, attrs = globalattrs)
## End: R-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] how to cumulate up times
Dear List,
given a vecor of times in 5,15 and 30 minutes and a start point
in time, lets say 09:30:00, how do I add up those times
to the start time getting a cumulative time sequence?
mt-times(c('00:05:00', '00:15:00', '00:30:00'))
mt wanted
00:05:00 09:35:00
00:15:00 09:50:00
00:30:00 10:20:00
Regards,
René
__
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] how to cumulate up times
Thanks Michael,
cumsum - yes of course!
Regards,
René
Zitat von R. Michael Weylandt michael.weyla...@gmail.com:
? cumsum
something like
library(chron) # Reporting packages you use is always considerate
mt - times(c('00:05:00', '00:15:00', '00:30:00')) # Spaces are legible!
times(09:30:00) + cumsum(mt)
Michael
On Tue, Apr 24, 2012 at 11:35 AM, René Mayer
ma...@psychologie.tu-dresden.de wrote:
Dear List,
given a vecor of times in 5,15 and 30 minutes and a start point
in time, lets say 09:30:00, how do I add up those times
to the start time getting a cumulative time sequence?
mt-times(c('00:05:00', '00:15:00', '00:30:00'))
mt wanted
00:05:00 09:35:00
00:15:00 09:50:00
00:30:00 10:20:00
Regards,
René
__
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] scatter3d: problem with spheres-color
Dear John and Duncan, thanks for your ideas! Unfortunatly, calling spheres from rgl did not resolve the problem on my machine. Both - spheres3d() and rgl.spheres() - behave the same: black spheres, all aqual colored. The only difference beeing the looking angle and thebackround color. Seems to me that the 'col=1:5' argument is completly ignored. René Zitat von Duncan Murdoch murdoch.dun...@gmail.com: On 12/04/2012 2:27 PM, John Fox wrote: Dear René, I've confirmed that the spheres aren't coloured correctly on my Ubuntu system (the first colour is used for all of the spheres), and I know that this works right on Windows, as you mentioned. I'm curious to try it on my Mac, but don't have that handy at the moment. I also looked at the code for scatter3d.default(), and that is pretty straightforward; scatterplot3d.default() draws the spheres with the command rgl.spheres(x, y, z, color = surface.col[as.numeric(groups)], radius = size) I'm copying this response to Duncan Murdoch (the coauthor and maintainer of the rgl package) in case he has any insight into the problem. Thank you for drawing this issue to my attention. Calling rgl.spheres looks dangerous to me: the rgl.* functions make permanent changes to material properties. Generally it's safer to call spheres3d, as all of the *3d versions of functions make local changes. But there should be no differences in that between Ubuntu and Windows. Can you put together a simple example that does give differences? For example, on Windows this gives 5 different colours: rgl.spheres(1:5, 1:5, 1:5, col=1:5, radius=(1:5)/10) My preferred version would be spheres3d(1:5, 1:5, 1:5, col=1:5, radius=(1:5)/10) Do they behave the same? Duncan Murdoch __ 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] scatter3d: problem with spheres-color
I get the rgb-a are the same rgb-a's: spheres3d(1:5, 1:5, 1:5, col=1:5, radius=(1:5)/10) rgl.ids() idtype 1 11 spheres rgl.attrib(11, color) r g b a [1,] 0 0.000 0 1 [2,] 1 0.000 0 1 [3,] 0 0.8039216 0 1 [4,] 0 0.000 1 1 [5,] 0 1.000 1 1 as all sheres are black it seems that only the first rgb-a row is used for all spheres. I remember that I could set different colors in different scatter3d calls, e.g., having all sheres 'red' or 'yellow' within one plot. As if only the first row is applied to all objects. René René Zitat von Duncan Murdoch murdoch.dun...@gmail.com: On 12-04-13 5:32 AM, René Mayer wrote: Dear John and Duncan, thanks for your ideas! Unfortunatly, calling spheres from rgl did not resolve the problem on my machine. Both - spheres3d() and rgl.spheres() - behave the same: black spheres, all aqual colored. The only difference beeing the looking angle and thebackround color. Seems to me that the 'col=1:5' argument is completly ignored. Could you try this: After using one of the commands that's not working (e.g. my spheres3d command from below), run the following: rgl.ids() This will give a dataframe of objects in the scene, something like idtype 1 6 spheres (There will be more things listed if you have plotted other objects, and the id number could be different.) Then run rgl.attrib(6, color) (The 6 is the id associated with the spheres object.) On my system after my spheres3d() call below, this gives rgl.attrib(6, color) r g b a [1,] 0 0.000 0 1 [2,] 1 0.000 0 1 [3,] 0 0.8039216 0 1 [4,] 0 0.000 1 1 [5,] 0 1.000 1 1 Can you tell me what you get? If you get the same thing as me, then you've got a rendering problem; if you only see black colors listed, then it's a problem at a higher level. Duncan MUrdoch René Zitat von Duncan Murdochmurdoch.dun...@gmail.com: On 12/04/2012 2:27 PM, John Fox wrote: Dear René, I've confirmed that the spheres aren't coloured correctly on my Ubuntu system (the first colour is used for all of the spheres), and I know that this works right on Windows, as you mentioned. I'm curious to try it on my Mac, but don't have that handy at the moment. I also looked at the code for scatter3d.default(), and that is pretty straightforward; scatterplot3d.default() draws the spheres with the command rgl.spheres(x, y, z, color = surface.col[as.numeric(groups)], radius = size) I'm copying this response to Duncan Murdoch (the coauthor and maintainer of the rgl package) in case he has any insight into the problem. Thank you for drawing this issue to my attention. Calling rgl.spheres looks dangerous to me: the rgl.* functions make permanent changes to material properties. Generally it's safer to call spheres3d, as all of the *3d versions of functions make local changes. But there should be no differences in that between Ubuntu and Windows. Can you put together a simple example that does give differences? For example, on Windows this gives 5 different colours: rgl.spheres(1:5, 1:5, 1:5, col=1:5, radius=(1:5)/10) My preferred version would be spheres3d(1:5, 1:5, 1:5, col=1:5, radius=(1:5)/10) Do they behave the same? Duncan Murdoch __ 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] scatter3d: problem with spheres-color
Dear List, I don't get scatter3d to color the sheres according to the '|' argument. library(car) scatter3d(prestige ~ income + education|type, data=Prestige) The spheres on my screen are all colored the same and they are not conditional on Prestige$type. On the other hand: Fit3d and Ellipse3d are colored according to the group argument. rgl_0.92.879 car_2.0-12 R version 2.15.0 i686-pc-linux-gnu (32-bit) I checked this under windows and: here they are colored according to 'Prestige$type'. Hmm? What goes wrong here, any ideas? thanks, René __ 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] how to get inflection point in binomial glm
Dear All, I have a binomial response with one continuous predictor (d) and one factor (g) (8 levels dummy-coded). glm(resp~d*g, data, family=binomial) Y=b0+b1*X1+b2*X2 ... b7*X7 how can I get the inflection point per group, e.g., P(d)=.5 I would be grateful for any help. Thanks in advance, René __ 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] how to get inflection point in binomial glm
Thanks David and Rubén! @David: indeed 15 betas I forgot the interaction terms, thanks for correcting! @Rubén: the re-parameterize would be done within nls()? how to do this practically with including the factor predictor? and yes, we can solve within each group for Y=0 getting 0=b0+b1*X |-b0 -b0=b1*X |/b1 -b0/b1=X but I was hoping there might a more general solution for the case of multiple logistic regression. HTH René Zitat von David Winsemius dwinsem...@comcast.net: On Dec 1, 2011, at 8:24 AM, René Mayer wrote: Dear All, I have a binomial response with one continuous predictor (d) and one factor (g) (8 levels dummy-coded). glm(resp~d*g, data, family=binomial) Y=b0+b1*X1+b2*X2 ... b7*X7 Dear Dr Mayer; I think it might be a bit more complex than that. I think you should get 15 betas rather than 8. Have you done it? how can I get the inflection point per group, e.g., P(d)=.5 Wouldn't that just be at d=1/beta in each group? (Thinking, perhaps naively, in the case of X=X1 that (Pr[y==1])/(1-Pr[y==1])) = 1 = exp( beta *d*(X==X1) ) # all other terms = 0 And taking the log of both sides, and then use middle school math to solve. Oh, wait. Muffed my first try on that for sure. Need to add back both the constant intercept and the baseline d coefficient for the non-b0 levels. (Pr[y==1])/(1-Pr[y==1])) = 1 = exp( beta_0 + beta_d_0*d + beta_n + beta_d_n *d*(X==Xn) ) And just (Pr[y==1])/(1-Pr[y==1])) = 1 = exp( beta_0 + beta_d_0*d ) # for the reference level. This felt like an exam question in my categorical analysis course 25 years ago. (Might have gotten partial credit for my first stab, depending on how forgiving the TA was that night.) I would be grateful for any help. Thanks in advance, René -- David Winsemius, MD West Hartford, CT __ 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] cca with repeated measures
Thanks again Gavin!, this is very clear and enlightening, vegan is an amazing package! It's a pleasure to use it. the best, René -- Dr. René Mayer Email: ma...@psychologie.tu-dresden.de Research Assistant Phone: +49-351-463-34568 Department of Psychology Fax: +49-351-463-33522 Dresden University of Technology Web: http://tu-dresden.de/en Zitat von Gavin Simpson gavin.simp...@ucl.ac.uk: On Fri, 2011-11-18 at 15:17 +0100, René Mayer wrote: Thanks a lot Gavin!, this was what I was looking for. Have I got this right that with no 'cyclic shifts *within* strata' you mean that I cannot define a nesting within animal, e.g., animal/year/season (speaking in regression-terms random-effects for the animal-specific season and year variation). In the permutation framework, you condition the permutation on animal (`strata = animal` in the anova() method), but that alone would say that observations are freely exchangeable within each animal, but not exchangeable between animals. As you have time series data then the observations are not really freely exchangeable within animal either. We can try to maintain the within-animal temporal dependence structure by using cyclic shift permutations: require(permute) shuffleSeries(1:10) [1] 4 5 6 7 8 9 10 1 2 3 t(replicate(10, shuffleSeries(1:10))) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1012345678 9 [2,]12345678910 [3,]3456789 101 2 [4,] 1012345678 9 [5,]789 1012345 6 [6,]89 10123456 7 [7,]9 101234567 8 [8,]3456789 101 2 [9,]456789 1012 3 [10,]56789 10123 4 cyclic shifts basically take a random starting point for the permuted timeseries and then keep sample in the order they appear in the data, bending the ends of the series round together into a circle. The above example: shuffleSeries(1:10) [1] 4 5 6 7 8 9 10 1 2 3 shows a single cyclic shift of the observations 1:10 in an equally spaced time series. Here the random starting point was observation 4 and note that we join the ends of the timeseries so after observation 10, we have observation 1. There are problems with this if there is a trend because you'll get a discontinuity when you join the two ends. As you only have two years of data, and include season as a fixed effect (main term in the model), you could assume that the residuals within animal might be free of temporal dependence. You can check that as I said by fitting the model and looking at the residuals within animal over time. Unless you are prepared to do some coding, the above discussion is moot as vegan doesn't possess the code to use these restricted permutations yet. If you want to do it yourself by hacking the anova.ccabyterm() function, take a look at the vignette that comes with the permute package for ideas on how to specify the permutation design you want and then generate the permutations using shuffle() or shuffleSet() instead of via permuted.index() that is used in vegan. Basically, with: mod - cca(food ~ season*sex, data) anova(mod, by = term, strata = data$animal) I guess you are allowing for a random intercept in animal only. For animal specific season and year effects you'll need to include animal and year in the model with appropriate interactions - we don't, as far as I know, have ways of dealing with more complex dependencies in the residuals via permutations built into vegan other than permuting within animal, or maintaining temporal structure within animal (if you hack the code yourself using functions from permute). HTH G thanks, René -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% 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] cca with repeated measures
Dear all, How can I run a constrained correspondence analysis with the following data: 15 animals were measured repeatedly month-wise (over to 2 years) according to ther diet composition (8 food categories). our data.frame looks like this: food 1 2 ... 8 sex season year animal freq 12 8 ... 1 0 summer 2011 1 freq 0 7 ... 0 1 winter 2011 1 ... freq 0 7 ... 0 1 spring 2011 15 We want to find out if season and sex influences diet composition. My experience with CCA is limited, but in repeated measures ANOVA, e.g. with aov() on has to define the between (animal) error term in order to deal with the pseudoreplication. Do I have to restructure or reshape the data in order to deal with pseudoreplication the data? Or do I have to define an error strata? I suspect I cannot simply run: library(vegan) model=cca(food ~ season*sex+year+animal, data) I would be grateful for any help. Thanks in advance, René __ 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] cca with repeated measures
Thanks a lot Gavin!, this was what I was looking for. Have I got this right that with no 'cyclic shifts *within* strata' you mean that I cannot define a nesting within animal, e.g., animal/year/season (speaking in regression-terms random-effects for the animal-specific season and year variation). thanks, René -- Dr. René Mayer Email: ma...@psychologie.tu-dresden.de Research Assistant Phone: +49-351-463-34568 Department of Psychology Fax: +49-351-463-33522 Dresden University of Technology Web: http://tu-dresden.de/en Zitat von Gavin Simpson gavin.simp...@ucl.ac.uk: On Fri, 2011-11-18 at 10:25 +0100, René Mayer wrote: Dear all, How can I run a constrained correspondence analysis with the following data: 15 animals were measured repeatedly month-wise (over to 2 years) according to ther diet composition (8 food categories). our data.frame looks like this: food 1 2 ... 8 sex season year animal freq 12 8 ... 1 0 summer 2011 1 freq 0 7 ... 0 1 winter 2011 1 ... freq 0 7 ... 0 1 spring 2011 15 We want to find out if season and sex influences diet composition. My experience with CCA is limited, but in repeated measures ANOVA, e.g. with aov() on has to define the between (animal) error term in order to deal with the pseudoreplication. Do I have to restructure or reshape the data in order to deal with pseudoreplication the data? Or do I have to define an error strata? I suspect I cannot simply run: library(vegan) model=cca(food ~ season*sex+year+animal, data) You could do that although the analysis would be i) focussed on those particular animals in those years, and ii) you could only test the terms season, sex and season:sex in a sequential manner (i.e. dependent upon how the terms enter the model), so season, then sex after season is in the model, then their interaction after both main terms are included in the model. ii) is done by adding `by = terms` to the call to the `anova method for cca objects; examples are in `?anova.cca` That corresponds to a fixed effects formulation of the ANOVA (assuming I have my terminology right). The alternative is to adjust the permutation scheme used to reflect the clustering in your data. In that case, using `strata = animal` would be OK. Ideally one would want to control for temporal dependence so you would want cyclic shifts *within* `strata = animal` but vegan can not yet do this sort of permutation. It is coming; the actual code to generate those permutations is available in the permute package (upon which vegan depends), but as yet we have not hooked this into the vegan ordination functions (it is on the TODO list). That said, season should be accounting for much of the temporal dependence, so I think you might get away with just specifying strata (the permutation test itself is permuting residuals from the model so if the model is capturing the seasonal variation then the simple permuting within animal should be OK, but you can check by extracting the residuals using the resid() method and plotting them out by time and animal). HTH G I would be grateful for any help. Thanks in advance, René __ 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.
[R] lattice-dotplot: resize axis
dear all, I want to make a dotplot with ratings from Items in 6 ItemsGroups. I reordered the items by rating within each group. I plotted the items by rating conditional on ItemGroup. The ordering works as I wanted but my y-aches labels (items) within each ItemGroup are now unequally spaced, e.g., in some panels there is a gap between one lower rated item and the next higher, to give a picture items=a,e,f,g ItemGroup=n - g| . f| . e| . | | | a| . - How can I correct this? What have I overlooked? # code i've used (from latticeExtra/utilities/resize panels) library(latticeExtra) mean.ratings$item.name - with(mean.ratings, reorder(reorder(item, rating), as.numeric(ItemGroup))) dpratings - dotplot(item.name ~ rating | reorder(ItemGroup, rating), data = mean.ratings, layout = c(1, 6), xlim=c(1,6), aspect = .1, scales = list(y = list(relation = free, cex=.5))) ## approximate resizePanels(dpratings, h = with(mean.ratings, table(reorder(ItemGroup, rating thanks, René __ 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] prcomp: results with reversed sign in output?
Dear All, when I'm running a PCA with prcomp(USArrests, scale = TRUE) I get the right principal components, but with the wrong sign infront Rotation: PC1 PC2 PC3 PC4 Murder 0.5358995 -0.4181809 0.3412327 0.64922780 Assault 0.5831836 -0.1879856 0.2681484 -0.74340748 UrbanPop 0.2781909 0.8728062 0.3780158 0.13387773 Rape 0.5434321 0.1673186 -0.819 0.08902432 instead of PC1 PC2 PC3 PC4 Murder -0.5358995 0.4181809 -0.3412327 0.64922780 Assault -0.5831836 0.1879856 -0.2681484 -0.74340748 UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773 Rape -0.5434321 -0.1673186 0.819 0.08902432 what is happening here? any ideas? thanks, René __ 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] prcomp: results with reversed sign in output?
thanks for pointing out Paul, but the thing which is annoying me in the first place IS this direction reversal. this makes no sense for me why could this be? Zitat von Paul Hiemstra paul.hiems...@knmi.nl: Hi, If all the signs are switched the PC's are still the same. The principal vectors are along the same axis, only in a different direction. So there is no problem :). hope this helps, Paul On 09/09/2011 09:01 AM, René Mayer wrote: Dear All, when I'm running a PCA with prcomp(USArrests, scale = TRUE) I get the right principal components, but with the wrong sign infront Rotation: PC1 PC2 PC3 PC4 Murder 0.5358995 -0.4181809 0.3412327 0.64922780 Assault 0.5831836 -0.1879856 0.2681484 -0.74340748 UrbanPop 0.2781909 0.8728062 0.3780158 0.13387773 Rape 0.5434321 0.1673186 -0.819 0.08902432 instead of PC1 PC2 PC3 PC4 Murder -0.5358995 0.4181809 -0.3412327 0.64922780 Assault -0.5831836 0.1879856 -0.2681484 -0.74340748 UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773 Rape -0.5434321 -0.1673186 0.819 0.08902432 what is happening here? any ideas? thanks, René __ 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. -- Paul Hiemstra, Ph.D. Global Climate Division Royal Netherlands Meteorological Institute (KNMI) Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39 P.O. Box 201 | 3730 AE | De Bilt tel: +31 30 2206 494 http://intamap.geo.uu.nl/~paul http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770 __ 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] prcomp: results with reversed sign in output?
thanks for explaining Duncan and Ted, Indeed, I did compare my results from a textbook and noticed that I consitenly get flipped signs and biplots. regards René Zitat von ted.hard...@wlandres.net: The point is that a principal component vector is a solution, say V, of a matrix equation A%*%V = L*V where A is the matrix and L is a scalar.. Since this equation can be written A%*%(-V) = L*(-V), the result is indeterminate with respect to its sign. If V is a solution, so is (-V), and vice versa. It is not a case of direction reversal, since neither L nor (-L) has a primary role -- they are equivalent, and you can adopt either one. Just make it clear which one you adopt -- or someone else might think that they disagree! You originally wrote I get the right principal components, but with the wrong sign in front. You did not get the wrong sign -- both are correct! It may be that you are comparing your result from R with the result from some other software (or from a textbook, or whatever) which produced the equivalent result but with the opposite sign. Again, both are correct. If, for some reason, you do not like the sign of the result you get, then change its sign. Hoping this helps, Ted. On 09-Sep-11 09:42:49, René Mayer wrote: thanks for pointing out Paul, but the thing which is annoying me in the first place IS this direction reversal. this makes no sense for me why could this be? Zitat von Paul Hiemstra paul.hiems...@knmi.nl: Hi, If all the signs are switched the PC's are still the same. The principal vectors are along the same axis, only in a different direction. So there is no problem :). hope this helps, Paul On 09/09/2011 09:01 AM, René Mayer wrote: Dear All, when I'm running a PCA with prcomp(USArrests, scale = TRUE) I get the right principal components, but with the wrong sign infront Rotation: PC1 PC2 PC3 PC4 Murder 0.5358995 -0.4181809 0.3412327 0.64922780 Assault 0.5831836 -0.1879856 0.2681484 -0.74340748 UrbanPop 0.2781909 0.8728062 0.3780158 0.13387773 Rape 0.5434321 0.1673186 -0.819 0.08902432 instead of PC1 PC2 PC3 PC4 Murder -0.5358995 0.4181809 -0.3412327 0.64922780 Assault -0.5831836 0.1879856 -0.2681484 -0.74340748 UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773 Rape -0.5434321 -0.1673186 0.819 0.08902432 what is happening here? any ideas? thanks, René -- Paul Hiemstra, Ph.D. Global Climate Division Royal Netherlands Meteorological Institute (KNMI) Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39 P.O. Box 201 | 3730 AE | De Bilt tel: +31 30 2206 494 http://intamap.geo.uu.nl/~paul http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770 E-Mail: (Ted Harding) ted.hard...@wlandres.net Fax-to-email: +44 (0)870 094 0861 Date: 09-Sep-11 Time: 11:04:16 -- XFMail -- __ 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] rgl how to plot a cylinder like arrow3d?
Dear List,
I'm trying to draw vector in XYZ with rgl under use of a cylinder3d.
Therefore I scale and rotate a basis-cylinder).
However, somehow the rotation is wrong as
verified by overplotting arrow3d().
Where is my mistake?
library(heplots)
library(rgl)
# ... 2 vectors
data=data.frame(row.names=c('X','Y','Z'), x1=c(2,1,5),y=c(4,3,2))
vector3D=function(start=c(0,0,0), end, mycol='green'){
# ... cylinder as basis-vector
c=cylinder3d(rbind( c(start), # start
c(0,0,1)), # end
radius = 0.2,
e1=cbind(0, 0, 1),
e2=cbind(1, 0, 0),
sides=10
)
# ... rotate cylinder horizontally and scale it
len=sqrt(sum(abs(end)*abs(end)))
c=scale3d(c,0.1,0.1,len)
c=rotate3d(c,-(pi/2),0,1,0)
# ... rotation angles
theta = atan2(end[2],end[1]) # angle in yx-plane for Z-achses rotation
phi = atan2(end[3],sqrt(end[1]^2+end[1]^2)) # angle in zx for Y-achses
rotation
# ... rotation
c=rotate3d(c,phi,0,1,0) # Y-achses rotation
c=rotate3d(c, -theta,0,0,1) # Z-achses rotation
shade3d(addNormals(c), col=mycol)
}
open3d()
vector3D(start=c(0,0,0),end=data[,1], 'red');
vector3D(start=c(0,0,0),end=data[,2],'green')
arrow3d(c(0,0,0),data[,1],color='red');arrow3d(c(0,0,0),data[,2],
color='green')
axes3d(c('x','y','z'));title3d('main','sub','X','Y','Z');box3d()
thanks,
René
__
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] rgl how to plot a cylinder like arrow3d?
Thanks Duncan!
yes! this works;
I paste the code in case someone wants to draw
3d-vectors with cylinders and cones.
René
vector3D=function(start=c(0,0,0), end, mycol='green', cone.length=0.1, ... ){
# ... cylinder as basis-vector
c=cylinder3d(rbind( c(start), # start
c(0,0,1)), # end
radius = 0.4,
e1=cbind(0, 0, 1),
e2=cbind(1, 0, 0),
sides=10
)
# ... rotate cylinder horizontally and scale it
len=sqrt(sum(abs(end)*abs(end)))
c=scale3d(c,0.1,0.1,len-cone.length)
# ... Duncan Murdoch's code snipe
rotation -
rgl:::GramSchmidt(end-start+c(1,0,0),end-start+c(0,1,0),end-start,
order=c(3,1,2))
c - rotate3d(c, matrix=rotation)
shade3d(addNormals(c), col=mycol)
# ... cone3d from rgl-demos
q1 - cone3d(qmesh=T,trans=diag(4)) # height=1,radius=1, base at (0,0,0)
q1=translate3d(scale3d(q1,0.15,0.15,1.4),0,0,len-cone.length)
q1=rotate3d(q1, matrix=rotation) # hinlegen ... negative rotation ==
zeigerrrichtung
shade3d(q1,col=mycol) # verschachtelt aufrufen
}
open3d()
vector3D(start=c(0,0,0),end=data[,1], 'red');
vector3D(start=c(0,0,0),end=data[,2],'green')
arrow3d(c(0,0,0),data[,1],color='red');arrow3d(c(0,0,0),data[,2],
color='green')
axes3d(c('x','y','z'));title3d('main','sub','X','Y','Z');box3d()
On 11-08-09 5:22 AM, René Mayer wrote:
Dear List,
I'm trying to draw vector in XYZ with rgl under use of a cylinder3d.
Therefore I scale and rotate a basis-cylinder).
However, somehow the rotation is wrong as
verified by overplotting arrow3d().
Where is my mistake?
I would guess it is in assuming that theta can be computed on the
original vector, not on the rotated one: rotating about the Y axis
changes end[1].
I usually use a Gram-Schmidt type calculation to do this sort of
thing. The undocumented rgl:::GramSchmidt function does this for
3x3 matrices, orthogonalizing by row. So you could replace all 3 of
your rotations (including the first one) by this:
rotation -
rgl:::GramSchmidt(end-start+c(1,0,0),end-start+c(0,1,0),end-start,
order=c(3,1,2))
c - rotate3d(c, matrix=rotation)
(I notice you also forgot to translate the cylinders to start; this
doesn't matter in your example, but would if you had a non-zero
origin.)
One problem with this kind of display is that it assumes the
coordinates are all of equal range. If you try plotting one of
these vectors when the range of x, y and z are drastically
different, the cylinders will look really strange.
Duncan Murdoch
library(heplots)
library(rgl)
# ... 2 vectors
data=data.frame(row.names=c('X','Y','Z'), x1=c(2,1,5),y=c(4,3,2))
vector3D=function(start=c(0,0,0), end, mycol='green'){
# ... cylinder as basis-vector
c=cylinder3d(rbind( c(start), # start
c(0,0,1)), # end
radius = 0.2,
e1=cbind(0, 0, 1),
e2=cbind(1, 0, 0),
sides=10
)
# ... rotate cylinder horizontally and scale it
len=sqrt(sum(abs(end)*abs(end)))
c=scale3d(c,0.1,0.1,len)
c=rotate3d(c,-(pi/2),0,1,0)
# ... rotation angles
theta = atan2(end[2],end[1]) # angle in yx-plane for Z-achses rotation
phi = atan2(end[3],sqrt(end[1]^2+end[1]^2)) # angle in zx for Y-achses
rotation
# ... rotation
c=rotate3d(c,phi,0,1,0) # Y-achses rotation
c=rotate3d(c, -theta,0,0,1) # Z-achses rotation
shade3d(addNormals(c), col=mycol)
}
open3d()
vector3D(start=c(0,0,0),end=data[,1], 'red');
vector3D(start=c(0,0,0),end=data[,2],'green')
arrow3d(c(0,0,0),data[,1],color='red');arrow3d(c(0,0,0),data[,2],
color='green')
axes3d(c('x','y','z'));title3d('main','sub','X','Y','Z');box3d()
thanks,
René
__
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] how to merge within range?
Hello, how can one merge two data frames when in the second data frame one column defines the start values and another defines the end value of the to be merged range. data.frame.1 time ... 13 24 35 46 55 ... data.frame.2 start end 24 37 ?h? ? ... should result in this 13 NA 24 ?h? 35 ?h? 46 NA 55 ? thanks, René __ 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] how to merge within range?
thanks David and Ian, let me make a better example as the first one was flawed df.1=data.frame(round((1:10)*100+rnorm(10)), value=NA) names(df.1) = c(time, value) df.1 time value 1 101NA 2 199NA 3 301NA 4 401NA 5 501NA 6 601NA 7 700NA 8 800NA 9 900NA 10 1000NA # from and to define ranges within time, # note that from and to may not match the numbers given in time df.2=data.frame(from=c(99,500,799),to=c(303,702,950), value=c(1,3,5)) df.2 from to value 1 99 303 1 2 500 702 3 3 799 950 5 what I want is: time value 1 1011 2 1991 3 3011 4 401NA 5 5013 6 6013 7 7003 8 8005 9 9005 10 1000NA @David I don't know what you mean by 2 merges, René Zitat von David Winsemius dwinsem...@comcast.net: On May 14, 2011, at 9:16 AM, Ian Gow wrote: If I assume that the third column in data.frame.2 is named val then in SQL terms it _seems_ you want SELECT a.time, b.val FROM data.frame.1 AS a LEFT JOIN data.frame.2 AS b ON a.time BETWEEN b.start AND b.end; Not sure how to do that elegantly using R subsetting/merge, Huh? It's just two merge()'s (... once you fix the error in the example.) -- David but you might try a package that allows you to use SQL, such as sqldf. On 5/14/11 8:03 AM, David Winsemius dwinsem...@comcast.net wrote: On May 14, 2011, at 8:12 AM, René Mayer wrote: Hello, how can one merge And what happened when you typed: ?merge two data frames when in the second data frame one column defines the start values and another defines the end value of the to be merged range. data.frame.1 time ... 13 24 35 46 55 ... data.frame.2 start end 24 37 ?h? ? ... should result in this 13 NA 24 ?h? 35 ?h? 46 NA 55 ? And _why_ would that be? thanks, René __ 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. David Winsemius, MD West Hartford, CT __ 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. David Winsemius, MD West Hartford, CT __ 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] how to merge within range?
sqldf is impressive - compiled it now; the trick with findInterval is nice, too. thanks guys!! Zitat von David Winsemius dwinsem...@comcast.net: On May 14, 2011, at 2:27 PM, William Dunlap wrote: You could use findInterval() along with a trick with c(rbind(...)): i - findInterval(x=df.1$time, vec=c(rbind(df.2$from, df.2$to))) i [1] 1 1 1 2 3 3 3 5 5 6 That's nice. I was working on a slightly different trick findInterval( df.1[,1],t(df.2[,1:2])) [1] 1 1 1 2 3 3 3 5 5 6 I was then trying to get the right indices with (.)'%%' 2 and (.) '%/%' 2 The even-valued outputs would map to NA's, the odds to value[(i+1)/2], but you can use the c(rbind(...)) trick again: c(rbind(df.2$value, NA))[i] [1] 1 1 1 NA 3 3 3 5 5 NA I'd like to understand that. Maybe, maybe... ah, got it. At first I didn't realize those were the final answers since they looked like indices. My t(.) trick doesn't generalize as well. My earlier suggestion tht two merges woul do it was based on my erroneous interpretation of the example, since I thought the task was to match on the end points of the intervals. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of René Mayer Sent: Saturday, May 14, 2011 11:06 AM To: David Winsemius Cc: r-help@r-project.org Subject: Re: [R] how to merge within range? thanks David and Ian, let me make a better example as the first one was flawed df.1=data.frame(round((1:10)*100+rnorm(10)), value=NA) names(df.1) = c(time, value) df.1 time value 1 101NA 2 199NA 3 301NA 4 401NA 5 501NA 6 601NA 7 700NA 8 800NA 9 900NA 10 1000NA # from and to define ranges within time, # note that from and to may not match the numbers given in time df.2=data.frame(from=c(99,500,799),to=c(303,702,950), value=c(1,3,5)) df.2 from to value 1 99 303 1 2 500 702 3 3 799 950 5 what I want is: time value 1 1011 2 1991 3 3011 4 401NA 5 5013 6 6013 7 7003 8 8005 9 9005 10 1000NA @David I don't know what you mean by 2 merges, René Zitat von David Winsemius dwinsem...@comcast.net: On May 14, 2011, at 9:16 AM, Ian Gow wrote: If I assume that the third column in data.frame.2 is named val then in SQL terms it _seems_ you want SELECT a.time, b.val FROM data.frame.1 AS a LEFT JOIN data.frame.2 AS b ON a.time BETWEEN b.start AND b.end; Not sure how to do that elegantly using R subsetting/merge, Huh? It's just two merge()'s (... once you fix the error in the example.) -- David but you might try a package that allows you to use SQL, such as sqldf. On 5/14/11 8:03 AM, David Winsemius dwinsem...@comcast.net wrote: On May 14, 2011, at 8:12 AM, René Mayer wrote: Hello, how can one merge And what happened when you typed: ?merge two data frames when in the second data frame one column defines the start values and another defines the end value of the to be merged range. data.frame.1 time ... 13 24 35 46 55 ... data.frame.2 start end 24 37 ?h? ? ... should result in this 13 NA 24 ?h? 35 ?h? 46 NA 55 ? And _why_ would that be? thanks, René __ 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. David Winsemius, MD West Hartford, CT __ 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. David Winsemius, MD West Hartford, CT __ 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. David Winsemius, MD West Hartford, CT __ 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] lmer - how to test correlations?
Dear mixed modelers, If I have interactions between a categorical covariate and a random-effects grouping factor. How can I test formally the correlations? For example data(Machines, package = MEMSS) fm2aM - lmer(score ~ Machine + (0 + Machine|Worker), Machines) Random effects: Groups Name Variance Std.Dev. Corr Worker MachineA 13.81574 3.71695 MachineB 61.94506 7.87052 0.805 MachineC 16.00492 4.00061 0.625 0.772 How do I formally test the significance for every single Corr(A,B); Corr(A,C); Corr(B,C)? thanks for any help, René __ 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] How to code repeated measures interaction in lmer?
Dear Everybody: I want to test an interaction of two repeated measures in lmer() I've response times (Y) from N subjects with two within-subject-factors. aov(Y ~ A*B + Error(subjects/A*B), ...) shows an insigificant interaction p(A:B) = .35 now lmer() lm1 = lmer(Y ~ 1+A+B+(1+A*B|subjects), ...) lm2 = lmer(Y ~ 1+A*B+(1+A*B|subjects), ...) p(A:B) = anova(lm1,lm2) is this right? and when do I have to write lm1 = lmer(Y ~ 1+A+B+(1|subjects)+(1|subjects:A:B), ...) lm2 = lmer(Y ~ 1+A*B+(1|subjects)+(1|subjects:A:B), ...) p(A:B) = anova(lm1,lm2) ? many thanks in advance, René __ 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] how to analyze repeated measures count data?
Dear R community, I've data-set with reaction times and count data (answers - yes, no) of N subjects under conditions A, B. For the analysis reaction time I used aov. fit.rt = aov(rt ~ A * B + Error(subjects/(A*B)), data = m ) But how do I analyze the frequencies correctly? example fable of frequencies from one subject: , , = A1 B1 B2 B3 yes 31 3619 no22 2710 , , = A2 B1 B2B3 yes 22 2710 no31 3619 Is a generalized linear model the right method? How do I specify the same model for the count data (frequencies) in glm? is this right: glm(count~A*B*answer+(1|subject),family=poisson)? Regards, René __ 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] interpolation
Dear R-users, I have a complex line by xy-values (ordered by z). And I would like to get interpolated y-values on the positions of x = 0:600. How do I get the correct points? x=c(790,790,790,790,790,786,783,778,778,766,763,761,761,761,715,628,521,350,160,134,134,129,108,101,93,111,161,249,288,243,139,45,7) y=c(606,606,606,606,606,612,617,627,627,640,641,641,641,641,689,772,877,1048,1240,1272,1272,1258,1242,1239,1239,1214,1122,959,770,479,273,133,45) z=c(0,29,58,87,116,145,174,203,232,261,290,319,348,377,406,435,464,493,522,551,580,609,638,667,696,725,754,783,812,841,870,899,928) plot(y,x,type=b) # this fails ? lines(approx(y,x),col=blue) # with xout = c(0:600) thanks in advance, René -- Dr. rer. nat. Dipl.-Psych. René Mayer Dresden University of Technology Department of Psychology Zellescher Weg 17 D-01062 Dresden Tel.: +49-351-4633-4568 Email: ma...@psychologie.tu-dresden.de __ 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] interpolation
My problem is that x values increas with y until some point then the pattern reverses. The whole line is a kind of U-shape with a right-buttom to middel-top diagonal at the end of it (a look at the plot makes it clearer). The interpolation (approx, spline) makes a zick-zack aut of it. What I need is to interpolate some points and to preserve the shape. thanks in andvance René Zitat von David Winsemius dwinsem...@comcast.net: On Jan 11, 2010, at 7:44 AM, René Mayer wrote: Dear R-users, I have a complex line by xy-values (ordered by z). And I would like to get interpolated y-values on the positions of x = 0:600. How do I get the correct points? x=c(790,790,790,790,790,786,783,778,778,766,763,761,761,761,715,628,521,350,160,134,134,129,108,101,93,111,161,249,288,243,139,45,7) y=c(606,606,606,606,606,612,617,627,627,640,641,641,641,641,689,772,877,1048,1240,1272,1272,1258,1242,1239,1239,1214,1122,959,770,479,273,133,45) z=c(0,29,58,87,116,145,174,203,232,261,290,319,348,377,406,435,464,493,522,551,580,609,638,667,696,725,754,783,812,841,870,899,928) plot(y,x,type=b) That would plot x as a function of y (the reverse of the usual convention. Is that what you want. If so, then the statement that these are ordered by z is misleading. Ordering by z would suggest that each series is a function of z. What sort of interpolation do you want and in how many dimensions? # this fails ? lines(approx(y,x),col=blue) # with xout = c(0:600) thanks in advance, René -- Dr. rer. nat. Dipl.-Psych. René Mayer Dresden University of Technology Department of Psychology Zellescher Weg 17 D-01062 Dresden Tel.: +49-351-4633-4568 Email: ma...@psychologie.tu-dresden.de __ 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. David Winsemius, MD Heritage Laboratories West Hartford, CT __ 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] SPSS repeated interaction contrast in R
dear all, i'm trying to reproduce an spss-anova in R. It is an 2x3x3 repeated measures desingn with repeated contrasts. In R i've coded a contrast matrix for all factors and made a split in the aov summary - but I can't get the repeated interaction contrasts. The output from SPSS looks like this: TaskSw * CongNow * CongBefore: SS df Mean SquareF Sig. 1 vs. 2 1 vs. 2 1 vs. 213944,50 1 13944,50 0.37 0.56 2 vs. 34278,12 1 4278,12 0.31 0.59 2 vs. 3 1 vs. 27140,12 1 7140,12 0.17 0.68 2 vs. 353301,131 1 53301,13 1.38 0.27 The output from R looks like this: Error: Subj:TaskSw:CongNow:CongBefore Df Sum Sq Mean Sq F value Pr(F) 4 188234706 1.8804 0.1417 rep vs. se.con vs. inc.con vs. inc 1 inc vs. neu.1 con vs. inc.inc vs. neu 1 inc vs. neu.1 Residuals 28 700692502 I've pasted the R and SPSS code I uesd. thanks in advance. Rene R code: # read in data TaskSwitch - factor(rep(c(0:1), 72), levels=c(0,1), labels=c(Repetition,Switch)) CongruenceNow - factor(rep(c(0,0,1,1,2,2), 24), levels=c(0,1,2), labels=c(Congruent,Incongruent,Neutral)) CongruenceBefore - factor(rep(c(0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2), 8), levels=c(0,1,2), labels=c(Congruent,Incongruent,Neutral)) # Subjects Subjects - factor(rep(c(1:8), each = 18)) # RT = Response Time RT - c(648,682,798,765,707,677,676,698,712,770,702,683,635,716,748,703,737, 657,1047,1140,992,1088,925,967,1022,1105,884,1103,930,1046,1034,1032, 948,956,987,993,753,782,837,923,820,846,835,884,778,873,774,760,848, 824,788,864,969,833,662,683,848,903,731,755,778,748,792,913,707,721, 736,766,771,898,874,801,762,824,956,865,864,908,828,877,804,1096,864, 878,799,833,883,1005,863,824,691,666,698,791,712,778,731,743,749,781, 717,769,732,730,703,771,743,709,470,489,546,574,541,523,536,556,532, 510,543,557,531,554,532,573,528,514,842,785,851,932,888,920,856,840, 1010,877,784,901,845,923,798,907,903,758) m - data.frame(RT, TaskSwitch,CongruenceNow,CongruenceBefore,Subjects) # defining repeated contrasts for all factors conCN = matrix(c(1,0, -1,1, 0,-1), nrow = 3, ncol = 2, byrow = TRUE) contrasts(m$CongruenceNow) - conCN conTS = matrix(c(1, -1), nrow = 2, ncol = 1, byrow = TRUE) contrasts(m$TaskSwitch) - conTS conCB = matrix(c(1,0, -1,1, 0,-1), nrow = 3, ncol = 2, byrow = TRUE) contrasts(m$CongruenceBefore) - conCB # ANOVA fit - aov(RT ~ TaskSwitch*CongruenceNow*CongruenceBefore + Error(Subjects/(TaskSwitch*CongruenceNow*CongruenceBefore)), data = m) summary(fit) summary(fit, split = list(CongruenceNow = list(con vs. inc = 1, inc vs. neu = 2), CongruenceBefore= list(con vs. inc = 1, inc vs. neu = 2), TaskSwitch= list(rep vs. se = 1))) # export to spss require(foreign) codefile-tempfile() write.foreign(m, paste(C:/m.sav,sep = ''), codefile,package=SPSS) SPSS: ## BEGIN: SPSS code ### # GET DATA # /TYPE=TXT /FILE='C:\m.sav' /DELCASE=LINE /DELIMITERS=, /ARRANGEMENT=DELIMITED # /FIRSTCASE=1/IMPORTCASE=ALL # /VARIABLES= # RT F4.0 TaskSwitch F1.0 CongruenceNow F1.0 CongruenceBefore F1.0 Subjects F1.0. # CACHE. # EXECUTE. # DATASET NAME DataSet1 WINDOW=FRONT. # # SORT CASES BY Subjects TaskSwitch CongruenceNow CongruenceBefore . # CASESTOVARS # /ID= Subjects # /INDEX=TaskSwitch CongruenceNow CongruenceBefore # /GROUPBY=VARIABLE. # # GLM RT.1.1.1 RT.1.1.2 RT.1.1.3 RT.1.2.1 RT.1.2.2 RT.1.2.3 RT.1.3.1 RT.1.3.2 RT.1.3.3 RT.2.1.1 # RT.2.1.2 RT.2.1.3 RT.2.2.1 RT.2.2.2 RT.2.2.3 RT.2.3.1 RT.2.3.2 RT.2.3.3 # /WSFACTOR=TaskSwitch 2 Repeated CongruenceNow 3 Repeated CongruenceBefore 3 Repeated # /METHOD=SSTYPE(3) # /PRINT=DESCRIPTIVE OPOWER TEST(SSCP) # /CRITERIA=ALPHA(.05) # /WSDESIGN=TaskSwitch CongruenceNow CongruenceBefore TaskSwitch*CongruenceNow # TaskSwitch*CongruenceBefore CongruenceNow*CongruenceBefore # TaskSwitch*CongruenceNow*CongruenceBefore. ## END: SPSS code ### -- Dr. rer. nat. Dipl.-Psych. René Mayer Dresden University of Technology Department of Psychology Zellescher Weg 17 D-01062 Dresden Tel.: +49-351
