Re: [R] getting all contrasts from glm
Jim, In the glm object I can find the contrasts of the main treats vs the first i.e. 2v1, 3v1 and 4v1 ... however I would like to get the complete set including 3v2, 4v2, and 4v3 ... along with the Std. Errors of all contrasts. Your best all round approach would be to use the multcomp package. In particular, look at ?glht ?summary.glht ?contrMat Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/getting-all-contrasts-from-glm-tp3007027p3007421.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 table and lmer
Hi Jim, The decomposition of the sum of squares should be the same regardless of whether block is treated as random of fixed. Should it? By whose reckoning? The models you are comparing are different. Simple consideration of the terms listed in the (standard) ANOVA output shows that this is so, so how could the sum-of-squares be the same? I noticed that other people have asked similar questions in the past, but I haven't seen a satisfactory explanation. Maybe, but it has been answered (by me, and surely by others). However, canonical would be Venables and Ripley's MASS (: 283--286). The models you need to compare are the following: ## Aov.mod - aov(Y ~ V * N + Error(B/V/N), data = oats) Lme.mod - lme(Y ~ V * N, random = ~1 | B/V/N, data = oats) Lmer.mod - lmer(Y~ V * N +(1|B)+(1|B:V)+(1|B:N), data = oats) summary(Aov.mod) anova(Lme.mod) anova(Lmer.mod) HTH, Mark Difford. -- View this message in context: http://r.789695.n4.nabble.com/ANOVA-table-and-lmer-tp3027546p3027662.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] changing the limits of a secondary y-axis in a barplot
Hi Anna, How can I change the barplot so that the left hand axis scales from 0 to 15 and the right hand axis from 0 to 5? Try this: par(mfrow=c(1,1), mai=c(1.0,1.0,1.0,1.0)) Plot1-barplot(rbind(Y1,Y2), beside=T, axes=T, names.arg=c(a,b), ylim=c(0,15), xlim=c(1,9), space=c(0,1), col=c(darkgray,white), yaxt=n) Plot2-barplot(Y3, add=T, beside=T, names.arg=c, col=c(darkgray,white), ylim=c(0,5), space=c(0,7), width=1, yaxt=n) axis(side=2, at=seq(0,15,3), labels=seq(0,15,3)) axis(side=4, at=seq(0,15,3), labels=seq(0,5,1)) Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/changing-the-limits-of-a-secondary-y-axis-in-a-barplot-tp3046117p3046283.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] Wilcoxon Rank Sum in R with a multiple testing correction
Hi Selthy, I'd like to use a Wilcoxon Rank Sum test to compare two populations of values. Further, I'd like to do this simultaneously for 114 sets of values. Well, you read your data set into R using: ## ?read.table ?read.csv There are other ways to bring in data. Save the import to a workspace object at the same time: myDat - read.csv (...) Do the Wilcoxon Rank Sum tests using the implementation of your choice (there are several): ## See the examples at foot of help page. Lacking data we will make some. ?wilcox.test pv1 - wilcox.test(rnorm(10), rnorm(10, 2), conf.int = TRUE)$p.value pv2 - wilcox.test(rnorm(10), rnorm(10, 2), conf.int = TRUE)$p.value pv3 - wilcox.test(rnorm(10), rnorm(10, 2), conf.int = TRUE)$p.value Eventually you will discover more elegant ways of assembling a vector (or some other type of storage object). Finally, you feed your p-values to: ## ?p.adjust pAdj - p.adjust (c(pv1, pv2, pv3), method = c(BH)) ## ?round ?sprintf cbind.data.frame (Uncorrected = c(pv1, pv2, pv3), BH_Corrected = pAdj) Eventually you will discover how to turn all of this into an elegant function. I really do hope that this is not a school assignment. If so Well, you still need to do some work to get this going. Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/Wilcoxon-Rank-Sum-in-R-with-a-multiple-testing-correction-tp3056557p3056878.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] Factor analysis and cfa with asymptotically distributed data
Jane, Does someone know how to do fa and cfa with strong skewed data? Your best option might be to use a robustly estimated covariance matrix as input (see packages robust/robustbase). Or you could turn to packages FAiR or lavaan (maybe also OpenMx). Or you could try soft modelling via package plspm. Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/Factor-analysis-and-cfa-with-asymptotically-distributed-data-tp3056387p3056944.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] Compute polychoric matrix for more than 2 variables
Hi Raquel, routine in R to compute polychoric matrix to more than 2 categorical variables? I tried polycor package, but it seems to be suited only to 2-dimensional problems. Bur surely ?hetcor (in package polycor) does it. Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/Compute-polychoric-matrix-for-more-than-2-variables-tp3064941p3065027.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] reading lmer table
Hi Nicola, In few word: does this row indicate a global effect of the predictor 'cat' or a more specific passage? It indicates a more specific passage. Use anova(m7) for global/omnibus. Check this for yourself by fitting the model with different contrasts. The default contrasts in R are treatment contrasts. ## m7 - lmer(log.second ~ Cond + cat + (1|subjID) + (1|Code), data = march.f, contrasts=list(Cond=contr.treatment, cat=contr.treatment)) m7s - lmer(log.second ~ Cond + cat + (1|subjID) + (1|Code), data = march.f, contrasts=list(Cond=contr.sum, cat=contr.sum)) summary(m7) summary(m7s) anova(m7) anova(m7s) Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/reading-lmer-table-tp2329521p2330809.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] Coinertia randtest
Hi Petar, I dunno why, but I cannot make randtes[t].coinertia() from ade4 package working. I have two nice distance matrices (Euclidean): Could anyone help with this? Yes (sort of). The test has not yet been implemented for dudi.pco, as the message at the end of your listing tells you. The result is: Error in randtest.coinertia(coi, nrepet = 1000) : Not yet available So far it has only been implemented for some types of dudi.pca, and for dudi.coa, dudi.fca, and dudi.acm. You might be lucky and find code to do what you want on the ade4 list. http://pbil.univ-lyon1.fr/ADE-4/home.php?lang=eng http://listes.univ-lyon1.fr/wws/info/adelist Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/Coinertia-randtest-tp2335696p2335748.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] path analysis
Guy, For a partial least squares approach look at packages plspm and pathmox. Also look at sem.additions. Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/path-analysis-tp2528558p2530207.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] Creating publication-quality plots for use in Microsoft Word
I'd prefer to stick with JPEG, TIFF, PNG, or the like. I'm not sure EPS would fly. Preferring to stick with bitmap formats (like JPEG, TIFF, PNG) is likely to give you the jagged lines and other distortions you profess to want to avoid. EPS (encapsulated postscript, which handles vector+bitmap) is one of the graphic file formats preferred by most quality journals. Surprisingly, not too many people seem to be aware of the fact that PDF really is a crippled form of postscript. Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/Creating-publication-quality-plots-for-use-in-Microsoft-Word-tp2540676p2540858.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] PCA or CA
Hi Paul, I have a data set for which PCA based between group analysis (BGA) gives significant results but CA-BGA does not. I am having difficulty finding a reliable method for deciding which ordination technique is most appropriate. Reliability really comes down to you thinking about and properly defining what _information_ you want to extract from your data set, which we know nothing about. PCA and CA are fundamentally different. The classical use of CA lies in the analysis of count-data (contingency tables), for which it remains a brilliant method. It is also widely applied to analyzing normal n x p matrices of the type usually analyzed by PCA. A doubly-centered PCA would get you close to a CA of the normal n x p matrix (i.e. of an analysis of the same matrix). This is a biggish area, so grab your specs, and perhaps start with Jolliffe (PCA) and Benzecri/Greenacre (CA). Regards, Mark. Paul Dennis-3 wrote: Dear all I have a data set for which PCA based between group analysis (BGA) gives significant results but CA-BGA does not. I am having difficulty finding a reliable method for deciding which ordination technique is most appropriate. I have been told to do a 1 table CA and if the 1st axis is2 units go for CA if not then PCA. Another approach is that described in the Canoco manual - perform DCA and then look at the length of the axes. I used decorana in vegan and it gives axis lengths. I assume that these are measured in SD units. Anyway the manual say if the axis length is 3 go for PCA,4 use CA and if intermediate use either. Are either of these approaches good/valid/recommended or is there a better method? Thanks Paul _ Get the best of MSN on your mobile http://clk.atdmt.com/UKM/go/147991039/direct/01/ __ 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/PCA-or-CA-tp25668667p25670451.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] trouble with html() in Hmisc
Hi Liviu,
tmp - latex(.object, cdec=c(2,2), title=)
class(tmp)
[1] latex
html(tmp)
/tmp/RtmprfPwzw/file7e72f7a7.tex:9: Warning: Command not found:
\tabularnewline
Giving up command: \...@hevea@amper
/tmp/RtmprfPwzw/file7e72f7a7.tex:11: Error while reading LaTeX:
This array/tabular column has no specification
This has nothing to do with Hmisc or hevea. In the header/preample of all
LyX files you will, for instance, find the following:
## -
%% LyX 1.6.2 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{article}
.
.
.
%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}
Regards, Mark.
Liviu Andronic wrote:
Dear all
On my system html() conversion of a `latex()' object fails. Follows a
dummy example:
require(Hmisc)
data(Angell)
.object - cor(Angell[,1:2], use=complete.obs)
tmp - latex(.object, cdec=c(2,2), title=)
class(tmp)
[1] latex
html(tmp)
/tmp/RtmprfPwzw/file7e72f7a7.tex:9: Warning: Command not found:
\tabularnewline
Giving up command: \...@hevea@amper
/tmp/RtmprfPwzw/file7e72f7a7.tex:11: Error while reading LaTeX:
This array/tabular column has no specification
Adios
I have hevea-1.10 installed on Debian (according to the help page, it
performs the conversion). Attached is teh produced .tex file. Is this
a bug or would there be a way to work around this behaviour?
Thank you
Liviu
sessionInfo()
R version 2.9.2 (2009-08-24)
x86_64-pc-linux-gnu
locale:
LC_CTYPE=en_GB.UTF-8;LC_NUMERIC=C;LC_TIME=en_GB.UTF-8;LC_COLLATE=en_GB.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_GB.UTF-8;LC_PAPER=en_GB.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_GB.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] datasets grDevices tcltk splines graphics utils stats
[8] methods base
other attached packages:
[1] fortunes_1.3-6 RcmdrPlugin.Export_0.3-0 Rcmdr_1.5-2
[4] car_1.2-15 hints_1.0.1-1boot_1.2-39
[7] relimp_1.0-1 xtable_1.5-5 Hmisc_3.7-0
[10] survival_2.35-7
loaded via a namespace (and not attached):
[1] cluster_1.12.0 grid_2.9.2 lattice_0.17-25 tools_2.9.2
system(uname -a)
Linux debian-liv 2.6.30-1-amd64 #1 SMP Sat Aug 15 18:09:19 UTC 2009
x86_64 GNU/Linux
--
Do you know how to read?
http://www.alienetworks.com/srtest.cfm
Do you know how to write?
http://garbl.home.comcast.net/~garbl/stylemanual/e.htm#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.
--
View this message in context:
http://www.nabble.com/trouble-with-html%28%29-in-Hmisc-tp25721612p25728656.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] trouble with html() in Hmisc
Liviu,
Does this mean that Hmisc produces LyX code (as opposed to LaTeX code)?
No. Sorry for confusing you: it means that html does not know what
tabularnewline means, it cannot interpret it. My reply showed where the
problem lies and how to fix it.
You need to add the following command to the preamble of the *.tex file:
\providecommand{\tabularnewline}{\\}
Regards, Mark.
Liviu Andronic wrote:
Hello
On 10/3/09, Mark Difford mark_diff...@yahoo.co.uk wrote:
This has nothing to do with Hmisc or hevea.
Although I have LyX installed, I don't quite understand where LyX
comes into play. The R code in the original e-mail takes a table-like
object and transforms it into LaTeX; then html() calls hevea to conver
.tex to .html and opens a browser to display the result.
%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}
If I print the `latex' object from the previous e-mail, I get the
offending tabularnewline.
data(Angell)
.object - cor(Angell[,1:2], use=complete.obs)
tmp - latex(.object, cdec=c(2,2), file=, title=)
% latex.default(.object, cdec = c(2, 2), file = , title = )
%
\begin{table}[!tbp]
\begin{center}
\begin{tabular}{lrr}\hline\hline
\multicolumn{1}{l}{}\multicolumn{1}{c}{moral}\multicolumn{1}{c}{hetero}\tabularnewline
\hline
moral$ 1.00$$-0.59$\tabularnewline
hetero$-0.59$$ 1.00$\tabularnewline
\hline
\end{tabular}
\end{center}
\end{table}
Does this mean that Hmisc produces LyX code (as opposed to LaTeX code)?
Liviu
__
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/trouble-with-html%28%29-in-Hmisc-tp25721612p25731240.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] Trendline for a subset of data
Hi Steve, However, I am finding that ... the trendline ... continues to run beyond this data segment and continues until it intersects the vertical axes at each side of the plot. Your best option is probably Prof. Fox's reg.line function in package car. ## library(car) ?reg.line reg.line Regards, Mark. smurray444 wrote: Dear all, I am using abline(lm ...) to insert a linear trendline through a portion of my data (e.g. dataset[,36:45]). However, I am finding that whilst the trendline is correctly displayed and representative of the data portion I've chosen, the line continues to run beyond this data segment and continues until it intersects the vertical axes at each side of the plot. How do I display the line so that it only runs between point 36 and 45 (as shown in the example above) as doesn't continue to display a line throughout the rest of the plot space? Many thanks, Steve _ View your other email accounts from your Hotmail inbox. Add them now. http://clk.atdmt.com/UKM/go/167688463/direct/01/ __ 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/Trendline-for-a-subset-of-data-tp25818425p25821972.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] What is the difference between prcomp and princomp?
Peng Yu wrote: Some webpage has described prcomp and princomp, but I am still not quite sure what the major difference between them is. The main difference, which could be extracted from the information given in the help files, is that prcomp uses the singular value decomposition [i.e. does not rely eigenanalysis], which is generally the preferred method for numerical accuracy. There is plenty of information on the web about the differences between R- and Q-mode PCA. Regards, Mark. Peng Yu wrote: Some webpage has described prcomp and princomp, but I am still not quite sure what the major difference between them is. Can they be used interchangeably? In help, it says 'princomp' only handles so-called R-mode PCA, that is feature extraction of variables. If a data matrix is supplied (possibly via a formula) it is required that there are at least as many units as variables. For Q-mode PCA use 'prcomp'. What are R-mode and Q-mode? Are they just too different numerical methods to compute PCA? __ 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/What-is-the-difference-between-prcomp-and-princomp--tp25952965p25955831.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] Re posting various problems with two-way anova, lme, etc.
Hi Michael, How do you control what is the (intercept) in the model returned by the lme function and is there a way to still be able to refer to all groups and timepoints in there without referring to intercept? Here is some general help. The intercept is controlled by the contrasts that you used when fitting your model. By default these are treatment (i.e. Dunnett) contrasts, but you can change this. ## ?options options(contrasts) ?contrasts options(contrasts=c(contr.sum, contr.poly)) options(contrasts) You can design your own contrasts or you can use those provided by base R: ## ?contr.treatment (There is also contr.sdif in MASS, which is a very useful type for environmental studies) See http://users.monash.edu.au/~murray/stats/BIO4200/LinearModels.pdf, which covers the main ready-made types in R. For ANOVA, sum-to-zero contrasts are commonly used. If you do use treatment contrasts you can (usually) change your reference level on-the-fly or you can set it when you make your factor. ## ?factor ?levels ?reorder ?relevel If you are using the multcomp package then look at: ?contrMat An example of how to use it is given under the glht function. Regards, Mark. Michael Schacht Hansen wrote: Hi, I posted the message below last week, but no answers, so I'm giving it another attempt in case somebody who would be able to help might have missed it and it has now dropped off the end of the list of mails. I am fairly new to R and still trying to figure out how it all works, and I have run into a few issues. I apologize in advance if my questions are a bit basic, I'm also no statistics wizard, so part of my problem my be a more fundamental lack of knowledge in the field. I have a dataset that looks something like this: Week Subj Group Readout 01 A 352.2 11 A 366.6 21 A 377.3 31 A 389.8 41 A 392.5 02 B 344.7 12 B 360.9 . . . . So basically I have a number of subjects that are divided into different groups receiving different interventions/treatments. Observations on these subjects are made on 5 occasions (Week 0-4). I would like to see if there is difference between the treatment groups and if the observations that we are making change in the individual groups over time. According to my very limited statistics training that means that I would have to do a two-way anova with repeated measures because the same subjects are observed on the different weeks. Now in R I can do something like this: MyData$fWeek - factor(MyData$Week) #Convert week to factor m1 - aov(Readout ~ Group*fWeek + Error(Subj/fWeek), data=MyData) My first naive question is: Is the Error(...) term correct for the design that I describe? I am not quite sure if I understand the syntax correctly. In any case, it actually seems to work fine, but I can't figure out how to do post hoc testing on this. TukeyHSD does not work for the aovlist which is returned, so I am kind of stuck. Is there a simple way to do the post hoc test based on the aovlist? I have been reading various questions on the list and I think that I have understood that I should be using lme from the nlme package and this is where I run into some problems. As I understand it, the lme command that I would need would look something like this: m1 - lme(Readout ~ Group*fWeek,random=~1|Subj/fWeek, data=MyData) Now this actually fails with an error message like: Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1 The reason (I believe) is that I have some weeks where there are no observations for a specific group. In this case, one of the groups had a lot of drop-out and at Week 4, there are no subjects left in one of the groups and that seems to be causing the problem. I can get it to run by excluding the last week with something like: m1 - lme(Readout ~ Group*fWeek,random=~1|Subj/fWeek, data=MyData[which(MyData$Week 4), ]) My next question is: Is there another way around that? I would still like to run the analysis for the remaining groups on the last time point and I believe that it should somehow be included into the entire analysis. I have also managed to find a few postings with similar problems, but no real solutions, so anything helpful comments would be much appreciated. My final issue is how do I do the post hoc testing on the model. I understand (I think) that I should use the glht function from the multcomp package. For Instance, I could do something like: summary(glht(m1,linfct=c(GroupB:fWeek1 - GroupC:fWeek1 = 0,GroupB:fWeek2 - GroupC:fWeek2 = 0))) And that actually works fine. My problem is that Group A and fWeek 0 seems to have turned into the (intercept), so I can't write something like: summary(glht(m1,linfct=c(GroupB:fWeek0 - GroupB:fWeek1 = 0))) To check for changes between baseline and week 1 in Group B
Re: [R] Missing data and LME models and diagnostic plots
Peter Flom wrote: I am puzzled by the performance of LME in situations where there are missing data. As I understand it, one of the strengths of this sort of model is how well it deals with missing data, yet lme requires nonmissing data. You are confusing missing data with an unbalanced design. A strengths of LME is that it copes with the latter, which aov() does not. Regards, Mark. Peter Flom-2 wrote: Hello Running R2.9.2 on Windows XP I am puzzled by the performance of LME in situations where there are missing data. As I understand it, one of the strengths of this sort of model is how well it deals with missing data, yet lme requires nonmissing data. Thus, m1.mod1 - lme(fixed = math_1 ~ I(year-2007.5)*TFC_, data = long, random = ~I(year-2007.5)|schoolnum) causes an error in na.fail.default, but adding na.action = na.omit makes a model with no errors. However, if I create that model, i.e. m1.mod1 - lme(fixed = math_1 ~ I(year-2007.5)*TFC_, data = long, random = ~I(year-2007.5)|schoolnum, na.action = na.omit) then the diagnostic plots suggested in Pinheiro Bates produce errors; e.g. plot(m1.mod1, schoolnum~resid(.), abline = 0) gives an error could not find function NaAct. Searching the archives showed a similar question from 2007, but did not show any responses. Thanks for any help Peter ) Peter L. Flom, PhD Statistical Consultant Website: www DOT peterflomconsulting DOT com Writing; http://www.associatedcontent.com/user/582880/peter_flom.html Twitter: @peterflom __ 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/Missing-data-and-LME-models-and-diagnostic-plots-tp25996584p25998546.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] Missing data and LME models and diagnostic plots
Hi Peter, See e.g. Hedeker and Gibbons, Longitudinal Data Analysis, which repeatedly stresses that mixed models provide good estimates if the data are missing at random. This may be true. However, one of the real strengths of LME is that it handles unbalanced designs, which is a different thing. The fact that it also gets good estimates when data are missing is a bonus. Regards, Mark. Peter Flom-2 wrote: I wrote I am puzzled by the performance of LME in situations where there are missing data. As I understand it, one of the strengths of this sort of model is how well it deals with missing data, yet lme requires nonmissing data. Mark Difford replied You are confusing missing data with an unbalanced design. A strengths of LME is that it copes with the latter, which aov() does not. Thanks for your reply, but I don't believe I am confused with respect to missing data in mixed models. See e.g. Hedeker and Gibbons, Longitudinal Data Analysis, which repeatedly stresses that mixed models provide good estimates if the data are missing at random. Regards Peter Peter L. Flom, PhD Statistical Consultant Website: www DOT peterflomconsulting DOT com Writing; http://www.associatedcontent.com/user/582880/peter_flom.html Twitter: @peterflom __ 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/Missing-data-and-LME-models-and-diagnostic-plots-tp25996584p26004465.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] Lost all script
Hi David, Now when I turn on R again the script is now completely blank. This happened to me about 4--5 months ago under Vista. I cannot quite remember what I did but I think I got the script working by opening it in another editor (a hex editor would do) and removing either the first few bytes or the last few bytes. If you still have the file try this route. Good luck, Mark. David Young-18 wrote: Hi all, I just had a rather unpleasant experience. After considerable work I finally got a script working and set it to run. It had some memory allocation problems when I came back so I used Windows to stop it. During that process it told me that the script had been changed and asked if I wanted to save it. Not being positive that I'd saved the very last changes I said yes. Now when I turn on R again the script is now completely blank. I guess my questions are: Is there a way to interrupt a program without using Windows? Is there anyway to recover my script? And a nice to know: Anybody know why it saved blank space as the new script? Thanks for any advice. A humble, and humbled, new R user. -- Best regards, David Young Marketing and Statistical Consultant Madrid, Spain +34 913 540 381 http://www.linkedin.com/in/europedavidyoung mailto:dyo...@telefonica.net __ 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/Lost-all-script-tp26096908p26101731.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] Logistic and Linear Regression Libraries
Hi Phil, So far for logistic regression I've tried glm(MASS) and lrm (Design) and found there is a big difference. Be sure that you mean what you say, that you are saying what you mean, and that you know what you mean when making such statements, especially on this list. glm is not in MASS, so perhaps you mean polr in package MASS. And no, there is no big difference. You are doing something wrong. ## Edited output from polr and lrm house.lrm - lrm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) house.lrm Logistic Regression Model lrm(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq) snip CoefS.E.Wald Z P y=Medium 0.4961 0.12485 3.97 0.0001 y=High-0.6907 0.12547 -5.50 0. Infl=Medium 0.5664 0.10465 5.41 0. Infl=High 1.2888 0.12716 10.14 0. Type=Apartment -0.5724 0.11924 -4.80 0. Type=Atrium-0.3662 0.15517 -2.36 0.0183 Type=Terrace -1.0910 0.15149 -7.20 0. Cont=High 0.3603 0.09554 3.77 0.0002 house.plr - polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) summary(house.plr) Re-fitting to get Hessian Call: polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq) Coefficients: Value Std. Error t value InflMedium 0.5663937 0.10465285 5.412120 InflHigh 1.2888191 0.12715641 10.135699 TypeApartment -0.5723501 0.11923824 -4.800055 TypeAtrium-0.3661866 0.15517362 -2.359851 TypeTerrace -1.0910149 0.15148617 -7.202076 ContHigh 0.3602841 0.09553587 3.771192 snip Regards, Mark. tdm wrote: Hi all, I'm trying to discover the options available to me for logistic and linear regression. I'm doing some tests on a dataset and want to see how different flavours of the algorithms cope. So far for logistic regression I've tried glm(MASS) and lrm (Design) and found there is a big difference. Is there a list anywhere detailing the options available which details the specific algorithms used? Thanks in advance, Phil -- View this message in context: http://old.nabble.com/Logistic-and-Linear-Regression-Libraries-tp26140248p26140375.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] a publication quality table for summary.zeroinfl()
Hi Chris, My ideal would be to gather the information onto the clipboard so I could paste it into Excel and do the formatting there, but any approach would be better than what I have now. I would never use Excel for this since there are far superior tools available. But it is very easy to assemble the two parts into a single table for further manipulation. Not sure about the clipboard route, but the following rather crude approach saves the object (a list) to a *.csv file in your working directory. ## Simple, crude approach, with theta replicated as last column in the saved *.csv data(bioChemists, package = pscl) fm_zinb - zeroinfl(art ~ . | 1, data = bioChemists, dist = negbin) Temptab - list(coef=rbind(summary(fm_zinb)$coefficients[[1]], summary(fm_zinb)$coefficients[[2]]), theta=fm_zinb$theta) getwd() write.csv(Temptab, file=Temptab.csv) read.csv(Temptab.csv) ## If you have a latex distribution library(Hmisc) latex(list(coef=rbind(summary(fm_zinb)$coefficients[[1]], summary(fm_zinb)$coefficients[[2]]), theta=fm_zinb$theta), cdec=c(3,3,3,4,3)) Regards, Mark. Chris Fowler wrote: I will swallow my pride and post to this list for the second time in 24 hours--I have a paper due to a reviewer and I am desperate. Has anybody written code to move the output from summary() called on the results of a zeroinfl() model (from the pscl package) into a form suitable for publication? When I hit send on this message I am going to begin hand typing stars into a spreadsheet. The problem is that the zero-inflated model has two parts: a count and a zero portion--its coefficients are stored in separate arrays and there is a Log(theta) that needs to be thrown in there that is in a completely separate place in the structure of the summary function. As a result the functions that I have found for outputting summaries of other linear models all give error messages when attempted on this summary object. My ideal would be to gather the information onto the clipboard so I could paste it into Excel and do the formatting there, but any approach would be better than what I have now. Thanks, Chris __ 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://old.nabble.com/a-publication-quality-table-for-summary.zeroinfl%28%29-tp26140199p26140542.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] Logistic and Linear Regression Libraries
tdm wrote: OK, I think I've figured it out, the predict of lrm didn't seem to pass it through the logistic function. If I do this then the value is similar to that of lm. Is this by design? Why would it be so? Please take some time to read the help files on these functions so that you at least understand that they fit different models. ?Design:::lrm ?Design:::predict.lrm ?Design:::datadist ?lm This really is just a plainer way of saying what I said before. HTH, Mark. tdm wrote: OK, I think I've figured it out, the predict of lrm didn't seem to pass it through the logistic function. If I do this then the value is similar to that of lm. Is this by design? Why would it be so? 1 / (1 + Exp(-1 * 3.38)) = 0.967 tdm wrote: Anyway, do you know why the lrm predict give me a values of 3.38? -- View this message in context: http://old.nabble.com/Logistic-and-Linear-Regression-Libraries-tp26140248p26141711.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] [Fwd: Re: Plotting log-axis with the exponential base to a plot with the default logarithm base 10]
Elisabeth, You should listen to Ted (Harding). He answered your question with: the vertical axis is scaled logarithmically with the numerical annotations corresponding to the *raw* values of Y, not to their log-transformed values. Therefore it does not matter what base of logarithms is used... And he has tried to help you further by asking you to answer the following question: How is it that you recognise that it is a logaritmic axis with the base of 10, as opposed to any other base? Maybe a little graphic demonstration (run the script below) will help to convince you. The top row of the figure re-poses Ted's question, viz How do you know which logarithmic base was used to plot which sub-figure? You will see from my script that the first is plotted using natural logs, whereas the next two use logs to base 10. (The last of them uses R's log=y argument.) What's the difference? The first two sub-figures on the bottom row show the scales that were hidden in the sub-figures above them (from the script you will see that all I did was turn off the annotation/labelling of the scale). The last sub-figure on the bottom row is identical in __appearance__ to that above it. However, as you will see from my script, it plots data that have been transformed to natural logarithms. For scale-annotation, however, I have used raw data values. That was Ted's first point: it doesn't matter which base of logarithm your data have been transformed to if you annotate the scale using (backtransformed) raw values. ## Show that log=y does the job you want done even if internally it uses logs to the base 10 ## rather than natural logarithms par(mfrow=c(2,3)) plot(log(1:10), yaxt=n, ylab=) plot(log10(1:10), yaxt=n, ylab=) plot((1:10), log=y) plot(log(1:10)) plot(log10(1:10)) plot(log(1:10), yaxt=n) axis(side=2, at=sapply(c(1,2,5,10), FUN=function(x) log(x)), labels=c(1,2,5,10)) Regards, Mark -- View this message in context: http://r.789695.n4.nabble.com/Fwd-Re-Plotting-log-axis-with-the-exponential-base-to-a-plot-with-the-default-logarithm-base-10-tp2172435p2173481.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] [Fwd: Re: Plotting log-axis with the exponential base to
Hi All, You can also add a line using lines() if you transform in the call using the same log-base---but not via R's log=y argument (because of what's stored in par(yaxp)). ## par(mfrow=c(1,3)) plot(1:10, log=y) lines(log10(1:10)) par(yaxp) plot(log10(1:10), yaxt=n) axis(side=2, at=sapply(c(1,2,5,10), FUN=function(x) log10(x)), labels=c(1,2,5,10)) lines(log10(1:10)) par(yaxp) plot(log(1:10), yaxt=n) axis(side=2, at=sapply(c(1,2,5,10), FUN=function(x) log(x)), labels=c(1,2,5,10)) lines(log(1:10)) par(yaxp) Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/Fwd-Re-Plotting-log-axis-with-the-exponential-base-to-a-plot-with-the-default-logarithm-base-10-tp2172435p2182338.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] Multiple comparisons; rank-based anova
Dear List-Members, Is the application of multiple comparison procedures (using the multcomp package) to the output from a rank-based ANOVA straightforward, or do I need to take heed ? That is, is it as simple as: glht( aov(rank(NH4) ~ Site, data=mydat), linfct=mcp(Site=Tukey) ) Thanks in advance for your help. Regards, Mark Difford. -- View this message in context: http://www.nabble.com/Multiple-comparisons--rank-based-anova-tf4516025.html#a12881037 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] Spacing between x axis and labels
Also look at the mgp option under ?par. This allows one to set the margin line for axis title, axis labels and axis line... Regards, Mark Difford. Marc Schwartz wrote: On Tue, 2007-09-25 at 15:39 +0200, Christian Schäfer wrote: Hi, my x-axis contains labels that consist of two lines (the actual label and below information about groupsize). Unfortunately, there is too little spacing between labels and tickmarks in this situation. Is there a parameter to tune spacing between axis and labels? Thanks, Chris See ?mtext and take note of the 'line' argument, which allows you to specify the margin line upon which to place text labels. The 'at' argument is used to define where on the axis the labels should be places. If you are not already, you will want to use either 'axes = FALSE' or 'xaxt = n' in your plot call, so that the default axis labels are not drawn. HTH, Marc Schwartz __ 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/Spacing-between-x-axis-and-labels-tf4516362.html#a12883395 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] mantel tests
Hi Simon, This is a general statistics question so I'm sorry if its outside the field of r help. Could anyone please provide a starting point? This doesn't directly answer your question, but one thing I would do is get ade4 and look at Chessel's co-inertia analysis method. It's a very general, and robust, method. It would involve treating your male and female trait tables separately, and then matching then (usually after a PCA). The method is based on (or uses) Robert Escoufier's RV-coefficient, and looks at the covariance between the two matrices. See: Dray, S., Chessel, D. and J. Thioulouse (2003) Co-inertia analysis and the linking of the ecological data tables. Ecology, 84, 11, 3078–3089 http://biomserv.univ-lyon1.fr/~dray/files/articles/SD162.pdf There is an excellent plot method, which should give you good insight, and there is a method for doing Monte-Carlo permutation tests of the match. I hope this is useful. Regards, Mark Difford. Simon Pickett wrote: This is a general statistics question so I'm sorry if its outside the field of r help. Anyway, I have a suite of female and male traits and I have made a matrix of correlation coefficients using rcorr(). This results in a 6 by 6 matrix like this.. [1] 0.11287990 0.20441361 0.23837442 0.04713234 0.04331637 0.01461611 [7] 0.22627981 0.11720108 0.14252307 0.19531625 0.29989953 0.09989502 [13] 0.03888750 0.11157971 0.02693303 0.01373447 0.08913154 0.06770636 [19] 0.01984838 0.10047978 0.05200218 0.16317234 0.2663 0.10412373 [25] 0.06269722 0.14366454 0.13123054 0.27550149 0.43863848 0.28909831 [31] 0.01454485 0.02551081 0.05645427 0.15819397 0.16508231 0.12399349 I want to test 2 hypotheses 1) is there a pattern to the matrix, does it differ from random? 2) do the top left and bottom right quadrants differ from the other 2 quadrants and if so, which one has the highest values of r2. I have read alot about permutation tests, bootstrapping and mantel tests but cant decide which is best in this situation? Could anyone please provide a starting point? Thankyou in advance, Simon [[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/mantel-tests-tf4743237.html#a13571385 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] structure vs. matrix
Hi Edna, When creating a matrix, is it better to use the structure function or the matrix function...? I hope you have a huge (empty) jar in the kitchen, and that your pantry is empty. R isn't too difficult, except if you're trying to do stats (and don't know what you are doing --- though I am not suggesting that you don't know how to do stats;). ## see ?structure ?matrix A structure is not a matrix. So, if you want to make a matrix, use matrix(). Easy. HTH, Mark. Edna Bell wrote: Hi R Gurus! When creating a matrix, is it better to use the structure function or the matrix function, please? Thanks, Edna Bell __ 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/structure-vs.-matrix-tf4745521.html#a13572508 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] Getting started help
Hi Rthoughts, I am currently discouraged by the use of r. I cannot figure out how to use it despite extensive searches. Can anyone help me with getting started? How can import a txt file with series... There are piles of documents that you could (and should) read. I am surprised that you haven't found them in your extensive searches. To set your thoughts free, go here: http://www.r-project.org/ http://cran.r-project.org/manuals.html## read the first one, and whatever else you need Then click on CRAN, and choose a mirror site (nearby), say: http://cran.za.r-project.org/ Then go: http://cran.za.r-project.org/other-docs.html and start at the top. HTH, Mark. PS: It might help you in your future enquiries if you used a real name, and if you called R R rather than r (since your keyboard clearly doesn't have a broken R ;). Rthoughts wrote: Hello, I need to learn to use r-software for my PhD research project involving long timelines of radon radiation variations. I am using windows. I am currently discouraged by the use of r. I cannot figure out how to use it despite extensive searches. Can anyone help me with getting started? How can import a txt file with series of numbers from the RAD7 machine? How can I open up and set directories with the imported file or to start a new r session? Thank you so much if you can help. -- View this message in context: http://www.nabble.com/Getting-started-help-tp15560581p15561419.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] Getting started help
Hi Rthoughts, It isn't clear what you mean. When you install R, the installation program usually puts an icon on your desktop that you can click on to run the program. So, if you don't have that, but have installed R, and what you mean is, How do I start the R program? or How do I run R? then do the following: Go to the directory where R has been installed. If necessary, use your file manager to find the file named Rgui.exe. When you have found it, execute it to start an R session. (Rgui.exe usually lives in something like: ...\Program Files\R\bin\) Later you will find out how to set this all up properly. Usually the installation program does a faultless job. HTH, Mark. Rthoughts wrote: Hi Mark, Thank you for your reply. There is one link I haven't come across, the last one. I have seen them but I couldn't find where 'how to start R' is explained for Windows platforms. I will look further into them. As for everyone else who sent e-mails, thank you. I have printed them out and will look into them. Mark Difford wrote: Hi Rthoughts, I am currently discouraged by the use of r. I cannot figure out how to use it despite extensive searches. Can anyone help me with getting started? How can import a txt file with series... There are piles of documents that you could (and should) read. I am surprised that you haven't found them in your extensive searches. To set your thoughts free, go here: http://www.r-project.org/ http://cran.r-project.org/manuals.html## read the first one, and whatever else you need Then click on CRAN, and choose a mirror site (nearby), say: http://cran.za.r-project.org/ Then go: http://cran.za.r-project.org/other-docs.html and start at the top. HTH, Mark. PS: It might help you in your future enquiries if you used a real name, and if you called R R rather than r (since your keyboard clearly doesn't have a broken R ;). Rthoughts wrote: Hello, I need to learn to use r-software for my PhD research project involving long timelines of radon radiation variations. I am using windows. I am currently discouraged by the use of r. I cannot figure out how to use it despite extensive searches. Can anyone help me with getting started? How can import a txt file with series of numbers from the RAD7 machine? How can I open up and set directories with the imported file or to start a new r session? Thank you so much if you can help. -- View this message in context: http://www.nabble.com/Getting-started-help-tp15560581p15561960.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] Getting started help
Hi Rthoughts, Yes, I see now that they truly are (just) Rthoughts;) but take courage, for we are getting closer (to the start). You still need to read the basic documentation, and you will get used to the command line. What I think you need is a package called Rcmdr. So, start R using your desktop icon, click on the Packages menu, then on Install package(s)..., and then scroll down and choose Rcmdr (you will probably first be asked to choose a download site). For you the TeachingDemos plugin is a good additional choice. Click OK. Once Rcmdr has been installed, type the following at R's command prompt: require(Rcmdr)## loads Rcmdr and runs it You can now switch windows between Rcmdr and R. Again, at the R prompt type the following commands: require(datasets) ?datasets data(ChickWeight) ## careful: R is case sensitive! ls() ## try ?ls Now, switch to the Rcmdr window and use the Data set: button to load ChickWeight into Rcmdr's workspace. Now mess with it; or load another data set, more to your taste. Use Rcmdr's Tools menu to load any plugins you installed. ls() str(ChickWeight)## try ?str summary(ChickWeight) ?rm rm(ChickWeight) ls() HTH, Mark. Rthoughts wrote: Hi Mark, Thank you for the reply. I meant the command prompts to start an R file. To be followed on by importint data I can then use to practise the software with. The installation did put an icon on teh desktop. I am a very skilled user of computers but command lines for many programs is something that cna throw me sometimes! Regrads, Seb. Mark Difford wrote: Hi Rthoughts, It isn't clear what you mean. When you install R, the installation program usually puts an icon on your desktop that you can click on to run the program. So, if you don't have that, but have installed R, and what you mean is, How do I start the R program? or How do I run R? then do the following: Go to the directory where R has been installed. If necessary, use your file manager to find the file named Rgui.exe. When you have found it, execute it to start an R session. (Rgui.exe usually lives in something like: ...\Program Files\R\bin\) Later you will find out how to set this all up properly. Usually the installation program does a faultless job. HTH, Mark. Rthoughts wrote: Hi Mark, Thank you for your reply. There is one link I haven't come across, the last one. I have seen them but I couldn't find where 'how to start R' is explained for Windows platforms. I will look further into them. As for everyone else who sent e-mails, thank you. I have printed them out and will look into them. Mark Difford wrote: Hi Rthoughts, I am currently discouraged by the use of r. I cannot figure out how to use it despite extensive searches. Can anyone help me with getting started? How can import a txt file with series... There are piles of documents that you could (and should) read. I am surprised that you haven't found them in your extensive searches. To set your thoughts free, go here: http://www.r-project.org/ http://cran.r-project.org/manuals.html## read the first one, and whatever else you need Then click on CRAN, and choose a mirror site (nearby), say: http://cran.za.r-project.org/ Then go: http://cran.za.r-project.org/other-docs.html and start at the top. HTH, Mark. PS: It might help you in your future enquiries if you used a real name, and if you called R R rather than r (since your keyboard clearly doesn't have a broken R ;). Rthoughts wrote: Hello, I need to learn to use r-software for my PhD research project involving long timelines of radon radiation variations. I am using windows. I am currently discouraged by the use of r. I cannot figure out how to use it despite extensive searches. Can anyone help me with getting started? How can import a txt file with series of numbers from the RAD7 machine? How can I open up and set directories with the imported file or to start a new r session? Thank you so much if you can help. -- View this message in context: http://www.nabble.com/Getting-started-help-tp15560581p15562350.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] Why does plot() ignore the data type for axis labels?
Hi Stiffler, I was wondering why the plot() command ignores the datatype when displaying axis labels... plot() doesn't ignore the datatype: x - as.integer(c(1,2,3)) y -x typeof(x) [1] integer mode(x) [1] numeric plot(x,y) calls xy.coords(), which recasts x as: x = as.double(x), which is fine, since x is (also/primarily) numeric. ??? See ?double, sub: Note on names. HTH, Mark. Stiffler wrote: Hello, I was wondering why the plot() command ignores the datatype when displaying axis labels. More specifically, if the data points are integers then the axis labels should intuitively also be integers, right? x - as.integer(c(1,2,3)) y -x typeof(x) [1] integer plot(x,y) The axis labels are 1.0, 1.5, 2.0, 2.5, 3.0 but if the integer type were taken into account they would be 1, 2, 3. PS what's the right way to get integer labels? Z. -- View this message in context: http://www.nabble.com/Why-does-plot%28%29-ignore-the-data-type-for-axis-labels--tp15562325p15567499.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] fitted values are different from manually calculating
Hi Yianni, This just proves that you should be using R as your calculator, and not the other one! Regards, Mark. gatemaze wrote: Hello, on a simple linear model the values produced from the fitted(model) function are difference from manually calculating on calc. Will anyone have a clue... or any insights on how fitted function calculates the values? Thank you. -- -- Yianni [[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/fitted-values-are-different-from-manually-calculating-tp15568106p15569205.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] plotDensity
Hi Conny, It still isn't clear what your question is, but a density plot simply shows you the distribution of your data, say a set of measurements of something. Think of it as a modern replacement for the histogram. See http://en.wikipedia.org/wiki/Density_estimation for greater insight. HTH, Mark. laptopcss wrote: Hallo, I just loaded the limma package. Conny Original-Nachricht Datum: Tue, 19 Feb 2008 06:41:18 -0800 Von: Henrik Bengtsson [EMAIL PROTECTED] An: [EMAIL PROTECTED] CC: [EMAIL PROTECTED] Betreff: Re: [R] plotDensity Is 'plotDensity' a specific function you are referring to? If so, in which package? Then we can give a straight answer... /Henrik On Feb 19, 2008 4:10 AM, [EMAIL PROTECTED] wrote: Hallo, I have a question to plotDensity and do not understand what stand behind. In a density plot the x-axis is the signal intensity but for what stands than density on the y-axis? Here I have the values 0.00-0.30 Can anyone discribe it by his own words? I do not understand the help. Thanks, Conny -- __ 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. -- View this message in context: http://www.nabble.com/plotDensity-tp15565443p15570676.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] Why does plot() ignore the data type for axis labels?
Hi Stiffler, Duncan Murdoch wrote: The point is that a scatterplot is a graph of real numbers, not of integers. Use a plot designed for a discrete data type if you don't want things displayed as their real values ... To drive home Duncan's point (not that it needs it) ...and to complete the loop, consider how R treats factors: x [1] 1 2 3 factor(x) [1] 1 2 3 Levels: 1 2 3 ## And: factor(x, labels=c(one, two, three)) [1] one two three Levels: one two three unclass(factor(x, labels=c(one, two, three))) [1] 1 2 3 attr(,levels) [1] one two three It does, in the end, make perfect sense. HTH, Mark. Duncan Murdoch-2 wrote: On 19/02/2008 5:40 PM, Stiffler wrote: Mark Difford wrote: I was wondering why the plot() command ignores the datatype when displaying axis labels... plot() doesn't ignore the datatype: [...] plot(x,y) calls xy.coords(), which recasts x as: x = as.double(x), which is fine, since x is (also/primarily) numeric. HTH, Mark. Thanks for the explanation Mark. If integers are being recast as doubles R is ignoring/overriding the user's choice of data-type. There may be good reasons for doing that internally (uniformity, code re-use etc) , but it is not what I'd expect as an end-user --- neither ?plot nor ?xy.coords seem to mention that coordinates need to be floating point numbers. They don't need to be floating point numbers, they are converted if not. The point is that a scatterplot is a graph of real numbers, not of integers. Use a plot designed for a discrete data type if you don't want things displayed as their real values (see for example barplot, stripchart, dotchart, etc.) 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. -- View this message in context: http://www.nabble.com/Why-does-plot%28%29-ignore-the-data-type-for-axis-labels--tp15562325p15584404.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] Mixed model Nested ANOVA
Hi Stephen, Hopefully you will get an answer from one of the experts on mixed models who subscribe to this list. However, you should know that both lme() and lmer() currently have anova() methods. The first will give you p-values (but no SS), and the second will give you SS (but no p-values). You can, however, get the latter using functions in the languageR package. This also has an aovlmer.fnc() that uses MCMC to get p-values. From what you have said about your data, and about what you want from them, you clearly should be using either lme() or lmer(). Further, objects from both functions work with glht() from the multcomp package, so you also have access to a full range of post-hoc tests. HTH, Mark. Stephen Cole-2 wrote: hello R help I am trying to analyze a data set that has been collected from a hierarchical sampling design. The model should be a mixed model nested ANOVA. The purpose of my study is to analyze the variability at each spatial scale in my design (random factors, variance components), and say something about the variability between regions (fixed factor, contrast of means). The data is as follows; region (fixed) Location (random) Site(random) site nested in location nested in region. I would like to run this as an ANOVA and then compute variance components. My question is when i use the aov command; mod1 - aov(density ~ region/location/site) is there any way to tell R which factor is random and fixed. I know i can specify fixed and random factors using lme or lmer but these methods do not allow me to extract an anova table (or do they?) I know that the data can be analyzed using a nested ANOVA because i have based my design on several papers in the marine biological literature (MEPS). Thank-you for any advice in advance. Stephen Cole [[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/Mixed-model-Nested-ANOVA-tp15639930p15641867.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] mixed model nested ANOVA (part two)
Hi Stephen, Also i have read in Quinn and Keough 2002, design and analysis of experiments for biologists, that a variance component analysis should only be conducted after a rejection of the null hypothesis of no variance at that level. Once again the caveat: there are experts on this list who really know about this stuff, and I am not one of them. Your general strategy would be to set up two models with the same fixed effects, one of which doesn't have random effects. You then test the two models using anova(mod.withRandom, modWithoutRandom). I haven't tried this using lmer/2(), but with lme() you do this by fitting your fixed+random effects model using lme() and your fixed-only effects model using lm(). If you are using weights to model heteroskedasticity, then it's better to use gls(), as this will accept the same weights argument as the call to lme(). Then you simply do anova(lme.model, lm/gls.model). This tells you about the significance of your random effects, i.e. whether you need a random-effects component. To recap: mod.rand - lme(fixed=y ~ x, random=~x|Site, data=...) mod,fix - lm(fixed=y ~ x, data=...) anova(mod.rand, mod.fix) HTH, Mark. Stephen Cole-2 wrote: First of all thank you for the responses. I appreciate the suggestions i have received thus far. Just to reiterate I am trying to analyze a data set that has been collected from a hierarchical sampling design. The model should be a mixed model nested ANOVA. The purpose of my study is to analyze the variability at each spatial scale in my design (random factors, variance components), and say something about the variability between regions (fixed factor, contrast of means). The data is as follows; region (fixed) Location (random) Site(random) site nested in location nested in region. Also i have read in Quinn and Keough 2002, design and analysis of experiments for biologists, that a variance component analysis should only be conducted after a rejection of the null hypothesis of no variance at that level. I have tried to implement mod1-lmer(density ~ 1 + (1|site) + (1|location) + (1|region)) However, as i understand it, this treats all my factors as random. Plus I do not know how to extract SS or MS from this model. anova(mod1) gives me Analysis of Variance Table Df Sum Sq Mean Sq and summary(mod1) gives me Linear mixed-effects model fit by REML Formula: density ~ 1 + (1 | site) + (1 | location) + (1 | region) AIC BIC logLik MLdeviance REMLdeviance 15658 15678 -7825 1566215650 Random effects: Groups NameVariance Std.Dev. site (Intercept) 22191 148.97 location (Intercept) 33544 183.15 region (Intercept) 41412 203.50 Residual 696189 834.38 number of obs: 960, groups: site, 4; location, 4; region, 3 Fixed effects: Estimate Std. Error t value (Intercept)261.3 168.7 1.549 from what i understand the variance in the penultimate column are my variance components. But how do i conduct my significance test? I have also tried mod1-lmer(density ~ region + (1|site) + (1|location)) Which i think is the correct mixed model for my design. However once again i do not know how to evaluate significance for the random factors. Thank-you again for any additional advice i receive Stephen Cole __ 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/mixed-model-nested-ANOVA-%28part-two%29-tp15665478p15669384.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] mixed model nested ANOVA (part two)
Hi Stephen, Slip of the dactylus: lm() does not, of course, take a fixed=arg. So you need To recap: mod.rand - lme(fixed=y ~ x, random=~x|Site, data=...) mod,fix - lm(y ~ x, data=...) ## or ##mod,fix - lm(formula=y ~ x, data=...) Bye. Mark Difford wrote: Hi Stephen, Also i have read in Quinn and Keough 2002, design and analysis of experiments for biologists, that a variance component analysis should only be conducted after a rejection of the null hypothesis of no variance at that level. Once again the caveat: there are experts on this list who really know about this stuff, and I am not one of them. Your general strategy would be to set up two models with the same fixed effects, one of which doesn't have random effects. You then test the two models using anova(mod.withRandom, modWithoutRandom). I haven't tried this using lmer/2(), but with lme() you do this by fitting your fixed+random effects model using lme() and your fixed-only effects model using lm(). If you are using weights to model heteroskedasticity, then it's better to use gls(), as this will accept the same weights argument as the call to lme(). Then you simply do anova(lme.model, lm/gls.model). This tells you about the significance of your random effects, i.e. whether you need a random-effects component. To recap: mod.rand - lme(fixed=y ~ x, random=~x|Site, data=...) mod,fix - lm(fixed=y ~ x, data=...) anova(mod.rand, mod.fix) HTH, Mark. Stephen Cole-2 wrote: First of all thank you for the responses. I appreciate the suggestions i have received thus far. Just to reiterate I am trying to analyze a data set that has been collected from a hierarchical sampling design. The model should be a mixed model nested ANOVA. The purpose of my study is to analyze the variability at each spatial scale in my design (random factors, variance components), and say something about the variability between regions (fixed factor, contrast of means). The data is as follows; region (fixed) Location (random) Site(random) site nested in location nested in region. Also i have read in Quinn and Keough 2002, design and analysis of experiments for biologists, that a variance component analysis should only be conducted after a rejection of the null hypothesis of no variance at that level. I have tried to implement mod1-lmer(density ~ 1 + (1|site) + (1|location) + (1|region)) However, as i understand it, this treats all my factors as random. Plus I do not know how to extract SS or MS from this model. anova(mod1) gives me Analysis of Variance Table Df Sum Sq Mean Sq and summary(mod1) gives me Linear mixed-effects model fit by REML Formula: density ~ 1 + (1 | site) + (1 | location) + (1 | region) AIC BIC logLik MLdeviance REMLdeviance 15658 15678 -7825 1566215650 Random effects: Groups NameVariance Std.Dev. site (Intercept) 22191 148.97 location (Intercept) 33544 183.15 region (Intercept) 41412 203.50 Residual 696189 834.38 number of obs: 960, groups: site, 4; location, 4; region, 3 Fixed effects: Estimate Std. Error t value (Intercept)261.3 168.7 1.549 from what i understand the variance in the penultimate column are my variance components. But how do i conduct my significance test? I have also tried mod1-lmer(density ~ region + (1|site) + (1|location)) Which i think is the correct mixed model for my design. However once again i do not know how to evaluate significance for the random factors. Thank-you again for any additional advice i receive Stephen Cole __ 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/mixed-model-nested-ANOVA-%28part-two%29-tp15665478p15669608.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] Hmisc xYplot won't do conditioning on factors?
Hi Ivan, It appears that xYplot, unlike standard xyplot (or coplot to that matter) does not accept factors as x variable in formula. To add to what you have said. It may not be too well documented in ?xYplot, but xYplot() is really designed to do a number of very useful things with two sets of _numeric_ variables, viz x and y. The help file for xYplot() does also mention the subfunction numericScale(), about which it says: numericScale is a utility function that facilitates using xYplot to plot variables that are not considered to be numeric but which can readily be converted to numeric using as.numeric(). This does what you want, and a little bit more (for timeDate variables). It also documents what you say is not documented;) HTH, Mark. Ivan Adzhubey wrote: Hi, Since nobody replied, I'd rather post this self-reply for the sake of documenting the solution. It appears that xYplot, unlike standard xyplot (or coplot to that matter) does not accept factors as x variable in formula. With x converted to numeric everything worked as expected. This small discrepancy was not documented on xYplot help page. --Ivan On Tuesday 26 February 2008 07:47:14 pm Ivan Adzhubey wrote: Hi, I am trying to replace (lattice) standard xyplot with xYplot variant from Hmisc package to be able to add error bars to my plots. However, this does not work, e.g: library(lattice) d - data.frame( SKU=gl(3, 1, 21, labels=c(a, b, c)), Weekday=gl(7, 3, 21), QCRate=runif(21)) xyplot(QCRate ~ Weekday | SKU, data=d) (this plots nice 3 panels as per a,b,c conditionals) library(Hmisc) xYplot(QCRate ~ Weekday | SKU, data=d) Error in Summary.factor(1:7, na.rm = TRUE) : range not meaningful for factors Is there a workaround? Thanks, Ivan __ 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/Hmisc-xYplot-won%27t-do-conditioning-on-factors--tp15703668p15760460.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] model R^2 and partial R^2 values
Hi Mirela, Are the relative R^2 values the CP values? No. CP is your complexity parameter. I’ve read that the R^2 = 1-rel error, so I am assuming that in my case this would be 1-0.64949. Is this correct? Yes. See ?rsq.rpart, and run the example, which I've copied below. ## par(ask=T) z.auto - rpart(Mileage ~ Weight, car.test.frame) rsq.rpart(z.auto) par(ask=F) The values plotted in graph-1 come from 1-rel.error and 1-xerror. HTH, Mark. Mirela Tulbure wrote: Dear R-list members, I am doing a CART analysis in R using the rpart function in the rpart package: Phrag.rpart=rpart(PhragDiff~., data = Phrag, method=anova, xval=10). I used the xerror values in the CP table to prune the tree to 4 nsplits: CPnsplit rel error xerrorxstd 1 0.098172 0 1.0 1.02867 0.12768 2 0.055991 3 0.70548 1.00823 0.12911 3 0.029306 4 0.64949 0.83275 0.12074 4 0.018943 5 0.62019 0.86994 0.12467 5 0.010503 6 0.60124 0.86975 0.12080 6 0.01 7 0.59074 0.87944 0.11757 I would like to get R^2 values for the model as well as partial R^2 values for each split. I’ve read that the R^2 = 1-rel error, so I am assuming that in my case this would be 1-0.64949. Is this correct? Are the relative R^2 values the CP values? Your help would be much appreciated. Many thanks in advance, Mirela Looking for la ols.search.yahoo.com/newsearch/category.php?category=shopping __ 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/model-R%5E2-and-partial-R%5E2-values-tp15772594p15773072.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] [OT] Make plots with GNUplot. Have anyone tried that?
Hi Everyone, Please don't denigrate the capabilities of GNUplot (Louise excluded). It can, in fact, do some truly awesome stuff. http://linuxgazette.net/133/luana.html The PDF is worth a shot. Cheers, Mark. Louise Hoffman-3 wrote: If you still want to then read ?write.table, that can export your data into a spreadsheet-like ascii format which can be used from GNUplot easily. Very interesting. So if I e.g. write: ts.sim - arima.sim(list(order = c(1,1,0), ar = 0.7), n = 200) ts.plot(ts.sim) How do I know the names of the rows to put in the data.frame() command? Btw, comparing the graphics capabilities of GNUplot and R, it is something like a three-wheel bicycle and a spaceship. Guess which is which. =) I know that I will most likely spend a lot of time on just making the plots, but I atleast (for now =) ) think it could be fun to try. __ 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/Make-plots-with-GNUplot.-Have-anyone-tried-that--tp15767196p15773081.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] The Claw Density and LOCFIT
Hi Jaap, Could anybody please direct me in finding an updated version of this document, or help me correct the code given in the file. The (out-of-date) code is as follows: You are not helping yourself, or anyone else, by not including the error messages you get when trying to execute your code. The matter is in fact quite straightfowards: ev does not accept an evaluation structure called grid. The current documentation for locfit does tell you this: ?locfit.raw (sub ev). When I executed your code I got fit1 - locfit( ~ claw54, deg = 0, kern = gauss, alpha = c(0, 0.315), ev = + grid(100, ll = -3.5, ur = 2.7)) Error in grid(100, ll = -3.5, ur = 2.7) : unused argument(s) (ll = -3.5, ur = 2.7) This explicitly tells you that the problem lies with the call to ev = grid(...), which should in fact be ev = lfgrid(...). The following works for me and should do so for you. fit1 - locfit( ~ claw54, deg = 0, kern = gauss, alpha = c(0, 0.315), ev = lfgrid(100, ll = -3.5, ur = 2.7)) HTH, Mark. R version 2.10.0 Under development (unstable) (2009-06-20 r48806) i386-pc-mingw32 locale: [1] LC_COLLATE=English_South Africa.1252 LC_CTYPE=English_South Africa.1252 [3] LC_MONETARY=English_South Africa.1252 LC_NUMERIC=C [5] LC_TIME=English_South Africa.1252 attached base packages: [1] splines stats graphics grDevices utils datasets methods base other attached packages: [1] locfit_1.5-4lattice_0.17-25 akima_0.5-2 DPpackage_1.0-7 ade4_1.4-11 Design_2.2-0 [7] survival_2.35-4 Hmisc_3.6-0 loaded via a namespace (and not attached): [1] cluster_1.12.0 grid_2.10.0tools_2.10.0 Van Wyk, Jaap wrote: I am trying to reproduce Figure 10.5 of Loader's book: Local Regression and Likelihood. The code provided in the book does not seem to work. I have managed (a while ago) to get the accompanied R-code for the figures in the book (file called lffigs.R) from somewhere - cannot find it on the web anymore. The code in the .R script file does not work either. Could anybody please direct me in finding an updated version of this document, or help me correct the code given in the file. The (out-of-date) code is as follows: data(claw54) fit1 - locfit( ~ claw54, deg = 0, kern = gauss, alpha = c(0, 0.315), ev = grid(100, ll = -3.5, ur = 2.7)) fit2 - locfit( ~ claw54, deg = 0, kern = gauss, alpha = c(0, 0.985), ev = grid(100, ll = -3.5, ur = 2.7)) x - seq(-3.5, 2.7, length.out = 200) y - dnorm(x, -1., 0.1) + dnorm(x, -0.5, 0.1) + dnorm(x, 0, 0.1) + dnorm(x, 0.5, 0.1) + dnorm(x, 1., 0.1) y - (y + 5 * dnorm(x))/10 plot(fit1, get.data = T, main = h=0.315, ylim = c(0, max(y))) lines(x, y, lty = 2) plot(fit2, get.data = T, main = h=0.985, ylim = c(0, max(y))) lines(x, y, lty = 2) THANKS FOR ANY ASSISTANCE. ps: This code differs from that in the book. I have tried both, without success. Even if I just use, for example, fit1 - locfit( ~ claw54, deg = 0, kern = gauss, alpha = c(0, 0.315)) I do not get the same result. Jacob L van Wyk, Dept. of Statistics, University of Johannesburg (APK), Box 524, Auckland Park, 2006. Office: +27 11 559 3080, Fax: +27 11 559 2499 This email and all contents are subject to the following disclaimer: http://www.uj.ac.za/UJ_email_legal_disclaimer.htm [[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/The-Claw-Density-and-LOCFIT-tp24218274p24220007.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] Correlation Network Diagram?
Hi Rory, There are several. Have a look at the gR Task Views. There you will also find a link to the statnet suite, where you will find links to a dedicated set of jstatsoft articles. Regards, Mark. Rory Winston wrote: Hi all On page 39 of this paper [1] by Andrew Lo there is a very interesting correlation network diagram (sorry I dont have a more convenient link to the type of diagram I'm talking about). does anyone know of any package in R that can generate these types of diagrams? Cheers -- Rory [1] http://web.mit.edu/alo/www/Papers/august07.pdf [[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/Correlation-Network-Diagram--tp24321104p24329465.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] Differing Variable Length Inconsistencies in Random Effects/Regression Models
Hi Aditi,
Parts of _your_ code for the solution offered by Jerome Goudet are wrong;
see my comments.
famfit-lmer(peg.no~1 + (1|family), na.action=na.omit, vcdf) ## use:
na.action=na.exclude
resfam-residuals(famfit)
for( i in 1:length(colms))
+ {
+ print (Marker, i)
+ regfam-abline(lm(resfam~i))## you need to use:
abline(lm(resfam~colms[,i]))
+ print(regfam)
+ }
Corrected code:
famfit-lmer(peg.no~1 + (1|family), na.action=na.exclude, vcdf)
resfam-residuals(famfit)
for( i in 1:length(colms))
{
print (Marker, i)
regfam-abline(lm(resfam~colms[,i]))
}
This should work.
Regards, Mark.
A Singh wrote:
Dear All,
I am quite new to R and am having a problem trying to run a linear model
with random effects/ a regression- with particular regard to my variable
lengths being different and the models refusing to compute any further.
The codes I have been using are as follows:
vc-read.table(P:\\R\\Testvcomp10.txt,header=T)
attach(vc)
family-factor(family)
colms-(vc)[,4:13] ## this to assign the 10 columns containing marker
datato a new variable, as column names are themselves not in any
recognisable sequence
vcdf-data.frame(family,peg.no,ec.length,syll.length,colms)
library(lme4)
for (c in levels(family))
+ {for (i in 1:length(colms))
+{ fit-lmer(peg.no~1 + (1|c/i), vcdf)
+}
+summ-summary(fit)
+av-anova(fit)
+print(summ)
+print(av)
+ }
This gives me:
Error in model.frame.default(data = vcdf, formula = peg.no ~ 1 + (1 + :
variable lengths differ (found for 'c')
I had posted a similar message on the R mixed model list a few days ago,
with respect to my fundamental methods, and Jerome Goudet had kindly
referred me to an alternative approach using residuals obtained from a
random effects model in lmer(), and then doing regressions using those
[residuals being the dependent variable and my marker data columns the
independent variable].
The code for that is as follows:
vc-read.table(P:\\R\\Text
Files\\Testvcomp10.txt,header=T,sep=,dec=.,na.strings=NA,strip.white=T)
attach(vc)
family-factor(family)
colms-(vc)[,4:13]
names(vc)
[1] male.parent family offspring.id P1L55P1L73
[6] P1L74P1L77P1L91P1L96P1L98
[11] P1L100 P1L114 P1L118 peg.no
ec.length
[16] syll.length
vcdf-data.frame(family, colms, peg.no, ec.length, syll.length)
library(lme4)
famfit-lmer(peg.no~1 + (1|family), na.action=na.omit, vcdf)
resfam-residuals(famfit)
for( i in 1:length(colms))
+ {
+ print (Marker, i)
+ regfam-abline(lm(resfam~i))
+ print(regfam)
+ }
This again gives me the error:
[1] Marker
Error in model.frame.default(formula = resfam ~ i, drop.unused.levels =
TRUE) :
variable lengths differ (found for 'i')
My variables all have missing values somewhere or the other. The missing
values are not consistent for all individuals, i.e some individuals have
some values missing, others have others,
and as much as I have tried to use na.action=na.omit and na.rm=T, the
differing variable length problem is dogging me persistently..
I also tried to isolate the residuals, save them in a new variable (called
'resfam' here), and tried to save that in the data.frame()-vcdf, that I
had created earlier.
The problem with that was that when the residuals were computed, lmer()
ignored missing data in 'peg.no' with respect to 'family', which is
obviously not the same data missing for say another variable E.g.
'ec.length'- leading again to an inconsistency in variable lengths.
Data.frame would then not accept that addition at all to the previous set.
This was fairly obvious right from the start, but I decided to try it
anyway. Didn't work.
I apologise if the solution to working with these different variable
lengths is obvious and I don't know it- but I don't know R that well at
all.
My data files can be downloaded at the following location:
http://www.filesanywhere.com/fs/v.aspx?v=896d6b88616173be71ab (excel-
.xlsx)
http://www.filesanywhere.com/fs/v.aspx?v=896d6b88616174a76e9e
(.txt file)
Any pointers would be greatly appreciated, as this is holding me up loads.
Thanks a ton for your help,
Aditi
--
A Singh
aditi.si...@bristol.ac.uk
School of Biological Sciences
University of Bristol
__
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/Differing-Variable-Length-Inconsistencies-in-Random-Effects-Regression-Models-tp24502146p24505495.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org
Re: [R] Differing Variable Length Inconsistencies in Random Effects/Regression Models
Perhaps I should have added the following: To see that it works, run the
following:
famfit-lmer(peg.no~1 + (1|family), na.action=na.exclude, vcdf)
resfam-residuals(famfit)
for( i in 1:length(colms))
{
print(coef(lm(resfam~colms[,i])))
}
Regards, Mark.
A Singh wrote:
Dear All,
I am quite new to R and am having a problem trying to run a linear model
with random effects/ a regression- with particular regard to my variable
lengths being different and the models refusing to compute any further.
The codes I have been using are as follows:
vc-read.table(P:\\R\\Testvcomp10.txt,header=T)
attach(vc)
family-factor(family)
colms-(vc)[,4:13] ## this to assign the 10 columns containing marker
datato a new variable, as column names are themselves not in any
recognisable sequence
vcdf-data.frame(family,peg.no,ec.length,syll.length,colms)
library(lme4)
for (c in levels(family))
+ {for (i in 1:length(colms))
+{ fit-lmer(peg.no~1 + (1|c/i), vcdf)
+}
+summ-summary(fit)
+av-anova(fit)
+print(summ)
+print(av)
+ }
This gives me:
Error in model.frame.default(data = vcdf, formula = peg.no ~ 1 + (1 + :
variable lengths differ (found for 'c')
I had posted a similar message on the R mixed model list a few days ago,
with respect to my fundamental methods, and Jerome Goudet had kindly
referred me to an alternative approach using residuals obtained from a
random effects model in lmer(), and then doing regressions using those
[residuals being the dependent variable and my marker data columns the
independent variable].
The code for that is as follows:
vc-read.table(P:\\R\\Text
Files\\Testvcomp10.txt,header=T,sep=,dec=.,na.strings=NA,strip.white=T)
attach(vc)
family-factor(family)
colms-(vc)[,4:13]
names(vc)
[1] male.parent family offspring.id P1L55P1L73
[6] P1L74P1L77P1L91P1L96P1L98
[11] P1L100 P1L114 P1L118 peg.no
ec.length
[16] syll.length
vcdf-data.frame(family, colms, peg.no, ec.length, syll.length)
library(lme4)
famfit-lmer(peg.no~1 + (1|family), na.action=na.omit, vcdf)
resfam-residuals(famfit)
for( i in 1:length(colms))
+ {
+ print (Marker, i)
+ regfam-abline(lm(resfam~i))
+ print(regfam)
+ }
This again gives me the error:
[1] Marker
Error in model.frame.default(formula = resfam ~ i, drop.unused.levels =
TRUE) :
variable lengths differ (found for 'i')
My variables all have missing values somewhere or the other. The missing
values are not consistent for all individuals, i.e some individuals have
some values missing, others have others,
and as much as I have tried to use na.action=na.omit and na.rm=T, the
differing variable length problem is dogging me persistently..
I also tried to isolate the residuals, save them in a new variable (called
'resfam' here), and tried to save that in the data.frame()-vcdf, that I
had created earlier.
The problem with that was that when the residuals were computed, lmer()
ignored missing data in 'peg.no' with respect to 'family', which is
obviously not the same data missing for say another variable E.g.
'ec.length'- leading again to an inconsistency in variable lengths.
Data.frame would then not accept that addition at all to the previous set.
This was fairly obvious right from the start, but I decided to try it
anyway. Didn't work.
I apologise if the solution to working with these different variable
lengths is obvious and I don't know it- but I don't know R that well at
all.
My data files can be downloaded at the following location:
http://www.filesanywhere.com/fs/v.aspx?v=896d6b88616173be71ab (excel-
.xlsx)
http://www.filesanywhere.com/fs/v.aspx?v=896d6b88616174a76e9e
(.txt file)
Any pointers would be greatly appreciated, as this is holding me up loads.
Thanks a ton for your help,
Aditi
--
A Singh
aditi.si...@bristol.ac.uk
School of Biological Sciences
University of Bristol
__
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/Differing-Variable-Length-Inconsistencies-in-Random-Effects-Regression-Models-tp24502146p24506118.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] Split plot analysis problems
Hi Jean-Paul, ... since R is not able to extract residuals? R can extract the residuals, but they are a hidden in models with an error structure ## str(aov(PH~Community*Mowing*Water + Error(Block))) residuals(aov(PH~Community*Mowing*Water + Error(Block))$Block) residuals(aov(PH~Community*Mowing*Water + Error(Block))$Within) Regards, Mark. Jean-Paul Maalouf wrote: Hello, I would be very grateful if someone could give me a hand with my split plot design problems. So here is my design : I am studying the crossed-effects of water (wet/dry) and mowing (mowed/not-mowed = nm) on plant height (PH) within 2 types of plant communities (Xerobromion and Mesobromion) : - Within each type of communities, I have localised 4 blocks - In each block, I have defined 4 plots in order to have the 4 possible treatments of both the water and mowing factors : nm/dry ; mowed/dry ; mowed/wet ; nm/wet. Here is my data table : Community Block Mowing WaterPH 1 Mesob1 Mowed Wet 7.40 2 Mesob1 nm Wet 13.10 3 Mesob1 Mowed Dry 5.55 4 Mesob1 nm Dry 10.35 5 Mesob2 nm Dry 10.70 6 Mesob2 Mowed Dry 6.38 7 Mesob2 nm Wet 9.75 8 Mesob2 Mowed Wet 6.35 9 Mesob3 nm Wet 9.60 10 Mesob3 Mowed Dry 5.10 11 Mesob3 nm Dry 10.05 12 Mesob3 Mowed Wet 6.25 13 Mesob4 nm Wet 9.00 14 Mesob4 Mowed Wet 6.50 15 Mesob4 nm Dry 7.75 16 Mesob4 Mowed Dry 5.90 17 Xerob5 nm Wet 7.69 18 Xerob5 Mowed Wet 8.11 19 Xerob5 nm Dry 3.98 20 Xerob5 Mowed Dry 3.69 21 Xerob6 nm Wet 5.24 22 Xerob6 Mowed Wet 4.22 23 Xerob6 nm Dry 6.55 24 Xerob6 Mowed Dry 4.40 25 Xerob7 Mowed Dry 3.79 26 Xerob7 nm Dry 3.91 27 Xerob7 nm Wet 9.00 28 Xerob7 Mowed Wet 8.50 29 Xerob8 Mowed Dry 3.33 30 Xerob8 nm Wet 6.25 31 Xerob8 Mowed Wet 8.00 32 Xerob8 nm Dry 6.33 I actually have 2 questions : I wrote my model in two different ways, and there were differences in P-Values according to the model written : First model : summary(aov(PH~Community*Mowing*Water + Error(Block))) Error: Block Df Sum Sq Mean Sq F value Pr(F) Community 1 42.182 42.182 24.407 0.002603 ** Residuals 6 10.370 1.728 --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 Error: Within Df Sum Sq Mean Sq F valuePr(F) Mowing 1 40.007 40.007 21.1747 0.0002215 *** Water 1 23.120 23.120 12.2370 0.0025673 ** Community:Mowing1 21.060 21.060 11.1467 0.0036554 ** Community:Water 1 6.901 6.901 3.6524 0.0720478 . Mowing:Water1 1.611 1.611 0.8527 0.3680090 Community:Mowing:Water 1 0.858 0.858 0.4542 0.5089331 Residuals 18 34.008 1.889 --- - Second model (assuming that Mowing*Water are nested inside the Block factor) : summary(aov(PH~Community*Mowing*Water + Error(Block/(Mowing*Water Error: Block Df Sum Sq Mean Sq F value Pr(F) Community 1 42.182 42.182 24.407 0.002603 ** Residuals 6 10.370 1.728 --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 Error: Block:Mowing Df Sum Sq Mean Sq F valuePr(F) Mowing1 40.007 40.007 37.791 0.0008489 *** Community:Mowing 1 21.060 21.060 19.893 0.0042820 ** Residuals 6 6.352 1.059 --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 Error: Block:Water Df Sum Sq Mean Sq F value Pr(F) Water1 23.1200 23.1200 6.0725 0.04884 * Community:Water 1 6.9006 6.9006 1.8125 0.22685 Residuals6 22.8439 3.8073 --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 Error: Block:Mowing:Water Df Sum Sq Mean Sq F value Pr(F) Mowing:Water1 1.6110 1.6110 2.0085 0.2062 Community:Mowing:Water 1 0.8581 0.8581 1.0697 0.3409 Residuals 6 4.8126 0.8021 Both models give me interesting (but different!) results. Which one would be the most appropriate? Second question : How can I verify preliminary assumptions (normality of residuals and variance homogeneity) in this kind of models? When I ask R to extract residuals, the answer is NULL: residuals(aov(PH~Community*Mowing*Water + Error(Block/(Mowing*Water NULL residuals(aov(PH~Community*Mowing*Water + Error(Block))) NULL A huge thanks to the one who will rescue (or at least try to rescue) my PhD! Sincerely, -- Jean-Paul Maalouf UMR 1202 BIOGECO Inra - Université Bordeaux 1 Bâtiment B8, Avenue des
Re: [R] Split plot analysis problems
Hi Jean-Paul, However, I've tried both solutions on my model, and I got different residuals :... What could be the difference between the two? There is no difference. You have made a mistake. ## tt - data.frame(read.csv(file=tt.csv, sep=)) ## imports your data set T.aov - aov(PH~Community*Mowing*Water + Error(Block), data=tt) summary(T.aov) ## This matches your output Df Sum Sq Mean Sq F valuePr(F) Mowing 1 40.007 40.007 21.1747 0.0002215 *** Water 1 23.120 23.120 12.2370 0.0025673 ** Community:Mowing1 21.060 21.060 11.1467 0.0036554 ** Community:Water 1 6.901 6.901 3.6524 0.0720478 . Mowing:Water1 1.611 1.611 0.8527 0.3680090 Community:Mowing:Water 1 0.858 0.858 0.4542 0.5089331 Residuals 18 34.008 1.889 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## length(residuals(T.aov - aov(PH~Community*Mowing*Water + Error(Block), data=tt)$Within)) [1] 24 length(proj(T.aov)[,Residuals]) [1] 24 ## Details residuals(T.aov - aov(PH~Community*Mowing*Water + Error(Block), data=tt)$Within) 9 10 11 12 13 14 15 16 -1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 0.215872118 -1.621627882 0.508372118 17 18 19 20 21 22 23 24 1.241311954 1.498811954 -0.616188046 0.483811954 -0.477827381 -1.660327381 2.684672619 1.924672619 25 26 27 28 29 30 31 32 -0.465198010 -1.735198010 1.502301990 0.839801990 -0.428515358 -0.751015358 0.836484642 1.181484642 proj(T.aov)[,Residuals] 9 10 11 12 13 14 15 16 17 18 19 20 21 -1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 0.215872118 -1.621627882 0.508372118 1.241311954 1.498811954 -0.616188046 0.483811954 -0.477827381 22 23 24 25 26 27 28 29 30 31 32 -1.660327381 2.684672619 1.924672619 -0.465198010 -1.735198010 1.502301990 0.839801990 -0.428515358 -0.751015358 0.836484642 1.181484642 Regards, Mark. Jean-Paul Maalouf wrote: Thanks Mark and Richard for your propositions on how to extract residuals. However, I've tried both solutions on my model, and I got different residuals : If we consider the within residuals : Mark's solution gave me a vector of 24 residuals : as.vector(residuals(aov(PH~Community*Mowing*Water + Error(Block))$Within)) [1] -1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 0.215872118 -1.621627882 0.508372118 1.241311954 1.498811954 -0.616188046 [12] 0.483811954 -0.477827381 -1.660327381 2.684672619 1.924672619 -0.465198010 -1.735198010 1.502301990 0.839801990 -0.428515358 -0.751015358 [23] 0.836484642 1.181484642 and Richard's solution gave me a vector 32 residuals : do.call(data.frame,proj(aov(PH~Community*Water*Mowing + Error(Block-proj proj$Within.Residuals [1] -0.216875 1.745625 -1.174375 -0.354375 0.800625 0.460625 -0.799375 -0.461875 -0.404375 -0.274375 0.695625 -0.016875 -0.541875 0.695625 -1.141875 [16] 0.988125 0.589375 0.846875 -1.268125 -0.168125 -1.095625 -2.278125 2.066875 1.306875 -0.500625 -1.770625 1.466875 0.804375 -0.638125 -0.960625 [31] 0.626875 0.971875 What could be the difference between the two? Regards JPM Quoting Richard M. Heiberger r...@temple.edu: Jean-Paul Maalouf wrote: Do you have any idea on how can I verify preliminary assumptions in this model (normality of the residuals and variance homogeneity), since R is not able to extract residuals? Of course, R extracts residuals. Use the proj() function. See ?proj for the example to get the projection of an aovlist object onto the terms of a linear model proj(npk.aovE) To get the projections into a simple data.frame, use tmpE - proj(npk.aovE) tmpE.df - do.call(data.frame, tmpE) tmpE.df Mark Difford's solution effectively retrieved columns 3 and 10 from tmpE.df Rich -- Jean-Paul Maalouf UMR 1202 BIOGECO Inra - Université Bordeaux 1 Bâtiment B8, Avenue des Facultés 33405 Talence, France Tel : 05 40008772 __ 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/Split-plot-analysis-problems-tp24589402p24632396.html Sent from the R
Re: [R] how to correlate nominal variables?
Hi Timo, I need functions to calculate Yule's Y or Cramérs Index... Are such functions existing? Also look at assocstats() in package vcd. Regards, Mark. Timo Stolz wrote: Dear R-Users, I need functions to calculate Yule's Y or Cramérs Index, in order to correlate variables that are nominally scaled? Am I wrong? Are such functions existing? Sincerely, Timo __ 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/how-to-correlate-nominal-variables--tp18441195p24688304.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] Adding picture to graph?
Hi Rainer, the question came up if it would be possible to add a picture (saved on the HDD) to a graph (generated by plot()), which we could not answer. Yes. Look at package pixmap and, especially, at the examples sub s.logo() in package ade4. Regards, Mark. Rainer M Krug-6 wrote: Hi while teaching R, the question came up if it would be possible to add a picture (saved on the HDD) to a graph (generated by plot()), which we could not answer. It might easily kill a clean graph, but: is there a way of doing this, even one should not do it? On a similar line of thought: is it possibe to define own symbols so that they can be used in the plot function with pch=? Rainer -- Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch University, South Africa __ 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/Adding-picture-to-graph--tp24714724p24716913.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] Selecting Bootstrap Method for Quantile Regression
Hi Tom, For example, if I want to use the xy-pair bootstrap how do I indicate this in summary.rq? The general approach is documented under summary.rq (sub se option 5). Shorter route is boot.rq, where examples are given. ## ?boot.rq y - rnorm(50) x - matrix(rnorm(100),50) fit - rq(y~x,tau = .4) summary(fit,se = boot, bsmethod= xy) Regards, Mark. Tom La Bone wrote: The help page and vignette for summary.rq(quantreg) mention that there are three different bootstrap methods available for the se=bootstrap argument, but I can't figure out how to select a particular method. For example, if I want to use the xy-pair bootstrap how do I indicate this in summary.rq? Tom -- View this message in context: http://www.nabble.com/Selecting-Bootstrap-Method-for-Quantile-Regression-tp24737299p24744179.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] i'm so stuck with text file and contour plot
Hannes, been trying to read a text file that contains heading in the first line in to R but cant. You want the following: ## TDat - read.csv(small.txt, sep=\t) TDat str(TDat) See ?read.csv Regards, Mark. hannesPretorius wrote: Ok i feel pretty stupid.. been trying to read a text file that contains heading in the first line in to R but cant. all i need to do is make a contour plot for a friend but after weeks i feel like giving up.. i included the first few lines of the file.. any help will be great Thanks Hannes http://www.nabble.com/file/p24777697/small.txt small.txt -- View this message in context: http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24780416.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] i'm so stuck with text file and contour plot
Hi David, I think he may also need to add the header=TRUE argument: No! The argument header= is not required in this case. ## TDat - read.csv(small.txt, sep=\t) str(TDat[,1:3]) 'data.frame': 10 obs. of 3 variables: $ Placename: Factor w/ 10 levels Aankoms,Aapieshoek,..: 1 2 3 4 5 6 7 8 9 10 $ X_coord : num 30.9 31.4 31.1 31.4 18.7 ... $ Y_coord : num -26.2 -27.4 -29 -29 -33.5 ... ?read.csv, sub header: If missing, the value is determined from the file format: header is set to TRUE if and only if the first row contains one fewer field than the number of columns. regards, Mark. David Winsemius wrote: I think he may also need to add the header=TRUE argument: tdat - read.csv(http://www.nabble.com/file/p24777697/small.txt;, header=TRUE, sep=\t) Note: read.table with those arguments should have worked as well. And then use names(tdat) - c(less bloated list of variable names) Perhaps along these lines: tdnames - names(tdat) tdnames #--don't paste-- [1] Placename [2] X_coord [3] Y_coord [4] Jan.to.Dec.2006.Stroke.Density.per.sq.km [5] Jan.to.Dec.2007.Stroke.Density.per.sq.km [6] Jan.to.Oct.2008.Stroke.Density.per.sq.km [7] Total.Strokes.per.sq.km.for.Jan.2006.to.Oct.2008 ###- names(tdat)[4:7] - c(Strk.dens.2006, Strk.dens.2007, Strk.dens. 2008, cumStrk.2006_8) # cannot use variable names that begin with numbers without special efforts tdat # now can be displayed more economically -- David On Aug 2, 2009, at 2:10 PM, Mark Difford wrote: Hannes, been trying to read a text file that contains heading in the first line in to R but cant. You want the following: ## TDat - read.csv(small.txt, sep=\t) TDat str(TDat) See ?read.csv Regards, Mark. hannesPretorius wrote: Ok i feel pretty stupid.. been trying to read a text file that contains heading in the first line in to R but cant. all i need to do is make a contour plot for a friend but after weeks i feel like giving up.. i included the first few lines of the file.. any help will be great Thanks Hannes http://www.nabble.com/file/p24777697/small.txt small.txt -- View this message in context: http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24780416.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. 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. -- View this message in context: http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24786217.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] i'm so stuck with text file and contour plot
And I meant to add, but somehow forgot, that the default for read.csv is header=TRUE (which is different from read.table, where it is FALSE). Regards, Mark. Mark Difford wrote: Hi David, I think he may also need to add the header=TRUE argument: No! The argument header= is not required in this case. ## TDat - read.csv(small.txt, sep=\t) str(TDat[,1:3]) 'data.frame': 10 obs. of 3 variables: $ Placename: Factor w/ 10 levels Aankoms,Aapieshoek,..: 1 2 3 4 5 6 7 8 9 10 $ X_coord : num 30.9 31.4 31.1 31.4 18.7 ... $ Y_coord : num -26.2 -27.4 -29 -29 -33.5 ... ?read.csv, sub header: If missing, the value is determined from the file format: header is set to TRUE if and only if the first row contains one fewer field than the number of columns. regards, Mark. David Winsemius wrote: I think he may also need to add the header=TRUE argument: tdat - read.csv(http://www.nabble.com/file/p24777697/small.txt;, header=TRUE, sep=\t) Note: read.table with those arguments should have worked as well. And then use names(tdat) - c(less bloated list of variable names) Perhaps along these lines: tdnames - names(tdat) tdnames #--don't paste-- [1] Placename [2] X_coord [3] Y_coord [4] Jan.to.Dec.2006.Stroke.Density.per.sq.km [5] Jan.to.Dec.2007.Stroke.Density.per.sq.km [6] Jan.to.Oct.2008.Stroke.Density.per.sq.km [7] Total.Strokes.per.sq.km.for.Jan.2006.to.Oct.2008 ###- names(tdat)[4:7] - c(Strk.dens.2006, Strk.dens.2007, Strk.dens. 2008, cumStrk.2006_8) # cannot use variable names that begin with numbers without special efforts tdat # now can be displayed more economically -- David On Aug 2, 2009, at 2:10 PM, Mark Difford wrote: Hannes, been trying to read a text file that contains heading in the first line in to R but cant. You want the following: ## TDat - read.csv(small.txt, sep=\t) TDat str(TDat) See ?read.csv Regards, Mark. hannesPretorius wrote: Ok i feel pretty stupid.. been trying to read a text file that contains heading in the first line in to R but cant. all i need to do is make a contour plot for a friend but after weeks i feel like giving up.. i included the first few lines of the file.. any help will be great Thanks Hannes http://www.nabble.com/file/p24777697/small.txt small.txt -- View this message in context: http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24780416.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. 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. -- View this message in context: http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24786283.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] i'm so stuck with text file and contour plot
Hannes,
When I read the entire text file in I get the following message
Then you have not followed the very simple instructions I gave you above,
which I repeat below. Or you have changed small.txt.
##
TDat - read.csv(small.txt, sep=\t)
TDat
str(TDat)
Mark.
hannesPretorius wrote:
When I read the entire text file in I get the following message
x - read.table('c:/small.txt', sep='\t', header=TRUE)
Warning message:
number of items read is not a multiple of the number of columns.
thanks.
hannesPretorius wrote:
Ok i feel pretty stupid.. been trying to read a text file that contains
heading in the first line in to R but cant. all i need to do is make a
contour plot for a friend but after weeks i feel like giving up.. i
included the first few lines of the file.. any help will be great
Thanks
Hannes http://www.nabble.com/file/p24777697/small.txt small.txt
--
View this message in context:
http://www.nabble.com/i%27m-so-stuck-with-text-file-and-contour-plot-tp24777697p24824996.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] Testing year effect in lm() ***failed first time, sending again
Emmanuel, somewhat incomplete help pages : what in h*ll are valid arguments to mcp() beyond Tukey ??? Curently, you'll have to dig in the source to learn that...). Not so: they are clearly stated in ?contrMat. Regards, Mark. Emmanuel Charpentier-3 wrote: Le jeudi 30 juillet 2009 à 16:41 -0600, Mark Na a écrit : Dear R-helpers, I have a linear model with a year effect (year is coded as a factor), i.e. the parameter estimates for each level of my year variable have significant P values (see some output below) and I am interested in testing: a) the overall effect of year; b) the significance of each year vis-a-vis every other year (the model output only tests each year against the baseline year). install.packges(multcomp) help.start() and use the vignettes ! They're good (and are a good complement to somewhat incomplete help pages : what in h*ll are valid arguments to mcp() beyond Tukey ??? Curently, you'll have to dig in the source to learn that...). HTH Emmanuel Charpentier __ 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/Testing-year-effect-in-lm%28%29-***failed-first-time%2C-sending-again-tp24748526p24832337.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] How do I get an inset graph (i.e. graph within a graph)?
Hi Arthur, This can be done quite easily using the appropriate arguments listed under ?par; and there are other approaches. Ready-made functions exist in several packages. I tend to use ?add.scatter from package ade4. It's a short function, so it's easy to customize it, but it works well straight out of the box. HTH, Mark. Arthur Roberts wrote: Hi, all, How do I get an inset graph (i.e. graph within a graph)? Your input is greatly appreciated. Best wishes, Art University of Washington Department of Medicinal Chemistry __ 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/How-do-I-get-an-inset-graph-%28i.e.-graph-within-a-graph%29--tp18796563p18798809.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] subset inside a lattice plot using panel.lines
Hi Michael,
Pulling my hair out here trying to get something very simple to work. ...
I can't quite see what you are trying to do [and I am not sure that you
clearly state it], but you could make things easier and simpler by (1)
creating a factor to identify your groups of rows more cleanly and (2) by
using the groups= argument in lattice.
## Something along these lines
copyDat - originalDat
copyDat$newFac - gl(2, 387, 774, labels=c(A,B))
xyplot(response ~ predictor|maybe.condition, data=copyDat, groups=newFac)
I hope this gives you some ideas.
Mark.
Michael Hopkins wrote:
Hi R people [duplicate - sorry, just posted HTML by mistake]
Pulling my hair out here trying to get something very simple to work.
Have a data frame of 774 rows and want to plot first and second half
on same axes with different colours. A variable is present call 'row'
created like this and checked to be OK:
row - seq( len=length( variable1 ) )
...so I just want to plot the two subsets where row = 387 and where
row 387. None of these work properly:
plots both over each other in correct colours
opt_plot2 - function( var_string, c, units ) {
print( xyplot( get(var_string) ~ variable1 | variable2, panel =
function( x, y, ... ) {
panel.xyplot( x, y, type =n, ... )
panel.grid( h=-1, v=-1, lty=3, lwd=1, col=grey )
panel.lines( x, y, col = red, subset = row = 387 )
panel.lines( x, y, col = dark green, subset = row 387 )
},
}
plots both just in red
... panel.lines( x[row = 387], y[row = 387], col = red )
panel.lines( x[row 387], y[row 387], col = dark
green )
first - (row = 387)
second - (row 387)
plots both over each other in correct colours
... panel.lines( x, y, col = red, subset = first )
panel.lines( x, y, col = dark green, subset = second )
plots both just in red
... panel.lines( x[first], y[first], col = red )
panel.lines( x[second], y[second], col = dark green )
I'm feeling frustrated and a bit stupid but should this be so
difficult? Any help or tips on what I am doing wrong greatly
appreciated.
TIA
Michael
__
Hopkins Research Touch the Future
__
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/subset-inside-a-lattice-plot-using-panel.lines-tp18813236p18813818.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] Is there any way to make pretty tables in R to pdf?
Hi Arthur, I was wondering if there was a package that can make pretty R tables to pdf. You got through TeX/LateX, but PDF could be your terminus. Package Hmisc: ? summary.formula and its various arguments and options. You can't get much better. http://cran.za.r-project.org/doc/contrib/Harrell-statcomp-notes.pdf HTH, Mark. Arthur Roberts wrote: Hi, all, All your comments have been very useful. I was wondering if there was a package that can make pretty R tables to pdf. I guess I could use xtable, but I would like something a little more elegant. Your input is greatly appreciated. Best wishes, Art __ 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/Is-there-any-way-to-make-pretty-tables-in-R-to-pdf--tp18818187p18818675.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] Is there any way to make pretty tables in R to pdf?
Hi Arthur, Sorry, sent you down the wrong track: this will help you to get there: http://biostat.mc.vanderbilt.edu/twiki/pub/Main/StatReport/summary.pdf Regards, Mark. Arthur Roberts wrote: Hi, all, All your comments have been very useful. I was wondering if there was a package that can make pretty R tables to pdf. I guess I could use xtable, but I would like something a little more elegant. Your input is greatly appreciated. Best wishes, Art __ 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/Is-there-any-way-to-make-pretty-tables-in-R-to-pdf--tp18818187p18818837.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] re quest for fine panel axis controls in lattice
Hi David, Specifically, within each panel, I want to set the limits for x and y equal to each other since it is paired data (using the max value of the two). In addition to the code Chuck Cleland sent you, you may want to square things up by adding the argument: aspect = iso before the closing bracket xyplot(..., aspect=iso) Regards, Mark. DavidChosid wrote: I'm trying to use fine axis controls in lattice for each panel. Specifically, within each panel, I want to set the limits for x and y equal to each other since it is paired data (using the max value of the two). Of course, I have no problems setting the limits for the entire plot but I am having trouble setting them for each specific panel. Could someone please provide me some guidance? Thanks in advance. David Chosid Massachusetts Division of Marine Fisheries 1213 Purchase St., 3rd Floor New Bedford, MA 02740 508-990-2860 x140 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. -- View this message in context: http://www.nabble.com/request-for-fine-panel-axis-controls-in-lattice-tp18830204p18831668.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] Correlation dichotomous factor, continous (numerical) and ordered factor
Hi Birgitle, ... my variables are dichotomous factors, continuous (numerical) and ordered factors. ... Now I am confused what I should use to calculate the correlation using all my variables and how I could do that in R. Professor Fox's package polycor will do this for you in a very nice way. Regards, Mark. Birgitle wrote: Hello R-User! I appologise in advance if this should also go into statistics but I am presently puzzled. I have a data.frame (about 300 rows and about 80 variables) and my variables are dichotomous factors, continuous (numerical) and ordered factors. I would like to calculate the linear correlation between every pair of my variables, because I would like to perform a logistic regression (glm()) without the correlation between variables. I thought I could use for the continous (numerical) and ordered factor a spearman correlation that is using the ranks. But I thought also that I have to use a contingency table for the dichotomous factors. I read also that it is possible to use a point-biserial correlation to calculate the correlation between dichotomous and continuous variables. Now I am confused what I should use to calculate the correlation using all my variables and how I could do that in R. Is it possible with cor(), rcorr(), cormat() or other R-functions using one of the available correlation-coefficients. I would be very happy if somebody could enlighten my darkness. Many thanks in advance. B. -- View this message in context: http://www.nabble.com/Correlation-dichotomous-factor%2C-continous-%28numerical%29-and-ordered-factor-tp18852158p18852399.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] defining the distance between axis label and axis
Hi Jörg, I haven't found anything in par()... No? Well don't bet your bottom $ on it (almost never in R). ?par (sub mgp). Mark. Jörg Groß wrote: Hi, How can I make the distance between an axis-label and the axis bigger? I haven't found anything in par()... __ 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/defining-the-distance-between-axis-label-and-axis-tp18858861p18858935.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] Problems using hetcor (polycor)
Hi Birgitle, You need to get this right if someone is going to spend their time helping you. Your code doesn't work: You have specified more columns in colClasses than you have in the provided data set. TestPart-read.table(TestPart.txt, header=TRUE,row.names=1, na.strings=NA ,colClasses = Classe81) Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, : line 1 did not have 82 elements TestPart-read.table(TestPart.txt, header=TRUE,row.names=1, na.strings=NA) Warning message: In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, : number of items read is not a multiple of the number of columns length(names(TestPart)) [1] 72 length(Classe81) [1] 82 That will never work. HTH, Mark. Birgitle wrote: Sorry if this post should be long but I tried to give you a piece of my data to reproduce my error message using hetcor: Fehler in result$rho : $ operator is invalid for atomic vectors Zusätzlich: Warning messages: 1: In polychor(x, y, ML = ML, std.err = std.err) : 1 row with zero marginal removed 2: In polychor(x, y, ML = ML, std.err = std.err) : the table has fewer than 2 rows Error in result$rho : $ operator is invalid for atomic vectors Additional: Warning message: 1: In polychor(x, y, ML = ML, std.err = std.err) : 1 row with zero marginal removed 2: In polychor(x, y, ML = ML, std.err = std.err) : the table has fewer than 2 rows Use tab delimited data at the end of the post. Copy in Texteditor and save as TestPart.txt Then use the following code library(methods) setClass(of) setAs(character, of, function(from) as.ordered(from)) Classe81-cclasses - c(rep(factor, 64), rep(numeric,6), rep (of, 12)) TestPart-read.table(TestPart.txt, header=TRUE,row.names=1, na.strings=NA ,colClasses = Classe81) str(TestPart) library(polycor) TestPart.hetcor-hetcor(TestPart, use=complete.obs) this will produce the above mentioned error message that I can not interprete. Would be great if somebody could help me to solve the problem. Thanks B. 1 2 3 4 7 8 9 10 12 13 14 15 16 17 18 19 21 22 23 25 27 28 29 30 31 33 34 3536 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 58 59 60 6162 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 AX1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 25 5 9 1 8.5 2.5 3 5 2 2 3 3 1 1 1 2 1 2 BX1 1 0 0 1 0 0 1 NA NA NA 0 0 0 0 1 0 0 1 0 0 1 0 NA NA NA NA NA NA NA NA0 0 0 1 0 NA NA NA NA NA NA NA 0 0 0 0 1 1 0 0 0 1 1 NA NA 6 1 3.25 2.25 5 5 2 2 3 3 1 1 1 1 1 1 CX1 1 0 0 1 0 0 1 1 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 15 3.5 6 1 5.5 5.5 5 5 2 2 1 2 1 1 1 1 2 2 DX1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 50 17.5 7.5 2.5 8.5 5 5 5 2 2 2 3 1 1 1 1 3 3 EX1 0 1
Re: [R] Problems using hetcor (polycor)
Hi Birgitle, It seems to be failing on those columns that have just a single entry (i.e = 1, with the rest as 0; having just 1, an NA, and then 0s gets you through). And there are other reasons for failure (in the call to get a positive definite matrix). The main problem lies in the calculation of standard errors for some combinations. You can get a result for all columns by turning this off. You can get standard errors for all the columns up to 60 by omitting columns 12, 23, and 44. You need to work out the rest yourself by columning forwards till you get a failure, then drop that column from the index. There probably isn't enough information to calculate SEs for the columns that cause an error. ## This works with default settings but without columns 12, 23 and 44 hetcor(TestPart[,c(1:11,13:22,24:43,45)], pd=T, std.err=T, use=complete.obs) ## The first fails; the second works hetcor(TestPart[,c(1:11,13:22,24:43,45:60)], pd=T, std.err=F) hetcor(TestPart[,c(1:72)], pd=F, std.err=F) HTH, Mark. Birgitle wrote: Sorry if this post should be long but I tried to give you a piece of my data to reproduce my error message using hetcor: Fehler in result$rho : $ operator is invalid for atomic vectors Zusätzlich: Warning messages: 1: In polychor(x, y, ML = ML, std.err = std.err) : 1 row with zero marginal removed 2: In polychor(x, y, ML = ML, std.err = std.err) : the table has fewer than 2 rows Error in result$rho : $ operator is invalid for atomic vectors Additional: Warning message: 1: In polychor(x, y, ML = ML, std.err = std.err) : 1 row with zero marginal removed 2: In polychor(x, y, ML = ML, std.err = std.err) : the table has fewer than 2 rows Use tab delimited data at the end of the post. Copy in Texteditor and save as TestPart.txt Then use the following code library(methods) setClass(of) setAs(character, of, function(from) as.ordered(from)) Classe81-cclasses - c(rep(factor, 64), rep(numeric,6), rep (of, 12)) TestPart-read.table(TestPart.txt, header=TRUE,row.names=1, na.strings=NA ,colClasses = Classe81) str(TestPart) library(polycor) TestPart.hetcor-hetcor(TestPart, use=complete.obs) this will produce the above mentioned error message that I can not interprete. Would be great if somebody could help me to solve the problem. Thanks B. 1 2 3 4 7 8 9 10 12 13 14 15 16 17 18 19 21 22 23 25 27 28 29 30 31 33 34 3536 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 58 59 60 6162 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 AX1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 25 5 9 1 8.5 2.5 3 5 2 2 3 3 1 1 1 2 1 2 BX1 1 0 0 1 0 0 1 NA NA NA 0 0 0 0 1 0 0 1 0 0 1 0 NA NA NA NA NA NA NA NA0 0 0 1 0 NA NA NA NA NA NA NA 0 0 0 0 1 1 0 0 0 1 1 NA NA 6 1 3.25 2.25 5 5 2 2 3 3 1 1 1 1 1 1 CX1 1 0 0 1 0 0 1 1 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 15 3.5 6 1 5.5 5.5 5 5 2 2 1 2 1 1 1 1 2 2 DX1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 1 0 0 0 0
Re: [R] Problems using hetcor (polycor)
Hi Birgitle, It seems than, that it is not possible to use all variables without somehow imputing missing values. It depends on what you are after. You can use the full data set if you set std.err=F and pd=F. Then exclude the columns that cause it to falter and redo with SEs turned on. You have the correlations; all you lack are SEs for ca. 5 columns. And I haven't tried your revised classification; that's for you to do. Regards, Mark. Birgitle wrote: Many Thanks Mark for your answer. It seems than, that it is not possible to use all variables without somehow imputing missing values. But I will try which variables I can finally use. Many thanks again. B. Mark Difford wrote: Hi Birgitle, It seems to be failing on those columns that have just a single entry (i.e = 1, with the rest as 0; having just 1, an NA, and then 0s gets you through). And there are other reasons for failure (in the call to get a positive definite matrix). The main problem lies in the calculation of standard errors for some combinations. You can get a result for all columns by turning this off. You can get standard errors for all the columns up to 60 by omitting columns 12, 23, and 44. You need to work out the rest yourself by columning forwards till you get a failure, then drop that column from the index. There probably isn't enough information to calculate SEs for the columns that cause an error. ## This works with default settings but without columns 12, 23 and 44 hetcor(TestPart[,c(1:11,13:22,24:43,45)], pd=T, std.err=T, use=complete.obs) ## The first fails; the second works hetcor(TestPart[,c(1:11,13:22,24:43,45:60)], pd=T, std.err=F) hetcor(TestPart[,c(1:72)], pd=F, std.err=F) HTH, Mark. -- View this message in context: http://www.nabble.com/Problems-using-hetcor-%28polycor%29-tp18867343p18870999.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] Where is the archive? Calling a function within a function.
Hi Kevin,
Where is the archive?
Start with this:
?RSiteSearch
HTH, Mark.
rkevinburton wrote:
I seem to remember this topic coming up before so I decided to look at the
archive and realized that I didn't know where it was. Is there a
searchable archive for this list? Thank you.
My question is calling a function from within a function. I have
smerge - function(d1,d2) {
temp - merge(d1,d2,all=TRUE,by.x=DayOfYear,by.y=DayOfYear)
return (xmerge(temp))
}
But I get an error message indicating that xmerge could not be found when
it is defined just above.
Thank you.
Kevin
__
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/Where-is-the-archive--Calling-a-function-within-a-function.-tp18874953p18875170.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] Tick marks that correspond with bars on barplot
Hi Megan, I would like to have an X-axis where the labels for the years line up after every two bars in the plot (there is one bar for hardwood, and another for softwood). It isn't clear to me from your description what you really want (I found no attachment)? What you seem to be trying to get is a tickmark/label at the centre/midpoint of each yearly grouping. If so, then look at the examples that the author of barplot() has provided, and _in particular_ read what the author of the method says under Value. You use the colMeans of the object returned by barplot() to position labels. Regards, Mark. Megan J Bellamy wrote: Hello all, I have created a barplot that shows change in hardwood/softwood density from 1965 to 2005 in 5 year periods (1965,1970, etc). I would like to have an X-axis where the labels for the years line up after every two bars in the plot (there is one bar for hardwood, and another for softwood). Below is my script: density-read.table(F:\\Megan\\Vtest.csv, header=TRUE, sep=,) attach (density) barplot(DENSITY,YEAR, col=c(blue, green, green, blue, blue, green, blue, green, green, blue, green, blue, blue, green , green, blue, blue, green)) legend(1,85,c(Softwood, Hardwood), fill=c(blue, green), bty=n) title(ylab=Density of Softwood and Hardwood by percent (%)) title(xlab=Year of measurement) title(main=Change in Softwood and Hardwood Density between 1965-2005 for PSP #80) axis(2, at=NULL) I have tried using the tck, axis.lty functions with no luck and have also tried axis(1, at=YEAR) # but only the first year (1965) comes up. Attached is the csv file. Any help with this would be greatly appreciated. Many thanks in advance, Megan __ 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/Tick-marks-that-correspond-with-bars-on-barplot-tp18890762p18892131.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] [lme4]Coef output with binomial lmer
Hi Tom, 1|ass%in%pop%in%fam This is non-standard, but as you have found, it works. The correct translation is in fact 1|fam/pop/ass and not 1|ass/pop/fam as suggested by Harold Doran. Dropping %, ass%in%pop%in%fam reads [means] as: nest ass in pop [= pop/ass], and then nest this in fam == fam/pop/ass HTH, Mark. T.C. Cameron wrote: Dear R users I have built the following model m1-lmer(y~harn+foodn+(1|ass%in%pop%in%fam),family = quasibinomial) where y-cbind(alive,dead) where harn and foodn are categorical factors and the random effect is a nested term to represent experimental structure e.g. Day/Block/Replicate ass= 5 level factor, pop= 2 populations per treatment factor in each assay, 7 reps per population The model can be family = quasibinomial or binomial My complete lack of understanding is in retrieving the coefficients for the fixed effects to back-transform the effects of my factors on proportional survival I get the following output: coef(m1) $`ass %in% pop %in% fam` (Intercept) harn1 harn2 foodn2 FALSE 1.0322375 -0.1939521 0.0310434 0.810084 TRUE0.5997679 -0.1939521 0.0310434 0.810084 Where FALSE and TRUE refer to some attribute of the random effect My hunch is that it refers to the Coefficients with (=TRUE) and without (=FALSE) the random effects? Any help appreciated Dr Tom C Cameron Genetics, Ecology and Evolution IICB, University of Leeds Leeds, UK Office: +44 (0)113 343 2837 Lab:+44 (0)113 343 2854 Fax:+44 (0)113 343 2835 Email: [EMAIL PROTECTED] Webpage: click here http://www.fbs.leeds.ac.uk/staff/profile.php?tag=Cameron_TC [[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/-lme4-Coef-output-with-binomial-lmer-tp18894407p18904468.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] Simple lme/lmer random effects questions
Hi Brandon, ...is it sufficient to leave the values as they are or should I generate unique names for all combinations of sleeve number and temperature, using something like data$sleeve.in.temp - factor(with(data, temp:sleeve)[drop=TRUE]) You might be luckier posting this on https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models You need to make a unique identifier for each of your five sets of sleeves (or you need to nest sleeve in temp) if you wish to take account of that structure (temp/sleeve) as a random effect. I could be mistaken, but I think there are subtle differences between ~ 1|temp/sleeve and ~ 1|data$sleeve.in.temp designs. confused on how to actually set up the random effects term for the models. Given my experimental setup, using the lme syntax... The difference between (1) random = ~ 1|sleeve, ... and (2) random = ~ 1+temp|sleeve is that the first will give a random intercept for each level of sleeve, but the slope is fixed to be the same as their is no random term for slope. In the second specification there is a random term for slope, viz temp, so both intercepts and slopes can --- and probably will --- vary. ## To get some idea of what's being done, look at the following example. ## random intercepts, parallel slopes mod1 - lme(distance ~ age, Orthodont, random = ~ 1 | Subject) ## random intercepts, separate slopes mod2 - lme(distance ~ age, Orthodont, random = ~ age | Subject) plot(augPred(mod1, primary=~age)) ## parallel slopes plot(augPred(mod2, primary=~age)) ## separate slopes HTH, Mark. Brandon Invergo wrote: Hello, I have two very rudimentary questions regarding the random effects terms in the lme and lmer functions. I apologize if this also partly strays into a general statistics question, but I'm a bit new to this all. So hopefully it'll be a quick problem to sort out... Here is my experimental setup: I raised butterflies in 5 different testing chambers all set to different temperatures. Within the testing chambers, the butterflies were held in 10 different sleeves, which were rotated daily to compensate for microenvironmental effects. I measured several traits of the butterflies and I am constructing models for each trait (unfortunately, multivariate analysis isn't possible). In my models, sex and temperature are fixed factors and the sleeve is a random effect. Most of the response variables are normally distributed, but there is one with a Gamma distribution (time until an event) and another with poisson distribution (counts), so some models use lme while others use lmer. I would like to determine if, despite the daily rotation, there are still random effects from the individual sleeves. My two questions (assuming I haven't already made grave errors in my description of the setup) are: 1) In my data file, the sleeve variable is just marked with a number 1 through 10; the temperature is noted in a different column, so the 50 sleeves do not have unique names, but rather there are 5 instances of each of the 10 sleeve numbers. If sleeve is to be properly included in the models as a random effect, is it sufficient to leave the values as they are or should I generate unique names for all combinations of sleeve number and temperature, using something like data$sleeve.in.temp - factor(with(data, temp:sleeve)[drop=TRUE]) 2) (this is the one that strays more into standard statistics territory, sorry) I'm a bit confused on how to actually set up the random effects term for the models. Given my experimental setup, using the lme syntax, should it be: model - lme(response ~ sex*temp, random=~temp|sleeve, data) or model - lme(response ~ sex*temp, random=~1|sleeve, data) or something else? I've searched and searched, but everything I find online seems to be significantly more advanced than what I'm doing, leaving me even more confused than when I started! Thank you very much for your help!! I want to be sure I do this analysis right Cheers, -brandon __ 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/Simple-lme-lmer-random-effects-questions-tp18933698p18939376.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] ylab with an exponent
what is the problem? A solution is: plot(1,2, ylab=expression(paste(insects , m^2))) The problem is very much more difficult to determine. stephen sefick wrote: plot(1,2, ylab= paste(insects, expression(m^2), sep= )) I get insects m^2 I would like m to the 2 what is the problem? -- 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. -- View this message in context: http://www.nabble.com/ylab-with-an-exponent-tp19003927p19004197.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] before-after control-impact analysis with R
Hi Nikolaos, My question again is: Why can't I reproduce the results? When I try a simple anova without any random factors: Lack of a right result probably has to do with the type of analysis of variance that is being done. The default in R is to use so-called Type I tests, for good reason. SAS, I think, still uses a type of test that many authorities consider to be meaningless, as a general, first-level, test. There is a long, long thread on this, going back to about approx April/May 1999, when a suitable Ripplyism was coined, which I believe found its way into the fortunes package. But RSiteSearch(type III) should do for a start. Also see ?drop1 library(car) ?Anova HTH, Mark. Nikolaos Lampadariou wrote: Hello everybody, In am trying to analyse a BACI experiment and I really want to do it with R (which I find really exciting). So, before moving on I though it would be a good idea to repeat some known experiments which are quite similar to my own. I tried to reproduce 2 published examples but without much success. The first one in particular is a published dataset analysed with SAS by McDonald et al. (2000). They also provide the original data as well as the SAS code. I don't know much about SAS and i really want to stick to R. So here follow 3 questions based on these 2 papers. Paper 1 (McDonald et al. 2000. Analysis of count data from before-after control-impact studies. J. Agric. Biol. Envir Stat., 5:262-279) Data were collected from 1984, 1989, 1990, 1991, 1993, and 1996. Samples from 1984 are considered as before an impact and the remaining years as after the impact. Each year, 96 transects were sampled (36 in the oiled area and 60 in the non-oiled area; 0 is for oiled and 1 for non-oiled). The authors compare 3 different ways of analysing the data including glmm. The data can be reproduced with the following commands (density is fake numbers but I can provide the original data since I've typed them in anyway): x-c(rep(0,36),rep(1,60)) oiled-rep(x,6) year-c(rep(1984,96), rep(1989,96), rep(1990,96), rep(1991,96), rep(1993,96), rep(1996,96)) density-runif(576, min=0, max=10) transect-c(rep(1:96,6)) oil-data.frame(oiled=oiled, transect=transect, year=year, density=density) Question 1: I can reproduce the results of the repeated measures anova with: oil.aov1-aov(density~factor(year)*factor(oiled)+Error(factor(transect)) But why is the following command not working? oil.aov2-aov(density~oiled*year + Error(oiled/transect), data=oil) After reading the R-help archive, as well as Chambers and Hasties (Statistical Models in S) and Pinheiro's and Bates (Mixed effects models in S and S-plus) I would expect that the correct model is the oil.aov2. As you might see from the data, transect is nested within oiled and I still get the wrong results when I try +Error(transect) or +Error(oiled:transect) and many other variants. Question 2 (on the same paper): The authors conclude that it is better to analyse the data with a Generalized Linear (Mixed) Model Technique. I tried lme and after reading Douglas Bates article (May 2005, vol. 5/1) also lmer as follows: oil.lme-lme(density~year*oiled, random=~1|oiled/transect) or harvest.lmer-lmer(abund~H*BA+BA/year+BA/year:H+site:BA+(1|H/site)) but again no luck. I will get always some error messages or some interactions missing. Paper 2 (Keough Quinn, 2000. Legislative vs. practical protection of an intertidal shoreline in southeastern Australia. Ecol. Appl. 10: 871-881) Data were collected from 1989, 1990, 1991, 1993, 1994, 1995, 1996, 1997). the 1989-1991 are the before impact years and the other 4 years are the after the impact years. Eight sites were samples (2 protected and 6 impacted). In this dataset, site is nested within harvest (H) and year is nested within before-after (BA). Also, site is considered as random by the authors. The data (fake again) can be produced with the following commands: site-c(rep(c(A1,A2, RR1, RR2, WT1, WT2, WT3, WT4),8)) H-c(rep(c(exp, exp, prot, pro, exp, exp, exp, exp), 8)) year-c(rep(1989,8), rep(1990,8), rep(1991,8), rep(1993,8), rep(1994,8), rep(1995,8), rep(1996,8), rep(1997,8)) BA-c(rep(bef,24), rep(after,40)) abund-runif(64, min=0, max=10) harvest-data.frame(abund, BA, H, site, year) Question 3. The authors use a complex anova design and here is part of their anova table which shows the design and the df. source.of.var df Harvesting(H) 1, 6 before-after(BA) 1, 6 H x BA1, 6 Site(H) 6, 31 Year(BA) 6, 31 Site x BA 6, 31 Year x H 6, 31 Resid.31 My question again is: Why can't I reproduce the results? When I try a simple anova without any random factors: harvest.lm-lm(abund~H*BA+H/site+BA/year+site:BA+year:H) I get close enought but the results are different (at least I think they are
Re: [R] before-after control-impact analysis with R
Hi ... Sorry, an e was erroneously elided from Ripley... Mark Difford wrote: Hi Nikolaos, My question again is: Why can't I reproduce the results? When I try a simple anova without any random factors: Lack of a right result probably has to do with the type of analysis of variance that is being done. The default in R is to use so-called Type I tests, for good reason. SAS, I think, still uses a type of test that many authorities consider to be meaningless, as a general, first-level, test. There is a long, long thread on this, going back to about approx April/May 1999, when a suitable Ripplyism was coined, which I believe found its way into the fortunes package. But RSiteSearch(type III) should do for a start. Also see ?drop1 library(car) ?Anova HTH, Mark. Nikolaos Lampadariou wrote: Hello everybody, In am trying to analyse a BACI experiment and I really want to do it with R (which I find really exciting). So, before moving on I though it would be a good idea to repeat some known experiments which are quite similar to my own. I tried to reproduce 2 published examples but without much success. The first one in particular is a published dataset analysed with SAS by McDonald et al. (2000). They also provide the original data as well as the SAS code. I don't know much about SAS and i really want to stick to R. So here follow 3 questions based on these 2 papers. Paper 1 (McDonald et al. 2000. Analysis of count data from before-after control-impact studies. J. Agric. Biol. Envir Stat., 5:262-279) Data were collected from 1984, 1989, 1990, 1991, 1993, and 1996. Samples from 1984 are considered as before an impact and the remaining years as after the impact. Each year, 96 transects were sampled (36 in the oiled area and 60 in the non-oiled area; 0 is for oiled and 1 for non-oiled). The authors compare 3 different ways of analysing the data including glmm. The data can be reproduced with the following commands (density is fake numbers but I can provide the original data since I've typed them in anyway): x-c(rep(0,36),rep(1,60)) oiled-rep(x,6) year-c(rep(1984,96), rep(1989,96), rep(1990,96), rep(1991,96), rep(1993,96), rep(1996,96)) density-runif(576, min=0, max=10) transect-c(rep(1:96,6)) oil-data.frame(oiled=oiled, transect=transect, year=year, density=density) Question 1: I can reproduce the results of the repeated measures anova with: oil.aov1-aov(density~factor(year)*factor(oiled)+Error(factor(transect)) But why is the following command not working? oil.aov2-aov(density~oiled*year + Error(oiled/transect), data=oil) After reading the R-help archive, as well as Chambers and Hasties (Statistical Models in S) and Pinheiro's and Bates (Mixed effects models in S and S-plus) I would expect that the correct model is the oil.aov2. As you might see from the data, transect is nested within oiled and I still get the wrong results when I try +Error(transect) or +Error(oiled:transect) and many other variants. Question 2 (on the same paper): The authors conclude that it is better to analyse the data with a Generalized Linear (Mixed) Model Technique. I tried lme and after reading Douglas Bates article (May 2005, vol. 5/1) also lmer as follows: oil.lme-lme(density~year*oiled, random=~1|oiled/transect) or harvest.lmer-lmer(abund~H*BA+BA/year+BA/year:H+site:BA+(1|H/site)) but again no luck. I will get always some error messages or some interactions missing. Paper 2 (Keough Quinn, 2000. Legislative vs. practical protection of an intertidal shoreline in southeastern Australia. Ecol. Appl. 10: 871-881) Data were collected from 1989, 1990, 1991, 1993, 1994, 1995, 1996, 1997). the 1989-1991 are the before impact years and the other 4 years are the after the impact years. Eight sites were samples (2 protected and 6 impacted). In this dataset, site is nested within harvest (H) and year is nested within before-after (BA). Also, site is considered as random by the authors. The data (fake again) can be produced with the following commands: site-c(rep(c(A1,A2, RR1, RR2, WT1, WT2, WT3, WT4),8)) H-c(rep(c(exp, exp, prot, pro, exp, exp, exp, exp), 8)) year-c(rep(1989,8), rep(1990,8), rep(1991,8), rep(1993,8), rep(1994,8), rep(1995,8), rep(1996,8), rep(1997,8)) BA-c(rep(bef,24), rep(after,40)) abund-runif(64, min=0, max=10) harvest-data.frame(abund, BA, H, site, year) Question 3. The authors use a complex anova design and here is part of their anova table which shows the design and the df. source.of.var df Harvesting(H) 1, 6 before-after(BA) 1, 6 H x BA1, 6 Site(H) 6, 31 Year(BA) 6, 31 Site x BA 6, 31 Year x H 6, 31 Resid.31 My question again is: Why can't I reproduce the results? When I try a simple anova without any random factors: harvest.lm-lm(abund~H*BA+H
Re: [R] Re lational Operators in Legend
Hi Lorenzo,
...but I would like to write that 5=k=15.
This is one way to do what you want
plot(1,1)
legend(topright, expression(paste(R[g]~k^{1/d[f]^{small}}~5=k, {}=15)))
HTH, Mark.
Lorenzo Isella wrote:
Dear All,
I am sure that what I am asking can be solved by less than a
one-liner, but I am having a hard time to get this right.
Inside a legend environment, I do not have any problem if I issue the
command:
expression(*R[g]*~* k^{1/d[f]^{small}}*, *k= {} *15)
but I would like to write that 5=k=15. It has to be trivial, but so
far attempts like:
expression(*R[g]*~* k^{1/d[f]^{small}}*, for 5 *={}* *k= {}
*15)
always return some error message and no improvement with all the
variations I tried out.
Does anyone know what I am doing wrong?
Many Thanks
Lorenzo
__
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/Relational-Operators-in-Legend-tp19032571p19033251.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] Re lational Operators in Legend
Hi Lorenzo,
I may (?) have left something out. It isn't clear what ~ is supposed to
mean; perhaps it is just a spacer, or perhaps you meant the following:
plot(1,1)
legend(topright, expression(paste(R[g] %~~% k^{1/d[f]^{small}},~5=k,
{}=15)))
HTH, Mark.
Mark Difford wrote:
Hi Lorenzo,
...but I would like to write that 5=k=15.
This is one way to do what you want
plot(1,1)
legend(topright, expression(paste(R[g]~k^{1/d[f]^{small}}~5=k,
{}=15)))
HTH, Mark.
Lorenzo Isella wrote:
Dear All,
I am sure that what I am asking can be solved by less than a
one-liner, but I am having a hard time to get this right.
Inside a legend environment, I do not have any problem if I issue the
command:
expression(*R[g]*~* k^{1/d[f]^{small}}*, *k= {} *15)
but I would like to write that 5=k=15. It has to be trivial, but so
far attempts like:
expression(*R[g]*~* k^{1/d[f]^{small}}*, for 5 *={}* *k= {}
*15)
always return some error message and no improvement with all the
variations I tried out.
Does anyone know what I am doing wrong?
Many Thanks
Lorenzo
__
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/Relational-Operators-in-Legend-tp19032571p19033652.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] How to solve this problem ?
Hi Daren, Small progress, ... m4 - list(m1=m1, m2=m2, m3=m3) boxplot(m4) It's always a good idea to have a look at your data first (assuming you haven't). This shows that the reliable instrument is m2. HTH, Mark. Daren Tan wrote: Small progress, I am relying on levene test to check for equality of variances. Is my understanding correct, the larger the p-value, the more likely the variances are the same ? trt [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 Levels: 1 2 3 4 levene.test(rep(rnorm(5), 4), trt, option=median) Modified Robust Brown-Forsythe Levene-type test based on the absolute deviations from the median data: rep(rnorm(5), 4) Test Statistic = 0, p-value = 1 From: [EMAIL PROTECTED] To: [EMAIL PROTECTED] Subject: How to solve this problem ? Date: Wed, 20 Aug 2008 13:33:23 + I have disabled html text editing mode, thanks to Prof. Ripley for the kind reminder. Given three geological devices that takes 5 readings at 4 environmental conditions (A to D). What will be the proper approach to select the most reliable device ? m1 - c(73,52,49,53,83,43,58,94,53,62,75,66,41,72,70,75,57,59,85,84) m2 - c(31,38,30,35,36,26,27,38,22,31,24,35,36,31,38,33,32,28,33,30) m3 - c(65,57,36,40,36,30,40,34,37,40,33,33,37,29,37,37,30,33,40,35) names(m1) - rep(LETTERS[1:4], each=5) names(m2) - rep(LETTERS[1:4], each=5) names(m3) - rep(LETTERS[1:4], each=5) Before writing this email, I have tried to compare the sd for each device at each condition, but ran into obstacle on how to formulate the null hypothesis. Alternative solution tried was ANOVA, I am unsure whether it can help, as it compares the differences in means of each group. Thanks _ Easily edit your photos like a pro with Photo Gallery. http://get.live.com/photogallery/overview __ 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/How-to-solve-this-problem---tp19069553p19084239.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] Test of Homogeneity of Variances
Have you read the documentation to either of the functions you are using?
?bartlett.test
Performs Bartlett's test of the null that the variances in each of the
groups (samples) are the same.
This explicitly tells you what is being tested, i.e. the null tested is that
var1 = var2.
?rnorm
Generates random deviates, so the answer [i.e. p-value] will (almost
certainly) differ each time. Only by setting a seed for the random number
generator to use will rnorm() generate the same number-set/distribution and
therefore the same p-value.
?set.seed
##
set.seed(7)
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
set.seed(7)
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
set(23)
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
HTH, Mark.
Daren Tan wrote:
I am testing whether the sample variances are equal. When p-value 0.05
(alpha), should accept null hypothesis (sample variances are equal) or
reject it ?
The two new examples with each having same sample variances also puzzle
me. Why are the p-values different ?
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
Bartlett test of homogeneity of variances
data: rep(rnorm(5), times = 4) and rep(1:5, each = 4)
Bartlett's K-squared = 0.8681, df = 4, p-value = 0.929
bartlett.test(rep(rnorm(5),times=4), rep(1:5, each=4))
Bartlett test of homogeneity of variances
data: rep(rnorm(5), times = 4) and rep(1:5, each = 4)
Bartlett's K-squared = 3.5599, df = 4, p-value = 0.4688
From: [EMAIL PROTECTED]
To: [EMAIL PROTECTED]; [EMAIL PROTECTED]
Date: Fri, 22 Aug 2008 11:25:36 -0400
Subject: RE: [R] Test of Homogeneity of Variances
What are your hypotheses? Once you state what they are, interpretation
should be straightforward.
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On Behalf Of Daren Tan
Sent: Friday, August 22, 2008 11:18 AM
To: [EMAIL PROTECTED]
Subject: [R] Test of Homogeneity of Variances
I am testing the homogeneity of variances via bartlett.test and
fligner.test. Using the following example, how should I interpret the
p-value in order to accept or reject the null hypothesis ?
set.seed(5)
x - rnorm(20)
bartlett.test(x, rep(1:5, each=4))
Bartlett test of homogeneity of variances
data: x and rep(1:5, each = 4)
Bartlett's K-squared = 1.7709, df = 4, p-value = 0.7778
fligner.test(x, rep(1:5, each=4))
Fligner-Killeen test of homogeneity of variances
data: x and rep(1:5, each = 4)
Fligner-Killeen:med chi-squared = 1.0819, df = 4, p-value = 0.8971
__
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 message, including any attachments, is for ...{{dropped:6}}
__
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/Test-of-Homogeneity-of-Variances-tp19109383p19112613.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] aov, lme, multcomp
Hi Richard, The tests give different Fs and ps. I know this comes up every once in a while on R-help so I did my homework. I see from these two threads: This is not so, or it is not necessarily so. The error structure of your two models is quite different, and this is (one reason) why the F- and p-values are different. For instance, try the following comparison: ## Example require(MASS) ## for oats data set require(nlme) ## for lme() require(multcomp) ## for multiple comparison stuff Aov.mod - aov(Y ~ N + V + Error(B/V), data = oats) Lme.mod - lme(Y ~ N + V, random = ~1 | B/V, data = oats) summary(Aov.mod) anova(Lme.mod) See: http://www.nabble.com/Tukey-HSD-(or-other-post-hoc-tests)-following-repeated-measures-ANOVA-td17508294.html#a17553029 The example itself is from MASS (Venables Ripley). HTH, Mark. Richard D. Morey wrote: I am doing an analysis and would like to use lme() and the multcomp package to do multiple comparisons. My design is a within subjects design with three crossed fixed factors (every participant sees every combination of three fixed factors A,B,C). Of course, I can use aov() to analyze this with an error term (leaving out the obvious bits): y ~ A*B*C+Error(Subject/(A*B*C)) I'd also like to use lme(), and so I use y ~ A*B*C, random= ~1|Subject The tests give different Fs and ps. I know this comes up every once in a while on R-help so I did my homework. I see from these two threads: http://www.biostat.wustl.edu/archives/html/s-news/2002-05/msg00095.html http://134.148.236.121/R/help/06/08/32763.html that this is the expected behavior because of the way grouping works with lme(). My questions are: 1. is this the correct random argument to lmer: anova(lme(Acc~A*B*C,random=list(Sub=pdBlocked(list( pdIdent(~1), pdIdent(~A-1), pdIdent(~B-1), pdIdent(~C-1,data=data)) 2. How much do the multiple comparisons depend on the random statement? 3. I'm also playing with lmer: Acc~A*B*C+(1|Sub) Is this the correct lmer call for the crossed factors? If not, can you point me towards the right one? 4. I'm not too concerned with getting correct Fs from the analyses (well, except for aov, where it is easy), I just want to make sure that I am fitting the same model to the data with all approaches, so that when I look at parameter estimates I know they are meaningful. Are the multiple comparisons I'll get out of lme and lmer meaningful with fully crossed factors, given that they are both tuned for nested factors? Thanks in advance. -- Richard D. Morey Assistant Professor Psychometrics and Statistics Rijksuniversiteit Groningen / University of Groningen __ 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/aov%2C-lme%2C-multcomp-tp19144362p19145003.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] Two envelopes problem
... To pick up on what Mark has said: it strikes me that this is related to the simplex, where the bounded nature of the vector space means that normal arithmetical operations (i.e. Euclidean) don't work---that is, they can be used, but the results are wrong. Covariances and correlations for instance, are wrong, something that Pearson noted more than a century ago. Taking logs of ratios of numbers (say a number divdided by geometric mean of the other numbers, as in Aitchison's set of transformations) in this space transfers everything to Euclidean space, so squaring the problem. This is why taking logs fixes things ?? Mark. statquant wrote: Duncan: I think I see what you're saying but the strange thing is that if you use the utility function log(x) rather than x, then the expected values are equal. Somehow, if you are correct and I think you are, then taking the log , fixes the distribution of x which is kind of odd to me. I'm sorry to belabor this non R related discussion and I won't say anything more about it but I worked/talked on this with someone for about a month a few years ago and we gave up so it's interesting for me to see this again. Mark -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Duncan Murdoch Sent: Tuesday, August 26, 2008 8:15 AM To: Jim Lemon Cc: r-help@r-project.org; Mario Subject: Re: [R] Two envelopes problem On 26/08/2008 7:54 AM, Jim Lemon wrote: Hi again, Oops, I meant the expected value of the swap is: 5*0.5 + 20*0.5 = 12.5 Too late, must get to bed. But that is still wrong. You want a conditional expectation, conditional on the observed value (10 in this case). The answer depends on the distribution of the amount X, where the envelopes contain X and 2X. For example, if you knew that X was at most 5, you would know you had just observed 2X, and switching would be a bad idea. The paradox arises because people want to put a nonsensical Unif(0, infinity) distribution on X. The Wikipedia article points out that it can also arise in cases where the distribution on X has infinite mean: a mathematically valid but still nonsensical possibility. 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-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/Two-envelopes-problem-tp19150357p19163195.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] random error with lme for repeated measures anova
Hi Jean-Pierre, A general comment is that I think you need to think more carefully about what you are trying to get out of your analysis. The random effects structure you are aiming for could be stretching your data a little thin. It might be a good idea to read through the archives of the R-sig-mixed-models mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models This should give you a better idea of some of the issues involved. If you can't source the source (PB) then there are other documents that might help. You could begin by locationing the directory where nlme package is installed and looking at the scripts in the scripts subdirectory. These are from PB. Baayen has the following draft document on his website. Though DRAFT'ed all over the place, it is well worth a read, even if you are not interested in linguistic analysis. I think it has now been published by Cambridge UP. http://www.ualberta.ca/~baayen/publications/baayenCUPstats.pdf Paul Bliese's document (Multilevel Modeling in R) also has some useful sections (find it under the contributed documents section on CRAN). HTH, Mark. Jean-Pierre Bresciani wrote: Hi, what is the appropriate syntax to get the random error correct when performing repeated measures anova with 'lme'. let's say i have 3 independent variables, with 'aov', i would write something like: aov(dep_var~(indep_var1*indep_var2*indep_var3) + Error(subject/(indep_var1*indep_var2*indep_var3)). With 'lme' however, i can't find the right formula. i tried things like: lme(dep_var~(indep_var1*indep_var2*indep_var3), random = ~1|subject) or nesting my independent variables in 'subject', but those are obviously wrong with my design. i'm quite clueless (and i haven't found any convincing piece of information about how to correctly use 'lme' or 'lmer'). So, any advice along that line is more than welcome. JP -- View this message in context: http://www.nabble.com/random-error-with-lme-for-repeated-measures-anova-tp19178239p19180027.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] convert princomp output to equation for plane?
Hi Bill, Since x, y,and z all have measurement errors attached, the proper way to do the fit is with principal components analysis, and to use the first component (called loadings in princomp output). The easiest way for you to do this is to use the pcr [principal component regression] function in the pls package. Be aware that unless you fit all components you will be carrying out a form of penalized regression. A small example follows (assumes that you have installed the pls package): ## lm.mod - lm(Ozone ~ Solar.R + Wind + Month, data=airquality) pc.mod - pcr(Ozone ~ Solar.R + Wind + Month, data=airquality) lm.mod coef(pc.mod, intercept = TRUE) coef(pc.mod, ncomp=1, intercept = TRUE) coef(pc.mod, ncomp=3, intercept = TRUE) Regards, Mark. William Simpson-2 wrote: I want to fit something like: z = b0 + b1*x + b2*y Since x, y,and z all have measurement errors attached, the proper way to do the fit is with principal components analysis, and to use the first component (called loadings in princomp output). My dumb question is: how do I convert the princomp output to equation coefficients in the format above? I guess another dumb question would be: how about getting the standard deviations of b0, b1, b2? Thanks very much for any help. Bill __ 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/convert-princomp-output-to-equation-for-plane--tp19182643p19197360.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] model II regression - how do I do it?
Hi Danilo, I need to do a model II linear regression, but I could not find out how!! The smatr package does so-called model II (major axis) regression. Regards, Mark. Danilo Muniz wrote: I need to do a model II linear regression, but I could not find out how!! I tryed to use the lm function, but I did not discovered how to specify the model (type I or type II) to the function... could you help me? -- Danilo Muniz [Gruingas Abdiel] [[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/model-II-regression---how-do-I-do-it--tp19224427p19225640.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] model II regression - how do I do it?
Hi Danilo, I need to do a model II linear regression, but I could not find out how!! The smatr package does so-called model II (major axis) regression. Regards, Mark. Danilo Muniz wrote: I need to do a model II linear regression, but I could not find out how!! I tryed to use the lm function, but I did not discovered how to specify the model (type I or type II) to the function... could you help me? -- Danilo Muniz [Gruingas Abdiel] [[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/model-II-regression---how-do-I-do-it--tp19224427p19225663.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] Density estimates in modelling framework
Hi Hadley, There is also locfit, which is very highly regarded by some authorities (e.g. Hastie, Tibs, and Friedman). Cheers, Mark. hadley wrote: Hi all, Do any packages implement density estimation in a modelling framework? I want to be able to do something like: dmodel - density(~ a + b, data = mydata) predict(dmodel, newdata) This isn't how sm or KernSmooth or base density estimation works. Are there other packages that do density estimation? Or is there some reason that this is a bad idea. Hadley -- http://had.co.nz/ __ 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/Density-estimates-in-modelling-framework-tp19225727p19226399.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] model II regression - how do I do it?
Hi Dylan, While this topic is fresh, are there any compelling reasons to use Model II regression? The fact that it is the type of regression used in principal component analysis makes it a compelling method. Compelling reason? It is used to take account of measurement errors in both y _and_ x. Usually practioners/analysts would sidestep this by putting what is measured with least error on the x-axis and regress y on that. For instance, in doing calibration curves in nutrient analysis where the analyte is measured using a spec. Then the spec reading goes on x. I can't tell you how often I have seen analysts put the (usually) inaccurately determined analyte on x and the spec reading on y. HTH, Mark. Dylan Beaudette-2 wrote: On Friday 29 August 2008, Mark Difford wrote: Hi Danilo, I need to do a model II linear regression, but I could not find out how!! The smatr package does so-called model II (major axis) regression. Regards, Mark. While this topic is fresh, are there any compelling reasons to use Model II regression? Cheers, Dylan -- Dylan Beaudette Soil Resource Laboratory http://casoilresource.lawr.ucdavis.edu/ University of California at Davis 530.754.7341 __ 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/model-II-regression---how-do-I-do-it--tp19224427p19226644.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] non-parametric Anova and tukeyHSD
Hi Stephen, See packages: coin nparcomp npmc There is also kruskalmc() in package pgirmess Regards, Mark. stephen sefick wrote: I have insect data from twelve sites and like most environmental data it is non-normal mostly. I would like to preform an anova and a means seperation like tukey's HSD in a nonparametric sense (on some sort of central tendency measure - median?). I am searching around at this time on the internet. Any suggestions, books, etc. would be greatly appreciated. -- Stephen Sefick Research Scientist Southeastern Natural Sciences Academy 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. -- View this message in context: http://www.nabble.com/non-parametric-Anova-and-tukeyHSD-tp19227579p19227986.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] ANCOVA/glm missing/ignored interaction combinations
Hi Lara, And I cant for the life of me work out why category one (semio1) is being ignored, missing etc. Nothing is being ignored Lara --- but you are ignoring the fact that your factors have been coded using the default contrasts in R, viz so-called treatment or Dunnett contrasts. That is, your intercept is semio1 and the other coefficients are calculated off it, hence the name treatment contrasts. In general this type of coding is recommended for GLM models (but this can raise deep divisions, as some don't consider them to be proper contrasts because they are not orthogonal). Perhaps you should read up on how dummy variables are coded. To really understand this you will have to move away from ?contrasts as a source. It might also help you to do the following on a simpler model: ?dummy.coef dummy.coef(saved.glm/lm.model) This shows how it is worked out. HTH, Mark. lara harrup (IAH-P) wrote: Hi I am using R version 2.7.2. on a windows XP OS and have a question concerning an analysis of covariance with count data I am trying to do, I will give details of a scaled down version of the analysis (as I have more covariates and need to take account of over-dispersion etc etc) but as I am sure it is only a simple problem but I just can't see how to fix it. I have a data set with count data as the response (total) and two continuous covariates and one categorical explanatory variable (semio). When I run the following lines, after loading the data and assigning 'semio' as a factor, taking into account that I want to consider two way interactions: model.a-glm(total~(temperature+humidity+semio)^2,family=poisson) summary(model.a) I get the output below. But not all interactions are listed there are 4 semio categories, 1,2,3 and 4 but only 2,3,and 4 are listed in the summary (semio2,semio3,semio4). And I cant for the life of me work out why category one (semio1) is being ignored, missing etc. Any help or suggestions would be most appreciated. Thanks in advance Lara [EMAIL PROTECTED] Call: glm(formula = total ~ (temperature + humidity + semio)^2, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -22.212 -5.132 -2.4843.200 18.793 Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) 23.848754 2.621291 9.098 2e-16 *** temperature -1.038284 0.150465 -6.901 5.18e-12 *** humidity -0.264912 0.033928 -7.808 5.81e-15 *** semio2 22.852664 1.291806 17.690 2e-16 *** semio33.699901 1.349007 2.743 0.0061 ** semio4 -1.851163 1.585997 -1.167 0.2431 temperature:humidity 0.012983 0.001983 6.545 5.94e-11 *** temperature:semio2 -0.870430 0.037602 -23.149 2e-16 *** temperature:semio3 -0.060846 0.038677 -1.573 0.1157 temperature:semio40.055288 0.046137 1.198 0.2308 humidity:semio2 -0.114718 0.013369 -8.581 2e-16 *** humidity:semio3 -0.031692 0.013794 -2.298 0.0216 * humidity:semio4 0.008425 0.016020 0.526 0.5990 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 10423.5 on 47 degrees of freedom Residual deviance: 2902.2 on 35 degrees of freedom AIC: 3086.4 Number of Fisher Scoring iterations: 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. -- View this message in context: http://www.nabble.com/ANCOVA-glm-missing-ignored-interaction-combinations-tp19285259p19286859.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] ANCOVA/glm missing/ignored interaction combinations
And perhaps I should also have added: fit your model without an intercept and look at your coefficients. You should be able to work it out from there quite easily. Anyway, you now have the main pieces. Regards, Mark. Mark Difford wrote: Hi Lara, And I cant for the life of me work out why category one (semio1) is being ignored, missing etc. Nothing is being ignored Lara --- but you are ignoring the fact that your factors have been coded using the default contrasts in R, viz so-called treatment or Dunnett contrasts. That is, your intercept is semio1 and the other coefficients are calculated off it, hence the name treatment contrasts. In general this type of coding is recommended for GLM models (but this can raise deep divisions, as some don't consider them to be proper contrasts because they are not orthogonal). Perhaps you should read up on how dummy variables are coded. To really understand this you will have to move away from ?contrasts as a source. It might also help you to do the following on a simpler model: ?dummy.coef dummy.coef(saved.glm/lm.model) This shows how it is worked out. HTH, Mark. lara harrup (IAH-P) wrote: Hi I am using R version 2.7.2. on a windows XP OS and have a question concerning an analysis of covariance with count data I am trying to do, I will give details of a scaled down version of the analysis (as I have more covariates and need to take account of over-dispersion etc etc) but as I am sure it is only a simple problem but I just can't see how to fix it. I have a data set with count data as the response (total) and two continuous covariates and one categorical explanatory variable (semio). When I run the following lines, after loading the data and assigning 'semio' as a factor, taking into account that I want to consider two way interactions: model.a-glm(total~(temperature+humidity+semio)^2,family=poisson) summary(model.a) I get the output below. But not all interactions are listed there are 4 semio categories, 1,2,3 and 4 but only 2,3,and 4 are listed in the summary (semio2,semio3,semio4). And I cant for the life of me work out why category one (semio1) is being ignored, missing etc. Any help or suggestions would be most appreciated. Thanks in advance Lara [EMAIL PROTECTED] Call: glm(formula = total ~ (temperature + humidity + semio)^2, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -22.212 -5.132 -2.4843.200 18.793 Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) 23.848754 2.621291 9.098 2e-16 *** temperature -1.038284 0.150465 -6.901 5.18e-12 *** humidity -0.264912 0.033928 -7.808 5.81e-15 *** semio2 22.852664 1.291806 17.690 2e-16 *** semio33.699901 1.349007 2.743 0.0061 ** semio4 -1.851163 1.585997 -1.167 0.2431 temperature:humidity 0.012983 0.001983 6.545 5.94e-11 *** temperature:semio2 -0.870430 0.037602 -23.149 2e-16 *** temperature:semio3 -0.060846 0.038677 -1.573 0.1157 temperature:semio40.055288 0.046137 1.198 0.2308 humidity:semio2 -0.114718 0.013369 -8.581 2e-16 *** humidity:semio3 -0.031692 0.013794 -2.298 0.0216 * humidity:semio4 0.008425 0.016020 0.526 0.5990 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 10423.5 on 47 degrees of freedom Residual deviance: 2902.2 on 35 degrees of freedom AIC: 3086.4 Number of Fisher Scoring iterations: 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. -- View this message in context: http://www.nabble.com/ANCOVA-glm-missing-ignored-interaction-combinations-tp19285259p19286930.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] Modality Test
Hi Amin, First, does R have a package that can implement the multimodality test, e.g., the Silverman test, DIP test, MAP test or Runt test. Jeremy Tantrum (a Ph.D. student of Werner Steutzle's, c. 2003/04) did some work on this. There is some useful code on Steutzle's website: http://www.stat.washington.edu/wxs/Stat593-s03/Code/jeremy-unimodality.R I used it last year when I was trying to solve the problem of how best to compare lots of density curves (age distributions of 3 spp. of tree euphorbias from about very different 35 sites). In particular I had to ensure that I wasn't creating spurious bimodality at a particular age range when combining sites. You might find it useful. Feel free to contact me off list if the code has gone, as I think I still have it (somewhere). Regards, Mark. Amin W. Mugera wrote: Dear Readers: I have two issues in nonparametric statistical analysis that i need help: First, does R have a package that can implement the multimodality test, e.g., the Silverman test, DIP test, MAP test or Runt test. I have seen an earlier thread (sometime in 2003) where someone was trying to write a code for the Silverman test of multimodality. Is there any other tests that can enable me to know how many modes are in a distribution? Second, i would like to test whether two distributions are equal. Does R have a package than can implement the Li (1996) test of the equality of two distributions? Is there any other test i can use rather than the Li test? Thank you in advance for your help. Amin Mugera Graduate Student AgEcon Dept. Kansas State 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/Modality-Test-tp19396085p19400095.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] Modality Test
Whoops! I think that should be Stuetzle --- though I very much doubt that he reads the list. Mark Difford wrote: Hi Amin, First, does R have a package that can implement the multimodality test, e.g., the Silverman test, DIP test, MAP test or Runt test. Jeremy Tantrum (a Ph.D. student of Werner Steutzle's, c. 2003/04) did some work on this. There is some useful code on Steutzle's website: http://www.stat.washington.edu/wxs/Stat593-s03/Code/jeremy-unimodality.R I used it last year when I was trying to solve the problem of how best to compare lots of density curves (age distributions of 3 spp. of tree euphorbias from about very different 35 sites). In particular I had to ensure that I wasn't creating spurious bimodality at a particular age range when combining sites. You might find it useful. Feel free to contact me off list if the code has gone, as I think I still have it (somewhere). Regards, Mark. Amin W. Mugera wrote: Dear Readers: I have two issues in nonparametric statistical analysis that i need help: First, does R have a package that can implement the multimodality test, e.g., the Silverman test, DIP test, MAP test or Runt test. I have seen an earlier thread (sometime in 2003) where someone was trying to write a code for the Silverman test of multimodality. Is there any other tests that can enable me to know how many modes are in a distribution? Second, i would like to test whether two distributions are equal. Does R have a package than can implement the Li (1996) test of the equality of two distributions? Is there any other test i can use rather than the Li test? Thank you in advance for your help. Amin Mugera Graduate Student AgEcon Dept. Kansas State 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/Modality-Test-tp19396085p19400138.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] Modality Test
Hi Amin, And I have just remembered that there is a function called curveRep in Frank Harrell's Hmisc package that might be useful, even if not quite in the channel of your enquiry. curveRep was added to the package after my struggles, so I never used it and so don't know how well it performs (quite well, I would think). Regards, Mark. Amin W. Mugera wrote: Dear Readers: I have two issues in nonparametric statistical analysis that i need help: First, does R have a package that can implement the multimodality test, e.g., the Silverman test, DIP test, MAP test or Runt test. I have seen an earlier thread (sometime in 2003) where someone was trying to write a code for the Silverman test of multimodality. Is there any other tests that can enable me to know how many modes are in a distribution? Second, i would like to test whether two distributions are equal. Does R have a package than can implement the Li (1996) test of the equality of two distributions? Is there any other test i can use rather than the Li test? Thank you in advance for your help. Amin Mugera Graduate Student AgEcon Dept. Kansas State 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/Modality-Test-tp19396085p19400426.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] Modality Test
Hi Bert,
## returns TRUE only if the distributions are the same ...
Thanks for the elegant code. Problem is, the result is elusive, constantly
slipping from view and then rolling into the dip. Perhaps it's broken in
this version of R (2.7.2.Patched). A fix would therefore be much
appreciated, as I may need to use it again in the near future. A little
example follows:
##
p1 - rnorm(100, 10, 5)
p2 - rnorm(100, 35, 5)
equaldist(p1, p2)
[1] FALSE
equaldist(c(p1,p2), c(p2,p1))
[1] FALSE
plot(density(c(p1, p2)), main=Twin Peaks)
lines(c(p1,p2))
?
Regards, Mark.
Bert Gunter wrote:
Here is a function that tests for equality of however many distributions
as
you like:
equaldist - function(...){
## ... numeric sample vectors from the possibly different distributions to
be tested
## returns TRUE only if the distributions are the same
FALSE
}
;-)
-- Bert Gunter
Genentech
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On
Behalf Of Mark Difford
Sent: Tuesday, September 09, 2008 1:23 PM
To: r-help@r-project.org
Subject: Re: [R] Modality Test
Hi Amin,
And I have just remembered that there is a function called curveRep in
Frank
Harrell's Hmisc package that might be useful, even if not quite in the
channel of your enquiry. curveRep was added to the package after my
struggles, so I never used it and so don't know how well it performs
(quite
well, I would think).
Regards, Mark.
Amin W. Mugera wrote:
Dear Readers:
I have two issues in nonparametric statistical analysis that i need
help:
First, does R have a package that can implement the multimodality test,
e.g., the Silverman test, DIP test, MAP test or Runt test. I have seen
an earlier thread (sometime in 2003) where someone was trying to write
a code for the Silverman test of multimodality. Is there any other
tests that can enable me to know how many modes are in a distribution?
Second, i would like to test whether two distributions are equal. Does R
have a package than can implement the Li (1996) test of the equality
of two distributions? Is there any other test i can use rather than the
Li test?
Thank you in advance for your help.
Amin Mugera
Graduate Student
AgEcon Dept. Kansas State 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/Modality-Test-tp19396085p19400426.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.
--
View this message in context:
http://www.nabble.com/Modality-Test-tp19396085p19412015.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] Greyed text in the background of a plot
Hi Agustin, Is there any way of having a greyed (ghosted) text... Yes! ## plot(1, type=n) text(1, Old Grey Whistle Test, col=slategray4, cex=2) text(1, y=1.2, OH!, col=grey95, cex=4) Then plot what you want on top. If you export or plot to a PDF/ps device the last-plotted items will overlie the first-plotted items. I think this is generally true. See: ?colors colors() ?text ?par Regards, Mark. Agustin Lobo-4 wrote: Hi! Is there any way of having a greyed (ghosted) text (i.e, 2006) in the background of a plot? I'm making a dynamic plot and would like to show the year of each time step as a big greyed text in the background. (the idea comes from Hans Rosling video: http://video.google.com/videoplay?docid=4237353244338529080sourceid=searchfeed ) Thanks Agus -- Dr. Agustin Lobo Institut de Ciencies de la Terra Jaume Almera (CSIC) LLuis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 email: [EMAIL PROTECTED] http://www.ija.csic.es/gt/obster __ 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/Greyed-text-in-the-background-of-a-plot-tp19456533p19457146.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] Greyed text in the background of a plot
Hi Yihui,
That's good, I like it! Very nice site.
Regards, Mark.
Yihui Xie wrote:
Well, his talk seems to have attracted a lot of people... You may
simply use gray text in your plot. Here is an example:
##
x = runif(10)
y = runif(10)
z = runif(10, 0.1, 0.3)
cl = rgb(runif(10), runif(10), runif(10), 0.5) # transparent colors!
par(mar = c(4, 4, 0.2, 0.2))
for (i in 1917:2007) {
x = x + rnorm(10, 0, 0.02)
y = y + rnorm(10, 0, 0.02)
z = abs(z + rnorm(10, 0, 0.05))
plot(x, y, xlim = c(0, 1), ylim = c(0, 1), type = n, panel.first = {
grid()
text(0.5, 0.5, i, cex = 5, col = gray) # here is the text!
})
symbols(x, y, circles = z, add = T, bg = cl, inches = 0.8)
box()
Sys.sleep(0.2)
}
##
Not difficult at all, right? :)
BTW, if you are interested in such animations, you may as well take a
look at my animation package:
http://cran.r-project.org/web/packages/animation/index.html
http://animation.yihui.name/
Regards,
Yihui
- Show quoted text -
On Fri, Sep 12, 2008 at 8:35 PM, Agustin Lobo [EMAIL PROTECTED] wrote:
Hi!
Is there any way of having a greyed (ghosted) text
(i.e, 2006) in the background of a plot?
I'm making a dynamic plot and would like to show the
year of each time step as a big greyed text in the background.
(the idea comes from Hans Rosling video:
http://video.google.com/videoplay?docid=4237353244338529080sourceid=searchfeed
)
Thanks
Agus
--
Dr. Agustin Lobo
Institut de Ciencies de la Terra Jaume Almera (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: [EMAIL PROTECTED]
http://www.ija.csic.es/gt/obster
--
Yihui Xie [EMAIL PROTECTED]
Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
Mobile: +86-15810805877
Homepage: http://www.yihui.name
School of Statistics, Room 1037, Mingde Main Building,
Renmin University of China, Beijing, 100872, China
__
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/Greyed-text-in-the-background-of-a-plot-tp19456533p19457541.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] Symbols on a capscale object plot
Hi Rodrigo,
I would like to use something like squares, triangles and circles (filled
and empty).
You would normally add this using points():
?points
##
plot(1:10, type=n)
points(1:5, pch=21:25, bg=1:5)
points(6:10, pch=21:25, bg=c(1,darkgrey,cyan, bisque))
points(6:10, y=rep(6,5), pch=1:5)
I don't know what vegan is doing at the moment but there used to be a points
command for doing this (sub ordiplot, I think). Vegan mostly still uses base
graphics so you can always get your scores [via: scores(ord.object): see
?scores] from your saved ordination object and add them as per above.
Regards, Mark.
Rodrigo Aluizio wrote:
Hi, I'm a beginner with R, but I'm getting excellent results with it.
Well I've got an capscale object (vegan package) and I want to made a
biplot
with symbols representing six groups of areas.
With the plot.cca function and some par attributes (like 'labels') I was
able to substitute the samples names for keyboard symbols
('x','o','#',ect),
but it's far from the ideal.
I've already search about it and find something using 'Hershey' fonts
symbols for a scatterplot, but it didn't work for the biplot (capscale
object).
I would like to use something like squares, triangles and circles (filled
and empty).
Does anyone have and idea to solve it?
Thank you for your attention and patience
Sorry if the English is not that good, I'm Brazilian.
Here is the script I'm using!
# The analysis
library(vegan)
library(xlsReadWrite)
PotiAbio-read.xls('PotiAbioCanoco.xls',sheet=1,rowNames=T)
PotiBio-read.xls('FatorialReplica.xls',sheet=8,rowNames=T)
attach(PotiAbio)
LogPotiBio-log(PotiBio+1)
dbRDA-capscale(t(LogPotiBio)~Environmental Variables,dist=bray,add=T)
dbRDA
#Preparing to generate and save the graphic
SymbolsRep-read.xls('Fatores.xls',sheet=2)
tiff('dbRDAPontos.tif',width=1250,height=1250,res=150)
plot.cca(dbRDA,type='none',display=c('bp','sites'))
text.cca(dbRDA,dis='cn',col=323232,cex=0.7,lwd=2,lty='dotted')
text.cca(dbRDA,dis='sites',col='black',cex=0.8,labels=FatoresSymbols$RepSimb)
dev.off()
___
MSc. Rodrigo Aluizio
Centro de Estudos do Mar/UFPR
Laboratório de Micropaleontologia
Avenida Beira Mar s/n - CEP 83255-000
Pontal do Paraná - PR - BRASIL
Fone: (0**41) 3455-1496 ramal 217
Fax: (0**41) 3455-1105
__
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/Symbols-on-a-capscale-object-plot-tp19469679p19470027.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] Symbols on a capscale object plot
Hi Rodrigo,
Maybe if I can define a factor that specifies these groups and use this
factor to assign the symbols.
Yes, that's the idea. Some ordination-oriented packages make this easier to
do than others. You should look through vegan carefully to see the full
offering. Nevertheless, it's pretty easy to do it from the ground up once
you get the knack.
You would do something like the following. Say you have a factor-column
called Site, the rows of which match the rows of the data frame (=mydf) you
used for the ordination. Then you could lasso your groups by doing:
points(ord.obj$scores[mydf$Site==MarkerA], pch=21, bg=red)
points(ord.obj$scores[1:10, ], pch=21, bg=red) ## by rownum()
HTH, Mark.
Rodrigo Aluizio wrote:
Thank You Mark,
With your tip now I'm able to chance the sites names for real symbols.
Well,
now I just have to find out how to change groups of sites names, that
compose a 'biofacies', with different symbols (actually six groups) in the
same plot.
Maybe if I can define a factor that specifies these groups and use this
factor to assign the symbols.
I don't know, I've to try, but I'm getting something, I hope.
Rodrigo.
--
From: Mark Difford [EMAIL PROTECTED]
Sent: Saturday, September 13, 2008 9:21 AM
To: r-help@r-project.org
Subject: Re: [R] Symbols on a capscale object plot
Hi Rodrigo,
I would like to use something like squares, triangles and circles
(filled
and empty).
You would normally add this using points():
?points
##
plot(1:10, type=n)
points(1:5, pch=21:25, bg=1:5)
points(6:10, pch=21:25, bg=c(1,darkgrey,cyan, bisque))
points(6:10, y=rep(6,5), pch=1:5)
I don't know what vegan is doing at the moment but there used to be a
points
command for doing this (sub ordiplot, I think). Vegan mostly still uses
base
graphics so you can always get your scores [via: scores(ord.object): see
?scores] from your saved ordination object and add them as per above.
Regards, Mark.
Rodrigo Aluizio wrote:
Hi, I'm a beginner with R, but I'm getting excellent results with it.
Well I've got an capscale object (vegan package) and I want to made a
biplot
with symbols representing six groups of areas.
With the plot.cca function and some par attributes (like 'labels') I
was
able to substitute the samples names for keyboard symbols
('x','o','#',ect),
but it's far from the ideal.
I've already search about it and find something using 'Hershey' fonts
symbols for a scatterplot, but it didn't work for the biplot (capscale
object).
I would like to use something like squares, triangles and circles
(filled
and empty).
Does anyone have and idea to solve it?
Thank you for your attention and patience
Sorry if the English is not that good, I'm Brazilian.
Here is the script I'm using!
# The analysis
library(vegan)
library(xlsReadWrite)
PotiAbio-read.xls('PotiAbioCanoco.xls',sheet=1,rowNames=T)
PotiBio-read.xls('FatorialReplica.xls',sheet=8,rowNames=T)
attach(PotiAbio)
LogPotiBio-log(PotiBio+1)
dbRDA-capscale(t(LogPotiBio)~Environmental
Variables,dist=bray,add=T)
dbRDA
#Preparing to generate and save the graphic
SymbolsRep-read.xls('Fatores.xls',sheet=2)
tiff('dbRDAPontos.tif',width=1250,height=1250,res=150)
plot.cca(dbRDA,type='none',display=c('bp','sites'))
text.cca(dbRDA,dis='cn',col=323232,cex=0.7,lwd=2,lty='dotted')
text.cca(dbRDA,dis='sites',col='black',cex=0.8,labels=FatoresSymbols$RepSimb)
dev.off()
___
MSc. Rodrigo Aluizio
Centro de Estudos do Mar/UFPR
Laboratório de Micropaleontologia
Avenida Beira Mar s/n - CEP 83255-000
Pontal do Paraná - PR - BRASIL
Fone: (0**41) 3455-1496 ramal 217
Fax: (0**41) 3455-1105
__
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/Symbols-on-a-capscale-object-plot-tp19469679p19470027.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.
--
View this message in context:
http://www.nabble.com/Symbols-on-a-capscale-object-plot-tp19469679p19473704.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r
Re: [R] Symbols on a capscale object plot
Hi Rogrido,
Sorry: The first points() call was missing a vital comma. It should have
been.
points(ord.obj$scores[mydf$Site==MarkerA, ], pch=21, bg=red)
See ?[
Mark Difford wrote:
Hi Rodrigo,
Maybe if I can define a factor that specifies these groups and use this
factor to assign the symbols.
Yes, that's the idea. Some ordination-oriented packages make this easier
to do than others. You should look through vegan carefully to see the full
offering. Nevertheless, it's pretty easy to do it from the ground up once
you get the knack.
You would do something like the following. Say you have a factor-column
called Site, the rows of which match the rows of the data frame (=mydf)
you used for the ordination. Then you could lasso your groups by doing:
points(ord.obj$scores[mydf$Site==MarkerA], pch=21, bg=red)
points(ord.obj$scores[1:10, ], pch=21, bg=red) ## by rownum()
HTH, Mark.
Rodrigo Aluizio wrote:
Thank You Mark,
With your tip now I'm able to chance the sites names for real symbols.
Well,
now I just have to find out how to change groups of sites names, that
compose a 'biofacies', with different symbols (actually six groups) in
the
same plot.
Maybe if I can define a factor that specifies these groups and use this
factor to assign the symbols.
I don't know, I've to try, but I'm getting something, I hope.
Rodrigo.
--
From: Mark Difford [EMAIL PROTECTED]
Sent: Saturday, September 13, 2008 9:21 AM
To: r-help@r-project.org
Subject: Re: [R] Symbols on a capscale object plot
Hi Rodrigo,
I would like to use something like squares, triangles and circles
(filled
and empty).
You would normally add this using points():
?points
##
plot(1:10, type=n)
points(1:5, pch=21:25, bg=1:5)
points(6:10, pch=21:25, bg=c(1,darkgrey,cyan, bisque))
points(6:10, y=rep(6,5), pch=1:5)
I don't know what vegan is doing at the moment but there used to be a
points
command for doing this (sub ordiplot, I think). Vegan mostly still uses
base
graphics so you can always get your scores [via: scores(ord.object):
see
?scores] from your saved ordination object and add them as per above.
Regards, Mark.
Rodrigo Aluizio wrote:
Hi, I'm a beginner with R, but I'm getting excellent results with it.
Well I've got an capscale object (vegan package) and I want to made a
biplot
with symbols representing six groups of areas.
With the plot.cca function and some par attributes (like 'labels') I
was
able to substitute the samples names for keyboard symbols
('x','o','#',ect),
but it's far from the ideal.
I've already search about it and find something using 'Hershey' fonts
symbols for a scatterplot, but it didn't work for the biplot (capscale
object).
I would like to use something like squares, triangles and circles
(filled
and empty).
Does anyone have and idea to solve it?
Thank you for your attention and patience
Sorry if the English is not that good, I'm Brazilian.
Here is the script I'm using!
# The analysis
library(vegan)
library(xlsReadWrite)
PotiAbio-read.xls('PotiAbioCanoco.xls',sheet=1,rowNames=T)
PotiBio-read.xls('FatorialReplica.xls',sheet=8,rowNames=T)
attach(PotiAbio)
LogPotiBio-log(PotiBio+1)
dbRDA-capscale(t(LogPotiBio)~Environmental
Variables,dist=bray,add=T)
dbRDA
#Preparing to generate and save the graphic
SymbolsRep-read.xls('Fatores.xls',sheet=2)
tiff('dbRDAPontos.tif',width=1250,height=1250,res=150)
plot.cca(dbRDA,type='none',display=c('bp','sites'))
text.cca(dbRDA,dis='cn',col=323232,cex=0.7,lwd=2,lty='dotted')
text.cca(dbRDA,dis='sites',col='black',cex=0.8,labels=FatoresSymbols$RepSimb)
dev.off()
___
MSc. Rodrigo Aluizio
Centro de Estudos do Mar/UFPR
Laboratório de Micropaleontologia
Avenida Beira Mar s/n - CEP 83255-000
Pontal do Paraná - PR - BRASIL
Fone: (0**41) 3455-1496 ramal 217
Fax: (0**41) 3455-1105
__
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/Symbols-on-a-capscale-object-plot-tp19469679p19470027.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.
--
View
Re: [R] Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?
Hi Roberto, but I can't figure out the /(Lobe*Tissue) part... This type of nesting is easier to do using lmer(). To do it using lme() you have to generate the crossed factor yourself. Do something like this: ## tfac - with(vslt, interaction(Lobe, Tissue, drop=T)) str(tfac); head(tfac) mod2-lme(Volume ~ Val*Lobe*Tissue, random = ~1|Subject/tfac, data = vslt) Pre-Scriptum: You can also use ?: but ?interaction is more flexible and powerful. Regards, Mark. roberto toro wrote: Hello, I'm using aov() to analyse changes in brain volume between males and females. For every subject (there are 331 in total) I have 8 volume measurements (4 different brain lobes and 2 different tissues (grey/white matter)). The data looks like this: Subject Sex LobeTissue Volume subect1 1 F g 262374 subect1 1 F w 173758 subect1 1 O g 67155 subect1 1 O w 30067 subect1 1 P g 117981 subect1 1 P w 85441 subect1 1 T g 185241 subect1 1 T w 83183 subect2 1 F g 255309 subect2 1 F w 164335 subect2 1 O g 71769 subect2 1 O w 31879 subect2 1 P g 120518 subect2 1 P w 90334 subect2 1 T g 168413 subect2 1 T w 75790 subect3 0 F g 243621 subect3 0 F w 167025 subect3 0 O g 65998 subect3 0 O w 29758 subect3 0 P g 118026 subect3 0 P w 91903 subect3 0 T g 156279 subect3 0 T w 82349 I'm trying to see if there is an interaction Sex*Lobe*Tissue. This is the command I use with aov(): mod1-aov(Volume~Sex*Lobe*Tissue+Error(Subject/(Lobe*Tissue)),data.vslt) Subject is a random effect, Sex, Lobe and Tissue are fixed effects; Sex is an outer factor (between subjects), and Lobe and Tissue are inner factors (within-subjects); and there is indeed a significant 3-way interaction. I was told, however, that the results reported by aov() may depend on the order of the factors (type I anova), and that is better to use lme() or lmer() with type II, but I'm struggling to find the right syntaxis... To begin, how should I write the model using lme() or lmer()?? I tried this with lme(): gvslt-groupedData(Volume~1|Subject,outer=~Val,inner=list(~Lobe,~Tissue),data=vslt) mod2-lme(Volume~Val*Lobe*Tissue,random=~1|Subject,data=gvslt) but I have interaction terms for every level of Lobe and Tissue, and 8 times the number of DF I should have... (around 331*8 instead of ~331). Using lmer(), the specification of Subject as a random effect is straightforward: mod2-lmer(Volume~Sex*Lobe*Tissue+(1|Subject),data.vslt) but I can't figure out the /(Lobe*Tissue) part... Thank you very much in advance! roberto __ 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/Help-please%21-How-to-code-a-mixed-model-with-2-within-subject-factors-using-lme-or-lmer--tp19479860p19480387.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] Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?
Hi Roberto, It's difficult to comment further on specifics without access to your data set. A general point is that the output from summary(aov.object) is not directly comparable with summary(lme.object). The latter gives you a summary of a fitted linear regression model, not an analysis of variance model, and what you see will depend on what contrasts were in place when the model was fitted. If you haven't changed these then they will be so-called treatment contrasts. What you are seeing for Lobe (which plainly is coded as a factor) in the output from summary(lme.object) are the regression coefficients for each level of Lobe relative to its reference/treatment/baseline level, which is your (Intercept). If you fitted your model with, say, Helmert or sum-to-zero contrasts then these values would change. To see what your current reference level is do levels(dataset$Lobe). See ?levels. What you want to look at to begin with is: anova(lme.object). HTH, Mark. roberto toro wrote: Thanks for answering Mark! I tried with the coding of the interaction you suggested: tfac-with(vlt,interaction(Lobe,Tissue,drop=T)) mod-lme(Volume~Sex*Lobe*Tissue,random=~1|Subject/tfac,data=vlt) But is it normal that the DF are 2303? DF is 2303 even for the estimate of LobeO that has only 662 values (331 for Tissue=white and 331 for Tissue=grey). I'm not sure either that Sex, Lobe and Tissue are correctly handled why are there different estimates called Sex:LobeO, Sex:LobeP, etc, and not just Sex:Lobe as with aov()?. Why there's Tissuew, but not Sex1, for example? Thanks again! roberto ps1. How would you code this with lmer()? ps2. this is part of the output of mod-lme: summary(mod) Linear mixed-effects model fit by REML Data: vlt AIC BIClogLik 57528.35 57639.98 -28745.17 Random effects: Formula: ~1 | Subject (Intercept) StdDev:11294.65 Formula: ~1 | tfac %in% Subject (Intercept) Residual StdDev:10569.03 4587.472 Fixed effects: Volume ~ Sex * Lobe * Tissue Value Std.Error DFt-value p-value (Intercept)245224.61 1511.124 2303 162.27963 0. Sex 2800.01 1866.312 3291.50029 0.1345 LobeO -180794.83 1526.084 2303 -118.46975 0. LobeP -131609.27 1526.084 2303 -86.23984 0. LobeT -73189.97 1526.084 2303 -47.95932 0. Tissuew-72461.05 1526.084 2303 -47.48168 0. Sex:LobeO-663.27 1884.789 2303 -0.35191 0.7249 Sex:LobeP -2146.08 1884.789 2303 -1.13863 0.2550 Sex:LobeT1379.49 1884.789 23030.73191 0.4643 Sex:Tissuew 5387.65 1884.789 23032.85849 0.0043 LobeO:Tissuew 43296.99 2158.209 2303 20.06154 0. LobeP:Tissuew 50952.21 2158.209 2303 23.60856 0. LobeT:Tissuew -15959.31 2158.209 2303 -7.39470 0. Sex:LobeO:Tissuew -5228.66 2665.494 2303 -1.96161 0.0499 Sex:LobeP:Tissuew -1482.83 2665.494 2303 -0.55631 0.5781 Sex:LobeT:Tissuew -6037.49 2665.494 2303 -2.26506 0.0236 __ 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/Help-please%21-How-to-code-a-mixed-model-with-2-within-subject-factors-using-lme-or-lmer--tp19480815p19481027.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] Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?
Hi Roberto, The other thing you can do --- if you don't wish to step across to lmer(), where you will be able to exactly replicate the crossed-factor error structure --- is stay with aov(... + Error()), but fit the factor you are interested in last. Assume it is Sex. Then fit your model as aov.model - aov(Volume ~ Lobe * Tissue * Sex + Error(Subject/(Lobe * Tissue)) This should give you a so-called Type-II test for Sex. You may verify this by fitting the model without the Error term and using Anova() from the car package (which does Type-II/III tests) to look at the SS and F values. I say should, because the only concern I have is whether this procedure is affected by the presence of an Error term in the model. Establishing this is beyond my capabilities. Regards, Mark. roberto toro wrote: Thanks for answering Mark! I tried with the coding of the interaction you suggested: tfac-with(vlt,interaction(Lobe,Tissue,drop=T)) mod-lme(Volume~Sex*Lobe*Tissue,random=~1|Subject/tfac,data=vlt) But is it normal that the DF are 2303? DF is 2303 even for the estimate of LobeO that has only 662 values (331 for Tissue=white and 331 for Tissue=grey). I'm not sure either that Sex, Lobe and Tissue are correctly handled why are there different estimates called Sex:LobeO, Sex:LobeP, etc, and not just Sex:Lobe as with aov()?. Why there's Tissuew, but not Sex1, for example? Thanks again! roberto ps1. How would you code this with lmer()? ps2. this is part of the output of mod-lme: summary(mod) Linear mixed-effects model fit by REML Data: vlt AIC BIClogLik 57528.35 57639.98 -28745.17 Random effects: Formula: ~1 | Subject (Intercept) StdDev:11294.65 Formula: ~1 | tfac %in% Subject (Intercept) Residual StdDev:10569.03 4587.472 Fixed effects: Volume ~ Sex * Lobe * Tissue Value Std.Error DFt-value p-value (Intercept)245224.61 1511.124 2303 162.27963 0. Sex 2800.01 1866.312 3291.50029 0.1345 LobeO -180794.83 1526.084 2303 -118.46975 0. LobeP -131609.27 1526.084 2303 -86.23984 0. LobeT -73189.97 1526.084 2303 -47.95932 0. Tissuew-72461.05 1526.084 2303 -47.48168 0. Sex:LobeO-663.27 1884.789 2303 -0.35191 0.7249 Sex:LobeP -2146.08 1884.789 2303 -1.13863 0.2550 Sex:LobeT1379.49 1884.789 23030.73191 0.4643 Sex:Tissuew 5387.65 1884.789 23032.85849 0.0043 LobeO:Tissuew 43296.99 2158.209 2303 20.06154 0. LobeP:Tissuew 50952.21 2158.209 2303 23.60856 0. LobeT:Tissuew -15959.31 2158.209 2303 -7.39470 0. Sex:LobeO:Tissuew -5228.66 2665.494 2303 -1.96161 0.0499 Sex:LobeP:Tissuew -1482.83 2665.494 2303 -0.55631 0.5781 Sex:LobeT:Tissuew -6037.49 2665.494 2303 -2.26506 0.0236 __ 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/Help-please%21-How-to-code-a-mixed-model-with-2-within-subject-factors-using-lme-or-lmer--tp19480815p19489323.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] Ward's Clustering Doubts
Hi Rodrigo, [apropos of Ward's method] ... we saw something like You must use it with Euclidean Distance... Strictly speaking this is probably correct, as Ward's method does an analysis of variance type of decomposition and so doesn't really make much sense (I think) unless Euclidean distance (i.e. least-squares) is used. However, there may be ways around this. First, because a distance metric is non-Euclidean does not mean that it is always non-Euclidean. You can test this using ?is.euclid in package ade4. You can also test your matrix by doing a principal co-ordinate analysis; then look for negative eigenvalues. If none are found, the matrix is Euclidean and it should be OK to use Ward's method on that data set. Probably a better approach is to make your distance matrix Euclidean. There are several functions in ade4 that will do this. The idea then is to present/compare the two solutions: the first using the uncorrected, non-Euclidean distance matrix, the second using the corrected version. You could use procrustes/co-inertia analysis to compare the two in an intermediate step. Regards, Mark. Rodrigo Aluizio wrote: Hi Everybody, Now I have a doubt that is more statistical than R's technical. Im working with ecology of recent Foraminifera. At the lab we used to perform cluster analysis using 1-Pearsons R and Wards method (we already saw it in bibliography of the area) which renders good results with our biological data. Recently, using R Software (vegan and Cluster packages) which allows the combination of any kind of distances matrix with any clustering method, we tried to used Bray Curtis + Wards (which seem to be more appropriate to a matrix with a lot of zeros) and it renders a better result. Furthermore, the results agree with our hypothesis and with the results we have got with the Distance-based Redundancy Analysis - dbRDA or CAP. It means, the analysis (Q-mode) clusters the stations according to the main physical, sedimentary and biological characteristics of the study area. We received some critical comments noticing that Wards Method accepts Euclidean Distance only. So, we made the analysis again using Euclidean Distance but we dont get the better results we had using 1-Pearsons R + Wards or Bray Curtis + Wards (actually any other distance + method combination rendered better results). Trying to find answers in the specialized literature we just got little more confused because in any moment we saw something like You must use it with Euclidean Distance and like I said above we already saw in some articles from respected journals, other kind of distance associated with the Ward's Clustering method. Is it wrong or is it non sense to do the analysis in the way we were doing? The results with Wards combined with 1-Pearsons R or Bray Curtis fit better with our hypothesis and have excellent agglomerative coefficients , but we dont want to make inappropriate statistical procedures. I'm starting to realize how powerful R is, but it doesn't justify doing nonsense statistics... I hope one of you may help us! Thank you in advance. Rodrigo. [[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/Ward%27s-Clustering-Doubts-tp19486028p19490991.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.
