Re: [R] Problems with lm()
In your data, subject is nested within sequence. Was that your intention? > a<-c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13,14,14) > b<-c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2) > c<-c(2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2) > d<-c(2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1) > e<-c(1739,1633,1481,1837,1780,2073,1374,1629,1555,1385,1756,1522,1566,1643, + 1939,1615,1475,1759,1388,1483,1127,1682,1542,1247,1235,1605,1598,1718 + ) > Data<-data.frame(subject=as.factor(a), + drug=as.factor(b), period=as.factor(c), + sequence=as.factor(d), Max=e) > Data subject drug period sequence Max 111 22 1739 212 12 1633 321 11 1481 422 21 1837 531 22 1780 632 12 2073 741 11 1374 842 21 1629 951 22 1555 10 52 12 1385 11 61 11 1756 12 62 21 1522 13 71 22 1566 14 72 12 1643 15 81 11 1939 16 82 21 1615 17 91 22 1475 18 92 12 1759 19 101 11 1388 20 102 21 1483 21 111 22 1127 22 112 12 1682 23 121 11 1542 24 122 21 1247 25 131 22 1235 26 132 12 1605 27 141 11 1598 28 142 21 1718 Andrew On Thu, Jun 19, 2008 at 04:29:16PM +0800, leeznar wrote: > Dear R-users: > > I am a new R-user and I have a question about lm > function. Here is my data. > a<-c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13,14,14) > b<-c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2) > c<-c(2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2) > d<-c(2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1) > e<-c(1739,1633,1481,1837,1780,2073,1374,1629,1555,1385,1756,1522,1566,1643,1939,1615,1475,1759,1388,1483,1127,1682,1542,1247,1235,1605,1598,1718 > ) > Data<-data.frame(subject=as.factor(a), > drug=as.factor(b), period=as.factor(c), > sequence=as.factor(d), Max=e) > > lm3<- lm(Max ~subject*sequence + sequence + period + > drug, data=Data) > print(lm3) > anova(lm3) > > When I use lm to fit the data, there are some problems > in ??subject*sequence??. I have use GLM in SPSS to > fit the same data, and it seems there is no problem. > > I don??t know where my problem is. How can I get the > same result with SPSS? How can I do? > > Best regards, > Hsin-Ya Lee > > > > > > __ > [[elided Yahoo spam]] > Content-Type: application/msword; name="Result_SPSS.doc" > Content-Transfer-Encoding: base64 > Content-Description: 3367377201-Result_SPSS.doc > Content-Disposition: attachment; filename="Result_SPSS.doc" > > > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] problem with NA and if
Hi Keld
you should read ?sum.
sum(c(1,2,NA), na.rm=TRUE)
Cheers
Andrew
On Fri, Jul 04, 2008 at 08:29:34AM +0200, Keld J?rn Simonsen wrote:
> Hi
>
> I would like to sum a number of time series, some of them having NA's
>
> Standard action is here that if I sum a value with a NA, then the result
> is NA. I would like it to just keep the value.
>
> I then try to:
>
> a = NA; if (a == NA) { a = 0}
>
> just to try it out, but it says
>
> Error in if (a == NA) { : missing value where TRUE/FALSE needed
>
> What is wrong, and can I do it smarter? I looked at na.action but I
> don't see how it affects addition of vectors, nor time series.
>
> Best regards
> keld
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Decomposing tests of interaction terms in mixed-effects models
Dear R colleagues,
a friend and I are trying to develop a modest workflow for the problem
of decomposing tests of higher-order terms into interpretable sets of
tests of lower order terms with conditioning.
For example, if the interaction between A (3 levels) and C (2 levels)
is significant, it may be of interest to ask whether or not A is
significant at level 1 of C and level 2 of C.
The textbook response seems to be to subset the data and perform the
tests on the two subsets, but using the RSS and DF from the original
fit.
We're trying to answer the question using new explanatory variables.
This approach (seems to) work just fine in aov, but not lme.
For example,
##
# Build the example dataset
set.seed(101)
Block <- gl(6, 6, 36, lab = paste("B", 1:6, sep = ""))
A <- gl(3, 2, 36, lab = paste("A", 1:3, sep = ""))
C <- gl(2, 1, 36, lab = paste("C", 1:2, sep = ""))
example <- data.frame(Block, A, C)
example$Y <- rnorm(36)*2 + as.numeric(example$A)*as.numeric(example$C)^2 +
3 * rep(rnorm(6), each=6)
with(example, interaction.plot(A, C, Y))
# lme
require(nlme)
anova(lme(Y~A*C, random = ~1|Block, data = example))
summary(aov(Y ~ Error(Block) + A*C, data = example))
# Significant 2-way interaction. Now we would like to test the effect
# of A at C1 and the effect of A at C2. Construct two new variables
# that decompose the interaction, so an LRT is possible.
example$AC <- example$AC1 <- example$AC2 <- interaction(example$A, example$C)
example$AC1[example$C == "C1"] <- "A1.C1" # A is constant at C1
example$AC2[example$C == "C2"] <- "A1.C2" # A is constant at C2
# lme fails (as does lmer)
lme(Y ~ AC1 + AC, random = ~ 1 | Block, data = example)
# but aov seems just fine.
summary(aov(Y ~ Error(Block) + AC1 + AC, data = example))
## AC was not significant, so A is not significant at C1
summary(aov(Y ~ Error(Block) + AC2 + AC, data = example))
## AC was significant, so A is significant at C2
## Some experimentation suggests that lme doesn't like the 'partial
## confounding' approach that we are using, rather than the variables
## that we have constructed.
lme(Y ~ AC, random = ~ 1 | Block, data = example)
lme(Y ~ A + AC, random = ~ 1 | Block, data = example)
######
Are we doing anything obviously wrong? Is there another approach to
the goal that we are trying to achieve?
Many thanks,
Andrew
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Decomposing tests of interaction terms in mixed-effects models
On Mon, Aug 04, 2008 at 10:17:38AM +0200, Peter Dalgaard wrote:
> Andrew Robinson wrote:
> >Dear R colleagues,
> >
> >a friend and I are trying to develop a modest workflow for the problem
> >of decomposing tests of higher-order terms into interpretable sets of
> >tests of lower order terms with conditioning.
> >
> >For example, if the interaction between A (3 levels) and C (2 levels)
> >is significant, it may be of interest to ask whether or not A is
> >significant at level 1 of C and level 2 of C.
> >
> >The textbook response seems to be to subset the data and perform the
> >tests on the two subsets, but using the RSS and DF from the original
> >fit.
> >
> >We're trying to answer the question using new explanatory variables.
> >This approach (seems to) work just fine in aov, but not lme.
> >
> >For example,
> >
> >##
> >
> ># Build the example dataset
> >
> >set.seed(101)
> >
> >Block <- gl(6, 6, 36, lab = paste("B", 1:6, sep = ""))
> >A <- gl(3, 2, 36, lab = paste("A", 1:3, sep = ""))
> >C <- gl(2, 1, 36, lab = paste("C", 1:2, sep = ""))
> >example <- data.frame(Block, A, C)
> >example$Y <- rnorm(36)*2 + as.numeric(example$A)*as.numeric(example$C)^2 +
> >3 * rep(rnorm(6), each=6)
> >
> >with(example, interaction.plot(A, C, Y))
> >
> ># lme
> >
> >require(nlme)
> >anova(lme(Y~A*C, random = ~1|Block, data = example))
> >
> >summary(aov(Y ~ Error(Block) + A*C, data = example))
> >
> ># Significant 2-way interaction. Now we would like to test the effect
> ># of A at C1 and the effect of A at C2. Construct two new variables
> ># that decompose the interaction, so an LRT is possible.
> >
> >example$AC <- example$AC1 <- example$AC2 <- interaction(example$A,
> >example$C)
> >example$AC1[example$C == "C1"] <- "A1.C1" # A is constant at C1
> >example$AC2[example$C == "C2"] <- "A1.C2" # A is constant at C2
> >
> ># lme fails (as does lmer)
> >
> >lme(Y ~ AC1 + AC, random = ~ 1 | Block, data = example)
> >
> ># but aov seems just fine.
> >
> >summary(aov(Y ~ Error(Block) + AC1 + AC, data = example))
> >
> >## AC was not significant, so A is not significant at C1
> >
> >summary(aov(Y ~ Error(Block) + AC2 + AC, data = example))
> >
> >## AC was significant, so A is significant at C2
> >
> >## Some experimentation suggests that lme doesn't like the 'partial
> >## confounding' approach that we are using, rather than the variables
> >## that we have constructed.
> >
> >lme(Y ~ AC, random = ~ 1 | Block, data = example)
> >lme(Y ~ A + AC, random = ~ 1 | Block, data = example)
> >
> >##
> >
> >Are we doing anything obviously wrong? Is there another approach to
> >the goal that we are trying to achieve?
> >
> This looks like a generic issue with lme/lmer, in that they are not
> happy with rank deficiencies in the design matrix.
>
> Here's a "dirty" trick which you are not really supposed to know about
> because it is hidden inside the "stats" namespace:
>
> M <- model.matrix(~AC1+AC, data=example)
> M2 <- stats:::Thin.col(M)
> summary(lme(Y ~ M2 - 1, random = ~ 1 | Block, data = example)
>
> (Thin.col(), Thin.row(), and Rank() are support functions for
> anova.mlm() et al., but maybe we should document them and put them out
> in the open.)
That is a neat idea, thanks, Peter, but it doesn't quite fit the bill.
The summary provides t-tests but I am hoping to find F-tests,
otherwise I'm not sure how to efficiently test A (3 levels) at the two
levels of C.
The anova.lme function doesn't help, sadly:
> anova(lme(Y ~ M2 - 1, random = ~ 1 | Block, data = example))
numDF denDF F-value p-value
M2 625 23.0198 <.0001
so I'm still flummoxed!
Andrew
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Decomposing tests of interaction terms in mixed-effects models
On Mon, Aug 04, 2008 at 02:51:48PM +0200, Peter Dalgaard wrote: > Andrew Robinson wrote: > >On Mon, Aug 04, 2008 at 10:17:38AM +0200, Peter Dalgaard wrote: > > > >>Andrew Robinson wrote: > >> > > > >That is a neat idea, thanks, Peter, but it doesn't quite fit the bill. > >The summary provides t-tests but I am hoping to find F-tests, > >otherwise I'm not sure how to efficiently test A (3 levels) at the two > >levels of C. > > > >The anova.lme function doesn't help, sadly: > > > > > >>anova(lme(Y ~ M2 - 1, random = ~ 1 | Block, data = example)) > >> > > numDF denDF F-value p-value > >M2 625 23.0198 <.0001 > > > >so I'm still flummoxed! > > > >Andrew > > > You do have to peek into M2 to know that the test is for whether the > last two coefs are zero, but how about > > > M3 <- M2[,2:4] > > M4 <- M2[,5:6] > > anova(lme(Y ~ M3+M4, random = ~ 1 | Block, data = example)) >numDF denDF F-value p-value > (Intercept) 125 10.66186 0.0032 > M3 325 55.31464 <.0001 > M4 225 1.27591 0.2967 Marvelous, many thanks, Peter. > Also, check out estimable() in the gmodels package. Will do. Actually had done, but will do again. Cheers Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Comparison of demographics between 2 study samples
On Thu, Aug 14, 2008 at 07:46:41PM +1000, Jim Lemon wrote: > On Wed, 2008-08-13 at 19:14 -0700, Mark Home wrote: > > Dear All: > > > > I have a clinical study where I would like to compare the demographic > > information for 2 samples in a study. The demographics include both > > categorical and continuous variables. I would like to be able to say > > whether the demographics are significantly different or not. > > > > The majority of papers that I have read use multiple techniques to achieve > > this (e.g., t-test for the continuous variables and either Fischer exact or > > Chi-square for categorical). I wonder whether this might lead to spurious > > differences due to multiple significance tests. Is there a better way to > > do this? > > > Hi Mark, > Most of these comparisons are uncorrected, as the aim is to demonstrate > that the samples come from the same population. Therefore, you aren't > worried about making a Type I error, but ignoring a sampling difference > that might bias your results. > > Jim Hi Mark, just following up on Jim's point, if your goal is to demonstrate that the samples come from the same population then you probably should take a look at equivalence testing. Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 JP, I suggest that you read the literature cited in the help file, most particularly Pinheiro, J.C., and Bates, D.M. (2000) "Mixed-Effects Models in S and S-PLUS", Springer. MASS (Modern Applied Statistics in S) also has some useful things in it. Cheers Andrew On Wed, Aug 27, 2008 at 03:23:58AM -0700, 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-tp19178239p19178239.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. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Sweave and/or beamer issue
Hi Michael,
I think that beamer needs that frames that contain any verbatim text
or R output to be declared as fragile.
Try
\begin{frame}[fragile]
...
\end{frame}
On Sat, Sep 06, 2008 at 06:22:41PM -0400, Michael Kubovy wrote:
> Dear Friends,
>
> I not sure whether this is an Sweave or a beamer problem.
>
> The Rnw file:
>
> \documentclass[compress,smaller]{beamer}
> %\documentclass{article}
> %\usepackage{beamerarticle}
>
> \usepackage{Sweave}
>
> \title{Psychophysics II}
> \date{September 9, 2008}
>
> \begin{document}
>
> \frame{
> \begin{Schunk}
> \begin{Sinput}
> > ro <- 0.2
> > c <- seq(from = -3, to = 4, by = 0.1)
> > fn <- 1 - pnorm(c)
> > fo <- 1 - pnorm(c, mean = 1)
> > h <- fo + ro - ro * fo
> > f <- fn
> > plot(h ~ f, type = "l", asp = 1, xlim = c(0.01, 1.01), ylim =
> c(0.01, 1.01), las = 1, xlab = "p(yes | old)",
> + ylab = "p(yes | new)", main = "Dual-process model,
> p(recollection) = 0.2")
> \end{Sinput}
> \end{Schunk}
> \includegraphics{20080909test-model}
>
> }
>
> \end{document}
>
> The resulting LaTeX
> ***FAILS: *
>
> Runaway argument?
> > ro <- 0.2 > c <- seq(from = -3, to = 4, by = 0.1) > fn <- 1 - pnorm
> \ETC.
> ! Paragraph ended before [EMAIL PROTECTED] was complete.
>
> \par
> l.27 }
>
> ?
>
> **
> But when the Rnw file starts:
> **
> %\documentclass[compress,smaller]{beamer}
> \documentclass{article}
> \usepackage{beamerarticle}
>
> ***IT DOES NOT FAIL: *
>
> _
> Professor Michael Kubovy
> University of Virginia
> Department of Psychology
> USPS: P.O.Box 400400Charlottesville, VA 22904-4400
> Parcels:Room 102Gilmer Hall
> McCormick RoadCharlottesville, VA 22903
> Office:B011+1-434-982-4729
> Lab:B019+1-434-982-4751
> Fax:+1-434-982-4766
> WWW:http://www.people.virginia.edu/~mk9y/
>
>
>
>
> [[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.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] removing a word, the following space and the next word
Hi Bob, I recommend doing some background reading on regular expressions[1] and using gsub(). Cheers Andrew [1] http://en.wikipedia.org/wiki/Regular_expressions On Sat, Sep 20, 2008 at 04:24:36PM +1000, Bob Green wrote: > > >Hello, > > I am hoping for advice as to how I could remove all words immediately > following the words 'Mr' or 'Mr.' in a csv file. For example, the > following phrases are included in lines of text (along with other Mr) > that could be anywhere in the file: "Mr Jones ate lunch" and "Mr > Smith was tied". > > > I want to remove the words Jones and Smith (etc) leaving the other > text intact. > > Any suggestions are appreciated, > > regards > > > Bob > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Deleting multiple variables
Mike, how about M_UC <- M_UC[,-(myvars[1]:myvars[2])] ? Andrew On Mon, Sep 22, 2008 at 11:04:34PM +0100, Michael Pearmain wrote: > Hi All, > i have searched the web for a simple solution but have been unable to find > one. Can anyone recommend a neat way of deleting multiple variable? > I see, i need to use dataframe$VAR<-NULL to get rid of one variable, > In my situation i need to delete all vars between two points. > > I've used the 'which' function to find these out and have assigned to myvar > >myvars > [1] 2 17 > > but i can't figure out how i should apply this? > > Should i loop through the values? (Psydo code below?) > > for (x in c(myvars[1]:myvars[2])) > (M_UC$x<-NULL)) > > Any help gratful > > Mike > > [[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. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] perl expression question
Hi Mark,
do you mean the regex to get the portion of the address after the
final slash? Something like
gsub(".*/([^/]*$)", "\\1", stock, fixed=FALSE)
Cheers
Andrew
On Mon, Sep 22, 2008 at 07:29:25PM -0500, [EMAIL PROTECTED] wrote:
> If I have the string below. does someone know a regular expression to
> just get the "BLC.NYSE". I bought the O'Reilley
> book and read it when I can and I study the solutions on the list but
> I'm still not self sufficient with these things. Thanks.
>
>
> stock<-"/opt/limsrv/mark/research/equity/projects/testDL/stock_data/fhdb/US/BLC.NYSE"
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] way to check if the evaluation of an expression returns an error?
?try Cheers, Andrew On Fri, Dec 14, 2007 at 02:10:52PM -0800, Thomas Pujol wrote: > Is there a recommended or "good" way to check if the evaluation of an > expression returns an error? > > e.g. > var(NA) > I wish var(NA) to return NA or "err", or some other value, even a > text-string, but not an error message. > > > I am using a loop to load many samples of data and to perform analysis on > each sample. On occasion, the function I execute returns an error message. > (e.g. when all data is NA, etc). I wish for the loop to continue to > execute, and to have the function return an NA or NULL rather then an error > message. > > > > - > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > > -- > This message has been scanned for viruses and > dangerous content by MailScanner, and is > believed to be clean. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Rbind-ing a list into one item
Hi Arjun,
try do.call()
Cheers
Andrew
On Wed, Dec 26, 2007 at 04:23:51PM -0500, Kondamani, Arjun (GMI - NY Corporate
Bonds) wrote:
> Hi,
>
> I am doing the following:
>
> 1. I have a list of files.. Files1=list.files("some
> directory",pattern="some pattern")
> 2. I define a list as res=vector("list", length(files1))
> 3. I read all the files into this list: res=lapply(files1, read.csv)
>
> I now want to rowbind all the items in the list into one big mass (all
> files have same number of columns). I tried lapply(res, rbind) but that
> did not work. Suggestions?
>
>
> This message w/attachments (message) may be privileged, confidential or
> proprietary, and if you are not an intended recipient, please notify the
> sender, do not use or share it and delete it. Unless specifically indicated,
> this message is not an offer to sell or a solicitation of any investment
> products or other financial product or service, an official confirmation of
> any transaction, or an official statement of Merrill Lynch. Subject to
> applicable law, Merrill Lynch may monitor, review and retain e-communications
> (EC) traveling through its networks/systems. The laws of the country of each
> sender/recipient may impact the handling of EC, and EC may be archived,
> supervised and produced in countries other than the country in which you are
> located. This message cannot be guaranteed to be secure or error-free. This
> message is subject to terms available at the following link:
> http://www.ml.com/e-communications_terms/. By messaging with Merrill Lynch
> you consent to the foregoing.
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
> --
> This message has been scanned for viruses and
> dangerous content by MailScanner, and is
> believed to be clean.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] groupedData function not found
Hi Andrea, did you try require(nlme) ? Adding a package is not necessarily the same as loading it. Cheers Andrew On Thu, Dec 27, 2007 at 10:42:27PM +, andrea previtali wrote: > > Hello, > I'm trying to use the groupedData function and R is giving me the message: > Error: can not find function "groupedData" > The code is coming from a textbook so I think it should be correct: > > pigs<-data.frame(cbind(pig.time,pig.id,pig.wt)) > pig.growth<-groupedData(pig.wt~pig.time|pig.id,data=pigs) > > I have added the package "nlme" and to confirm that it was installed > correctly I requested the list of functions included in the package > (library(help=nlme)) and I do see groupData in the list. > I am using R 2.6.1 in Windows XP. > > I'll appreciate your help. > Thanks, > Andrea Previtali > Post-doc fellow > Dept. of Biology, > Univ. of Utah, SLC, UT. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/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 message has been scanned for viruses and > dangerous content by MailScanner, and is > believed to be clean. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] [Off Topic] searching for a quote
Dear community, I'm trying to track down a quote, but can't recall the source or the exact structure - not very helpful, I know - something along the lines that: 80% of [applied] statistics is linear regression ... ? Does this ring a bell for anyone? Thanks, Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Grouping data
Hi Kimmo, try cut() to create a factor with levels according to the range of values, and (among other options) table() to make the table. Cheers Andrew. On Wed, Jan 16, 2008 at 10:06:23PM +0200, K. Elo wrote: > Hi, > > I am quite new to R (but like it very much!), so please apologize if > this is a too simple question. > > I have a large data frame consisting of data from a survey. There is, > for example, information about age and education (a numeric value from > 1-9). Now I would like to extract the total amount of each type of > education within different age groups (e.g. from 18 to 25, from 25 to > 35 etc.). How could I achieve this? (I have been thinking about > using 'subset', but if there are better ideas they are welcome :) ) > > An example might clarify my point. Let's assume the following data: > # age edu > 1 25 2 > 2 33 5 > 3 22 3 > 4 19 1 > 5 21 3 > 6 30 4 > 7 32 4 > 8 31 1 > > What I want to have is: > > edu 18-25 25-35 ... > 1 1 1 > 2 1 0 > 3 2 0 > 4 0 2 > 5 0 1 > > Thanks in advance & kind regards, > Kimmo > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/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 message has been scanned for viruses and > dangerous content by MailScanner, and is > believed to be clean. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] aaMI
Hi Navish, did you run require(aaMI) ? Cheers Andrew On Thu, Jan 17, 2008 at 02:17:12AM -0800, navish wrote: > > hi > i am new to R language. I want to use aaMI package which calculates the > amino acid mutual interaction for a given protein sequence. I had installed > the package but when i run the program it gives me the error could not find > function "aaMI". can anyone tell me what might be the problem.. > -- > View this message in context: > http://www.nabble.com/aaMI-tp14915744p14915744.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. > > -- > This message has been scanned for viruses and > dangerous content by MailScanner, and is > believed to be clean. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] [OT] - standard errors for parameter estimates under ridge regression and lasso?
Dear R community, I'm curious to know how people go about estimating standard errors for parameter estimates after model selection by ridge regression and the lasso. Do you have any practical or theoretical advice? Warmly, Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 Sebpe, the analysis of the data that you describe could be a complex and lengthy process, in which decisions that you are confronted by are affected by previous decisions that you have made. I recommend obtaining the assistance of a statistician, preferably a local one whose door you can knock on. If you are unable to do so then I suggest that you borrow/buy a copy of the Pinheiro and Bates book, which documents lme() and its friends, and study it carefully, especially the worked examples. Good luck! Andrew On Fri, Dec 12, 2008 at 03:13:06PM +, Sebpe De Smedt wrote: > Hi, > > > I'm working on leaf characteristics of trees. Each tree is characterised by > about 10 leaf traits. > The trees were sampled at 9 different locations (about 20 to 30 > trees/location, NOT balanced), grouped in 3 different climatic zones > (Sahelian, Soudanian and Guinean) (NOT balanced). > Further, each tree is characterised by some degree of human pressure > (mutilation degree), in total 4 different degrees were defined (NOT > balanced). > > In the dryer zones, the trees are under a much higher human pressure than in > the more humid climatic zones, "zone" and "mutilation degree" are thus > strongly correlated. > > I want to know how "zones" (fixed effects, climate interests me) and > "locations" (nested in "zones", random effects, location doesn't interests > me) are influencing the leaf traits (say for example "SLA"). Further, also > human pressure is affecting leaf traits so I want to characterise the > influence of "mutilation degree" (fixed effects) on "SLA". > > I found some interesting information, but still, I am not be able to analyse > the data properly. I think I have to use the function lme() or lme(). > > Can anyone tell me which function and command I have to use? And how I can > produce an ANOVA table? > > > Thanks in advance, > Sebpe De Smedt > > [[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. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Functions in R like lincom and nlcom of Stata
estimable in the gmodels package provides point estimates, standard errors and confidence intervals for arbitrary linear combinations of model parameters. I don't know for non-linear combinations, though. Cheers Andrew On Sat, Dec 13, 2008 at 11:33:12PM -0600, Stas Kolenikov wrote: > Those commands provide point estimates, standard errors and confidence > intervals based on linear combination of parameters or > linearization/delta-method, respectively. R's contrasts appear to be > limited to a single factor and combinations that sum up to zero. > > I am too so used to this Stata's concept, I now think it's odd R does > not seem to have it readily identifiable in two-three search commands. > And I would not believe R does not have this functionality, it must be > hiding somewhere! :)) > > On 12/13/08, David Winsemius wrote: > > > > On Dec 12, 2008, at 11:14 AM, Marc Mar? Dell'Olmo wrote: > > > > > > > Hello all, > > > > > > Does anyone know if there exists any function in R that resembles the > > > "lincom" and "nlcom" of STATA?. These functions computes point > > > estimates, standard errors, significance levels, confidence intervals, > > > etc. for linear and non linear combinations of previous estimated > > > parameters. Down here you've got links to descriptions of the > > > functions of STATA > > > > > > nlcom: > > > http://www.stata.com/help.cgi?nlcom > > > lincom: > > > http://www.stata.com/help.cgi?lincom > > > > > > > I did not find a description of the mathematical operations that let me > > understand exactly what lincom is doing, but suspect that you should be > > looking at how R handles contrasts. The help pages reference ch 2 of > > "Statistical Models in S". The search at the console prompt would be: > > > > ?C > > ?contrasts > > ?se.contrast > > ?model.tables > > > > -- > Stas Kolenikov, also found at http://stas.kolenikov.name > Small print: I use this email account for mailing lists only. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] p(H0|data) for lm/lmer-objects R
Dear Leo, > Dear R-List, > > I am interested in the Bayesian view on parameter estimation for multilevel > models and ordinary regression models. You might find Gelman & Hill's recent book to be good reading, and there is a book in the Use-R series that focuses on using R to perform Bayesian analyses. > AFAIU traditional frequentist p-values they give information about > p(data_or_extreme|H0). AFAIU it further, p-values in the Fisherian > sense are also no alpha/type I errors and therefor give no > information about future replications. I don't think that the last comment is necessarily relevant nor is it necessarily true. > However, p(data_or_extreme|H0) is not really interesting for social science > research questions (psychology). Much more interesting is > p(H0|data). That's fine, but first you have to believe that the statement has meaning. > Is there a way or formula to calculate these probabilities of the H0 > (or another hypothesis) from lm-/lmer objects in R? See the books above. Note that in order to do so, you will need to nominate a prior distribution of some kind. > Yes I know that multi-level modeling as well as regression can be done in a > purely Bayesian way. However, I am not capable of Bayesian statistics, > therefor I ask that question. I am starting to learn it a little bit. No offense, but it sounds to me like you want to have the Bayesian omelette without breaking the Bayesian eggs ;). Certain kinds of multi-level models are mathematically identical to certain kinds of Empirical Bayes models, but that does not make them Bayesian (despite what some people say). I caution against your implied goal of obtaining Bayesian statistics without performing a Bayesian analysis. Good luck, Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Problem with subset() function?
Steven, check the class of the objects that you are creating. Cheers, Andrew On Wed, January 21, 2009 10:02 am, Steven McKinney wrote: > Hi all, > > Can anyone explain why the following use of > the subset() function produces a different > outcome than the use of the "[" extractor? > > The subset() function as used in > > density(subset(mydf, ht >= 150.0 & wt <= 150.0, select = c(age))) > > appears to me from documentation to be equivalent to > > density(mydf[mydf$ht >= 150.0 & mydf$wt <= 150.0, "age"]) > > (modulo exclusion of NAs) but use of the former yields an > error from density.default() (shown below). > > > Is this a bug in the subset() machinery? Or is it > a documentation issue for the subset() function > documentation or density() documentation? > > I'm seeing issues such as this with newcomers to R > who initially seem to prefer using subset() instead > of the bracket extractor. At this point these functions > are clearly not exchangeable. Should code be patched > so that they are, or documentation amended to show > when use of subset() is not appropriate? > >> ### Bug in subset()? > >> set.seed(123) >> mydf <- data.frame(ht = 150 + 10 * rnorm(100), > +wt = 150 + 10 * rnorm(100), > +age = sample(20:60, size = 100, replace = TRUE) > +) > > >> density(subset(mydf, ht >= 150.0 & wt <= 150.0, select = c(age))) > Error in density.default(subset(mydf, ht >= 150 & wt <= 150, select = > c(age))) : > argument 'x' must be numeric > > >> density(mydf[mydf$ht >= 150.0 & mydf$wt <= 150.0, "age"]) > > Call: > density.default(x = mydf[mydf$ht >= 150 & mydf$wt <= 150, "age"]) > > Data: mydf[mydf$ht >= 150 & mydf$wt <= 150, "age"] (29 obs.); Bandwidth > 'bw' = 5.816 > >xy > Min. : 4.553 Min. :3.781e-05 > 1st Qu.:22.776 1st Qu.:3.108e-03 > Median :41.000 Median :1.775e-02 > Mean :41.000 Mean :1.370e-02 > 3rd Qu.:59.224 3rd Qu.:2.128e-02 > Max. :77.447 Max. :2.665e-02 > > >> sessionInfo() > R version 2.8.0 Patched (2008-11-06 r46845) > powerpc-apple-darwin9.5.0 > > locale: > C > > attached base packages: > [1] stats graphics grDevices datasets utils methods base > > loaded via a namespace (and not attached): > [1] Matrix_0.999375-16 grid_2.8.0 lattice_0.17-15 > lme4_0.99875-9 > [5] nlme_3.1-89 >> > > > > > > > Steven McKinney > > Statistician > Molecular Oncology and Breast Cancer Program > British Columbia Cancer Research Centre > > email: smckinney +at+ bccrc +dot+ ca > > tel: 604-675-8000 x7561 > > BCCRC > Molecular Oncology > 675 West 10th Ave, Floor 4 > Vancouver B.C. > V5Z 1L3 > Canada > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > Andrew Robinson Senior Lecturer in Statistics Tel: +61-3-8344-6410 Department of Mathematics and StatisticsFax: +61-3-8344 4599 University of Melbourne, VIC 3010 Australia Email: a.robin...@ms.unimelb.edu.auWebsite: http://www.ms.unimelb.edu.au __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 incomplete lmer output?
Hi Maria, the output is not incomplete. See the FAQ here: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-are-p_002dvalues-not-displayed-when-using-lmer_0028_0029_003f and the lengthy discussions about p-values and lmer (try Google). Andrew On Tue, Apr 15, 2008 at 10:51:46AM +0200, Maria Meyen wrote: > Dear all, > > I have problems with the lmer output (see > below)(2.6.2, WindowsXP), it is incomplete (no DFs, no > p-values)and also generating an anova-table doesn?t > work out. What is the fault? Does anybody know? > > This example is about plant height within different > treatments (each with 5 replicates) measured several > times. I?m interested in the treatment effect. > > > pap<-groupedData(heig~meas|box,outer=~treat,pap) > > mod<-lmer(heig~treat+(meas|box),pap) > > summary(mod) > Linear mixed-effects model fit by REML > Formula: heig ~ treat + (meas | box) >Data: pap >AIC BIC logLik MLdeviance REMLdeviance > 134.1 156.6 -58.05 97.14116.1 > Random effects: > Groups NameVariance Std.Dev. Corr > box (Intercept) 0.621046 0.78806 > meas0.790570 0.88914 -1.000 > Residual 0.037141 0.19272 > number of obs: 90, groups: box, 30 > > Fixed effects: > Estimate Std. Error t value > (Intercept) 1.794730.08353 21.487 > treatB1 0.233340.11813 1.975 > treatB2 0.626770.11813 5.306 > treatD1 0.079110.11813 0.670 > treatD2 0.439420.11813 3.720 > treatU -0.119330.11813 -1.010 > > Correlation of Fixed Effects: > (Intr) tretB1 tretB2 tretD1 tretD2 > treatB1 -0.707 > treatB2 -0.707 0.500 > treatD1 -0.707 0.500 0.500 > treatD2 -0.707 0.500 0.500 0.500 > treatU -0.707 0.500 0.500 0.500 0.500 > > > anova(mod) > Analysis of Variance Table > Df Sum Sq Mean Sq > treat 5 2.11106 0.42221 > > Thanks a lot in advance > > > > Maria Meyen > Department of Conservation Biology > Faculty of Biology > Philipps-University of Marburg > Karl-von-Frisch-Str.8 > 35043 Marburg > > > > __ > > Der Lieblings-Mailbox der Welt. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Errors using nlme's gls with autocorrelation
You're running out of RAM. Your options are 1) run code on a machine with more RAM 2) try the model on fewer observations 3) try a simpler model. Andrew On Fri, May 02, 2008 at 12:37:15PM -0700, zerfetzen wrote: > > Hi, > I am trying out a generalized least squares method of forecasting that > corrects for autocorrelation. I downloaded daily stock data from Yahoo > Finance, and am trying to predict Close (n=7903). I have learned to use > date functions to extract indicator variables for Monday - Friday (and > Friday is missing in the model to prevent it from becoming full rank). When > I run the following code... > > > library(nlme) > > MyModel <- gls(Close ~ Monday + Tuesday + Wednesday + Thursday, > correlation=corARMA(p=2), data=MyData, method="ML") > > ...I get the following error... > > Error in corFactor.corARMA(object) : > Calloc could not allocate (62457409 of 8) memory > > ...Does anybody know what I'm doing wrong? I appreciate any help. Thanks. > -- > View this message in context: > http://www.nabble.com/Errors-using-nlme%27s-gls-with-autocorrelation-tp17026417p17026417.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. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] GLMM and data manipulation (2nd try)
On Fri, May 02, 2008 at 03:06:31PM -0500, Giovanni Petris wrote:
>
> Hello,
>
> I posted a question yesterday but I got no replies, so I'll try to
> reformulate it in a more concise way.
>
> I have the following data, summarizing approval ratings on two
> different surveys for a random sample of 1600 individuals:
>
> > ## Example: Ratings of prime minister (Agresti, Table 12.1, p.494)
> > rating <- matrix(c(794, 86, 150, 570), 2, 2)
> > dimnames(rating) <- list(First = c("approve", "disapprove"),
> + Second = c("approve", "disapprove"))
> > rating
> Second
> Firstapprove disapprove
> approve794150
> disapprove 86570
>
> I would like to fit a logit model with approve/disapprove as response,
> survey (first/second) as a fixed effect, and subject as a random
> effect.
>
> 1) Is it possible to fit such a model directly using "lmer"?
>
> or
>
> 2) Should I unroll the table above into a dataframe containing also
>fictitious subject id's? If this is the case, what is a clean way
>of doing it?
Unroll it.
Asking for a "clean" way to do something is a disincentive because it
implies that you know how to do it but not cleanly. In the future I
would suggest that you do one of two things
a) post your on dirty version and ask for a cleaner one, or
b) just ask for something that works.
Something like ...
ratings <-
data.frame(
response = c(
rep(c(1,1), 794),
rep(c(1,0), 150),
rep(c(0,1), 86),
rep(c(0,0), 570)),
time = rep(c(1,2), 1600),
subject=rep(1:1600, each=2))
test.lmer <- lmer(response ~ time + (1|subject), data=ratings,
family=binomial)
but I don't know if you think that's clean or not.
Andrew
> Thank you in advance,
> Giovanni Petris
>
> --
>
> Giovanni Petris <[EMAIL PROTECTED]>
> Associate Professor
> Department of Mathematical Sciences
> University of Arkansas - Fayetteville, AR 72701
> Ph: (479) 575-6324, 575-8630 (fax)
> http://definetti.uark.edu/~gpetris/
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] GLMM and data manipulation (2nd try)
On Sat, May 03, 2008 at 06:43:22AM +1000, Andrew Robinson wrote:
> On Fri, May 02, 2008 at 03:06:31PM -0500, Giovanni Petris wrote:
> >
> > Hello,
> >
> > I posted a question yesterday but I got no replies, so I'll try to
> > reformulate it in a more concise way.
> >
> > I have the following data, summarizing approval ratings on two
> > different surveys for a random sample of 1600 individuals:
> >
> > > ## Example: Ratings of prime minister (Agresti, Table 12.1, p.494)
> > > rating <- matrix(c(794, 86, 150, 570), 2, 2)
> > > dimnames(rating) <- list(First = c("approve", "disapprove"),
> > + Second = c("approve", "disapprove"))
> > > rating
> > Second
> > Firstapprove disapprove
> > approve794150
> > disapprove 86570
> >
> > I would like to fit a logit model with approve/disapprove as response,
> > survey (first/second) as a fixed effect, and subject as a random
> > effect.
> >
> > 1) Is it possible to fit such a model directly using "lmer"?
> >
> > or
> >
> > 2) Should I unroll the table above into a dataframe containing also
> >fictitious subject id's? If this is the case, what is a clean way
> >of doing it?
>
>
> Unroll it.
>
> Asking for a "clean" way to do something is a disincentive because it
> implies that you know how to do it but not cleanly. In the future I
> would suggest that you do one of two things
>
> a) post your on dirty version and ask for a cleaner one, or
>
> b) just ask for something that works.
>
> Something like ...
Or (possibly better) make time a factor.
ratings <-
data.frame(
response = c(
rep(c(1,1), 794),
rep(c(1,0), 150),
rep(c(0,1), 86),
rep(c(0,0), 570)),
time = factor(rep(c(1,2), 1600)),
subject=rep(1:1600, each=2))
test.lmer <- lmer(response ~ time + (1|subject), data=ratings,
family=binomial)
Andrew
> but I don't know if you think that's clean or not.
>
> Andrew
>
> > Thank you in advance,
> > Giovanni Petris
> >
> > --
> >
> > Giovanni Petris <[EMAIL PROTECTED]>
> > Associate Professor
> > Department of Mathematical Sciences
> > University of Arkansas - Fayetteville, AR 72701
> > Ph: (479) 575-6324, 575-8630 (fax)
> > http://definetti.uark.edu/~gpetris/
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> --
> Andrew Robinson
> Department of Mathematics and StatisticsTel: +61-3-8344-6410
> University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
> http://www.ms.unimelb.edu.au/~andrewpr
> http://blogs.mbs.edu/fishing-in-the-bay/
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] improvement of Ancova analysis
Tobi, I think that it would be easier to provide advice if you were more explicit on what the model will be used for, and what is the structure of the data. Is there only one measurement for each marsupial? Is the goal to a) produce a model to predict marsupial weight given other variables, and if so, why, b) produce a model to estimate the effect of introduction on weight, with the other variables being nuisance variables, and if so, why, c) something else (and why) because all these factors affect the position that you could adopt about the questions that you have. Andrew On Sun, May 04, 2008 at 04:00:16AM +0200, Tobias Erik Reiners wrote: > Dear Helpers, > > I just started working with R and I'm a bit overloaded with information. > > My data is from marsupials reindroduced in a area. I have weight(wt), > hind foot > lenghts(pes) as continues variables and origin and gender as categorial. > condition is just the residuals i took from the model. > > >names(dat1) > [1] "wt" "pes" "origin" "gender" "condition" > > my model after model simplification so far: > model1<-lm(log(wt)~log(pes)+origin+gender+gender:log(pes)) > -->six intercepts and two slopes > > the problem is i have some things I can't include in my analysis: > 1.Very different sample sizes for each of the treatments > >tapply(log(wt),origin,length) > captivesitewild > 119 149 19 > 2.Substantial differences in the range of values taken by the > covariate (leg length) between treatments > >tapply(pes,origin,var) > captive site wild > 82.43601 71.2 60.42544 > >tapply(pes,origin,mean) > captive site wild > 147.3261 144.8698 148.2895 > > 4.Outliers > 5.Poorly behaved residuals > > thanks for the answer I am open minded to any different kind of analysis. > > Tobi > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] lme() with two random effects
Hi Mike,
I strongly suggest that you study Pinheiro and Bates (2000) to help
you make good decisions on the model specification and subsequent
steps.
In the meantime, you might find that exploring
lme(
x ~ gender*con*int*tone*cue
, random = ~ age | sub
, data = a
)
is helpful.
Cheers
Andrew
On Fri, May 09, 2008 at 04:09:07PM -0300, Mike Lawrence wrote:
> Hi all,
>
> I have collected response time data from 178 participants ('sub') for
> each combination of 4 within-Ss factors ('con','int','tone','cue').
> Additionally, I have recorded the gender of each participant, so this
> forms a between-Ss factor ('gender'). Normally this would be analyzed
> using aov:
>
> a=read.table('http://tinyurl.com/4pv5mh')
> a$sub = factor(a$sub)
> aov( #this may take a while!
> x~gender*con*int*tone*cue+Error(sub/(con*int*tone*cue))
> ,data=a
> )
>
> However, I'd like to also look at any effects of age, recorded in
> months ('age'). Since subjects were randomly sampled, their age should
> be a random effect and it is furthermore an unbalanced factor (N
> varies across levels of age). So I'm seeking the proper formulas to
> use in lme(). The following attempt yields the following error:
>
> library('nlme')
> lme(
> x ~ gender*con*int*tone*cue
> , random = ~ 1 | sub*age
> , data = a
> )
>
> Error in getGroups.data.frame(dataMix, groups) :
> Invalid formula for groups
>
>
> I would be very grateful for any suggestions.
>
> Cheers,
>
> Mike
>
> --
> Mike Lawrence
> Graduate Student, Department of Psychology, Dalhousie University
>
> Website: http://memetic.ca
>
> Public calendar subscribe link for iCal users:
> webcal://icalx.com/public/informavore/Public.ics
>
> "The road to wisdom? Well, it's plain and simple to express:
> Err and err and err again, but less and less and less."
> - Piet Hein
>
> ______
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lme nesting/interaction advice
On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote: > > The main point of my question is, having a 3 way anova (or ancova, if > you prefer), with *no* nesting, 2 fixed effects and 1 random effect, > why is it so boneheaded difficult to specify a bog standard fully > crossed model? I'm not talking about some rarified esoteric model > here, we're talking about stuff tought in a first year Biology Stats > course here[1]. That may be so, but I've never needed to use one. If it's bog-standard and yet boneheaded difficult, then presumably someone else would have had this problem before you. Perhaps a search of the archives will help? If you try, you will find many qualifiers to the effect that "lme isn't very well set up for crossed random effects". > Now, to avoid any chances of being misunderstood in my use of the > words 'fully crossed model', what I mean is a simple > > y ~ effect1 * effect2 * effect3 > > with effect3 being random (all all the jazz that comes from this > fact). I fully apprecciate that the only reasonable F-tests would be > for effect1, effect2 and effect1:effect2, but there is no way I can > use lme to specify such simple thing without getting the *wrong* > denDF. I need light on this topic and I'd say it's a general enough > question not to need much more handholding than this. Perhaps there are some circumstances unique to your situation. > I fully apprecciate that R is developed for love, not money, ... as is the R-help community ... > and if I > knew how to write an user friendly frontend for nlme and lme4 (and I > knew how to actually get the model I want) I'd be pretty happy to do > so and submit it as a library. In any case, I feel my complaint is > pefectly valid, because specifying such basic model should ideally > not such a chore, and I think the powers that be might actually find > some use from user feedback. This is not feedback. It is a compliant. But, the complaint boils down to the fact that you don't know what you're doing, and you show no evidence of having searched the R-help archives. How is that helpful? > Once I have sorted how to specify such trivial model I'll face the > horror of the nesting, in any case I attach a toy dataset I created > especially to test how to specify the correct model (silly me). Well, these data seem to differ. Is replica block? If not, then how can we reproduce your results? And, if I assume that it is, then the output df differ from what you sent in your original mail. So, I find this confusing. Then, from your original mail, > The easiest model ignores the nested random effects and uses just > selection, males and replica and the relative interactions. The > model > >lme(y ~ selection * males, random = ~1|replica/selection/males, mydata) forgive me, but I seem to see nesting in the random statement. That is what happens when we separate factors with a '/'; they are nested. We would expect that statement to not provide the correct df for the bog-standard fully crossed design. Perhaps if you were to comply with the request at the bottom of each R-help email, and provide commented, minimal, self-contained, reproducible code, that actually ran, ideally with fewer value judgements, you might get more attention from the people who are smarter than you and me, but have less time than either of us. Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] lme nesting/interaction advice
On Mon, May 12, 2008 at 10:34:40AM +1200, Rolf Turner wrote: > > On 12/05/2008, at 9:45 AM, Andrew Robinson wrote: > > >On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote: > >> > >>The main point of my question is, having a 3 way anova (or ancova, if > >>you prefer), with *no* nesting, 2 fixed effects and 1 random effect, > >>why is it so boneheaded difficult to specify a bog standard fully > >>crossed model? I'm not talking about some rarified esoteric model > >>here, we're talking about stuff tought in a first year Biology Stats > >>course here[1]. > > > >That may be so, but I've never needed to use one. > > So what? This is still a standard, common, garden-variety > model that you will encounter in exercises in many (if not > all!) textbooks on experimental design and anova. To reply in similar vein, so what? Why should R-core or the R community feel it necessary to reproduce every textbook example? How many times have *you* used such a model in real statistical work, Rolf? > >If it's bog-standard and yet boneheaded difficult, then presumably > >someone else would have had this problem before you. Perhaps a search > >of the archives will help? If you try, you will find many qualifiers > >to the effect that "lme isn't very well set up for crossed random > >effects". > > But that avoids the question as to *why* it isn't very well > set up for crossed random effects? What's the problem? > What are the issues? The model is indeed bog-standard. > It would seem not unreasonable to expect that it could be > fitted in a straightforward manner, and it is irritating to > find that it cannot be. If SAS and Minitab can do it at > the touch of a button, why can't R do it? Bates has made no secret of the fact that lme was intended first and foremost for nested designs, and that support for crossed designs is not promised. He has said so on many occasions, as a search would find. He is now working on lme4, which will support crossed designs. It's not done yet. > >and you show > >no evidence of having searched the R-help archives. How is that > >helpful? > > It doesn't seem to me to be a complaint as such. It is a > request for insight. I too would like some insight as to > what on earth is going on. And why do you say Federico > shows no evidence of having searched the archives? One can > search till one is blue in the face and come away no wiser > on this issue. At least one can know that there is an issue, which apparently Federico previously did not. Warm wishes Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] lme nesting/interaction advice
On Mon, May 12, 2008 at 10:50:03AM +0100, Federico Calboli wrote: > > On 12 May 2008, at 01:05, Andrew Robinson wrote: > > >On Mon, May 12, 2008 at 10:34:40AM +1200, Rolf Turner wrote: > >> > >>On 12/05/2008, at 9:45 AM, Andrew Robinson wrote: > >> > >>>On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote: > >>>> > >>>>The main point of my question is, having a 3 way anova (or > >>>>ancova, if > >>>>you prefer), with *no* nesting, 2 fixed effects and 1 random > >>>>effect, > >>>>why is it so boneheaded difficult to specify a bog standard fully > >>>>crossed model? I'm not talking about some rarified esoteric model > >>>>here, we're talking about stuff tought in a first year Biology > >>>>Stats > >>>>course here[1]. > >>> > >>>That may be so, but I've never needed to use one. > >> > >>So what? This is still a standard, common, garden-variety > >>model that you will encounter in exercises in many (if not > >>all!) textbooks on experimental design and anova. > > > >To reply in similar vein, so what? Why should R-core or the R > >community feel it necessary to reproduce every textbook example? How > >many times have *you* used such a model in real statistical work, > >Rolf? > > There is a very important reason why R (or any other stats package) > should *easily* face the challenge of bog standard models: because it > is a *tool* for an end (i.e. the analysis of data to figure out what > the heck they tell us) rather than a end in itself. But a tool that mostly (entirely?) appears in textbooks. > Bog standard models are *likely* to be used over and over again > because they are *bog standard*, and they became such by being used > *lots*. Well. I have documentation relevant to nlme that goes back about 10 years. I don't know when it was first added to S-plus, but I assume that it was about then. Now, do you think that if the thing that you want to do was really bog standard, that noone would have raised a fuss or solved it within 10 years? > If someone with a relatively easy model cannot use R for his job s/he > will use something else, and the R community will *not* increase in > numbers. Since R is a *community driven project*, you do the math on > what that would mean in the long run. Fewer pestering questions? ;) Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Cumulative average
Jacques, you should be able to construct a solution from cumsum(). Cheers Andrew On Wed, May 21, 2008 at 08:48:29PM -0500, Jacques Wagnor wrote: > Dear List, > > Does there exist a function that calculates a cumulative average? > Neither running() from library(gregmisc) nor running.mean() from > library(igraph) seems to be able to give a cumulative average. > > Any help or pointers would be greatly appreciated. > > Jacques > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Some problems with Sweave
Martin, try omitting the results=tex argument. Andrew On Fri, May 23, 2008 at 10:16:33AM +0200, [EMAIL PROTECTED] wrote: > Dear R users, > I'm working in a brief R-tutorial to a group of students. To make that I'm > using Sweave but I've got two problems: > > First, I want show how R operates with the matrix type but, I write in the > .rnw document the code > > <>= > matriz <- matrix(vector,nrow=3,ncol=6) > matriz > @ > > and after compilating the LaTex document I obtain in the pdf the next text > > > matriz <- matrix(vector, nrow = 3, ncol = 6) > > matriz > [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 4 1 4 1 4 [2,] 2 5 2 5 2 5 [3,] 3 6 3 > 6 3 6 > > My question is, How must I do To obtain in the pdf somethin near to > > > matriz <- matrix(vector, nrow = 3, ncol = 6) > > matriz > [,1] [,2] [,3] [,4] [,5] [,6] > [1,]141414 > [2,]252525 > [3,]363636 > > I`ve tought in xtable but the aspect is not as the R console. > > On the other hand I want show the list type R-treatment, the problem here > is in the LaTex compilation > I write in the .rnwd document > > <>= > lista <- list(cadena='String',vector=c(1,1,1),logica=TRUE) > lista$cadena > @ > > the Sweave call is ok, but when I compile the .tex document, It produces > errors caused by the $ simbols. Anybody knows how save this problem? > > Thanks in advance, > Mart?n Gast?n > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Need help for nlme
My advice is to try to simplify this as much as possible. When it is
as simple as possible, it will either work or not work. If it works,
then you have learned something useful. If it does not work, then
send your question again. Right now there is a great deal of detail
that may or may not be extraneous.
Andrew.
On Sun, May 25, 2008 at 01:36:49PM -0500, [EMAIL PROTECTED] wrote:
> Hi everyone,
> I try to write a module based on nlme however R always shows me the
> error message
> Error in eval(expr, envir, enclos) : object "y" not found. Does anyone
> know how to solve this? There is no problem in nls at all.
>
> require(nlme)
> AMPmixed<-function(y, x, S1=c("asymptotic","logistic"),
> S2=c("asymptotic","logistic"), data, start,random)
> {
> logist.logist<-function(x,alpha,delta,psi.l,tau.l,gamma,h){
>
> delta+(alpha-delta+gamma*(x-(h-1))/exp(x))/(1+exp(-(x-tau.l)/psi.l))}
> logist.asymp<-function(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h){
>
> delta+(alpha-delta)/(1+exp(-(x-tau.l)/psi.l))+(gamma*(x-(h-1))/exp(x))*exp(-exp(1/lpsi.a)*x)}
> asymp.asymp<-function(x,alpha,delta,lpsi.a,gamma,h){
>
> delta+(alpha-delta)*exp(-exp(1/lpsi.a)*x)+(gamma*(x-(h-1))/exp(x))*exp(-exp(1/lpsi.a)*x)}
> asymp.logist<-function(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h){
>
> delta+(alpha-delta)*exp(-exp(1/lpsi.a)*x)+(gamma*(x-(h-1))/exp(x))/(1+exp(-(x-tau.l)/psi.l))}
>
>(logistic.logistic<-function(y, x, data, start, random){
>
> nlme.out<-nlme(y~logist.logist(x,alpha,delta,psi.l,tau.l,gamma,h),
> data=data, start=start,
> fixed=alpha+delta+psi.l+tau.l+gamma+h~1,
> random=random)
>list(nlme.out=summary(nlme.out))
>})
>(logistic.asymptotic<-function(y, x, data, start, random){
>
> nlme.out<-nlme(y~logist.asymp(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h),
> data=data, start=start,
>
> fixed=alpha+delta+psi.l+tau.l+lpsi.a+gamma+h~1, random=random)
> list(nlme.out=summary(nlme.out))
>})
>(asymptotic.logistic<-function(y, x, data, start,random){
>
> nlme.out<-nlme(y~asymp.logist(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h),
> data=data, start=start,
>
> fixed=alpha+delta+psi.l+tau.l+lpsi.a+gamma+h~1, random=random)
> list(nlme.out=summary(nlme.out))
>})
>(asymptotic.asymptotic<-function(y, x, data, start, random){
>
> nlme.out<-nlme(y~asymp.asymp(x,alpha,delta,lpsi.a,gamma,h), data=data,
> start=start,
> fixed=alpha+delta+lpsi.a+gamma+h~1,random=random)
>list(nlme.out=summary(nlme.out))
>})
>
>if(S1=="logistic" && S2=="logistic")
> {(AMPmixed=logistic.logistic(y, x, data, start, random))}
>else if(S1=="logistic" &&
> S2=="asymptotic"){(AMPmixed=logistic.asymptotic(y, x, data, start,
> random))}
>else if(S1=="asymptotic" &&
> S2=="logistic"){(AMPmixed=asymptotic.logistic(y, x, data, start,
> random))}
>else if(S1=="asymptotic" &&
> S2=="asymptotic"){(AMPmixed=asymptotic.asymptotic(y, x, data, start,
> random))}
>}
>
> #
> con rep biomass
> 10.00 1 1.126
> 20.32 1 1.096
> 31.00 1 1.163
> 43.20 1 0.985
> 5 10.00 1 0.716
> 6 32.00 1 0.560
> 7 100.00 1 0.375
> 80.00 2 0.833
> 90.32 2 1.106
> 10 1.00 2 1.336
> 11 3.20 2 0.754
> 12 10.00 2 0.683
> 13 32.00 2 0.488
> 14 100.00 2 0.344
>
> iso<-read.table(file="E:\\Hormesis\\data\\isobutanol.txt", header=T)
> aa<-groupedData(biomass~con|rep, data=iso)
> van2<-AMPmixed(y=biomass, x=con, S1="asymptotic", S2="asymptotic", data=aa,
> random=pdDiag(alpha+delta+lpsi.a+gamma+h~1),
>start=c(alpha= 0.7954, delta= 0.3231, lpsi.a=-0.2738,
> gamma= 1.0366, h=0.8429))
> van2
>
> Thank you very much in advance.
> Chunhao
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented,
Re: [R] How to remove NAs and lme function
Jen, try na.action = na.exclude Andrew On Wed, May 28, 2008 9:26 pm, Jen_mp3 wrote: > > I am working on a project to find a model for the concentration of > dissolved > oxygen in the river clyde. Ive fitted a linear mixed model as > lme(DOW~Temperature+Salinity+Year+factor(Station)*factor(Depth), > random~1|id), where id is an identifier of the day over 20 years defined > as > Day*1 + Month*100 + (1900 - Year). > Anyway, there are some NAs for the concentration of dissolved oxygen in > the > water so I know you add in na.action = na.omit and that omits the NAs so > there are 9008 observations in the model, but it doesnt do it for the > whole > data set where there are 10965 including observations with NAs. I would > like > to plot the residuals from the model against the Salinity, Temperature and > Year, but when I try, it seems to want to take the observations of these > variables from the full data set and the residuals from the model which of > course doesnt work. I have tried using > data1 <- data[data$DOW != "NA",] on the whole data set but it doesnt work. > How can I remove the NAs from a data set? > > -- > View this message in context: > http://www.nabble.com/How-to-remove-NAs-and-lme-function-tp17510564p17510564.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. > Andrew Robinson Senior Lecturer in Statistics Tel: +61-3-8344-6410 Department of Mathematics and StatisticsFax: +61-3-8344 4599 University of Melbourne, VIC 3010 Australia Email: [EMAIL PROTECTED]Website: http://www.ms.unimelb.edu.au __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] can I do this with R?
On Wed, May 28, 2008 at 03:47:49PM -0700, Xiaohui Chen wrote: > Frank E Harrell Jr ??: > >Xiaohui Chen wrote: > >>step or stepAIC functions do the job. You can opt to use BIC by > >>changing the mulplication of penalty. > >> > >>I think AIC and BIC are not only limited to compare two pre-defined > >>models, they can be used as model search criteria. You could > >>enumerate the information criteria for all possible models if the > >>size of full model is relatively small. But this is not generally > >>scaled to practical high-dimensional applications. Hence, it is often > >>only possible to find a 'best' model of a local optimum, e.g. > >>measured by AIC/BIC. > > > >Sure you can use them that way, and they may perform better than other > >measures, but the resulting model will be highly biased (regression > >coefficients biased away from zero). AIC and BIC were not designed to > >be used in this fashion originally. Optimizing AIC or BIC will not > >produce well-calibrated models as does penalizing a large model. > > > Sure, I agree with this point. AIC is used to correct the bias from the > estimations which minimize the KL distance of true model, provided the > assumed model family contains the true model. BIC is designed for > approximating the model marginal likelihood. Those are all > post-selection estimating methods. For simutaneous variable selection > and estimation, there are better penalizations like L1 penalty, which is > much better than AIC/BIC in terms of consistency. Xiaohui, Tibshirani (1996) suggests that the quality of the L1 penalty depends on the structure of the dataset. As I recall, subset selection was preferred for finding a small number of large effects, lasso (L1) for finding a small to moderate number of moderate-sized effects, and ridge (L2) for many small effects. Can you provide any references to more up-to-date simulations that you would recommend? Cheers, Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Bimodal Distribution
Hi Mike, if you can decompose the bimodal distribution into (eg two) known forms, then you could try a stepwise approach, eg: If uniform < 0.5 then double it and use it to draw from the inverse cdf of A, else, double (uniform - 0.5) and use it to draw from the inverse cdf of B. You can change the cutoff and the weights to suit your need. It's best to double-check by plotting an empirical density of the numbers generated. I hope that this helps, Andrew On Thu, May 29, 2008 at 04:05:29PM -0600, Mike Williams wrote: > Hello R Users, > > I am doing a Latin Hypercube type simulation. I have found the > improvedLHS function and have used it to generate a bunch of properly > distributed uniform probabilities. Now I am using functions like qlnorm > to transform that into the appropriately lognormal or triangularly > distributed parameters for my modes. However I have a parameter which I > believe is bimodally distributed, could someone please point me at an > appropriate function equivalent to qlnorm which I can use, because for > some reason I have been unable to find one. It occurs to me that maybe > one doesn't exist, in which case could someone give me some other idea > of how to accomplish this goal? > > Thanks, > > Mike > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] "False convergence" in LME
These are great tips from Spencer.
The other thing that I have found useful is to use a different
optimizing algorithm. You can do this by:
control=lmeControl(opt = "optim")
Good luck!
Andrew
On Sat, Jun 14, 2008 at 09:45:22AM -0700, Spencer Graves wrote:
> This is a common problem, for which solutions are poorly documented.
>
> 1. Have you tried fitting simpler models, in the hopes of
> getting something that converges without complaint, then use 'update' to
> try more complicated models? It sometimes works, though often not.
>
> 2. Also, have you tried something like 'lme(...,
> control=lmeControl(returnObject=TRUE))'? If no, I suggest you try this
> first. With the error message you report, I would expect this to work,
> though it may not.
>
> 3. Finally, have you tried something like 'lme(...,
> control=lmeControl(singular.ok=TRUE))' OR 'lme(...,
> control=lmeControl(singular.ok=TRUE, returnObject=TRUE))'? I'm not
> sure, but I believe this may work when "returnObject=TRUE" does not.
>
> 4. Have you tried the following:
>
>library(FinTS)
>package.dir('nlme')
>
>I tried this just now and learned that the 'nlme'
> packages contains two non-standard subdirectories named "mlbook" and
> "scripts". The second contains files names 'ch01.R', 'ch02.R', etc.,
> which contain the R commands required to reproduce virtually all the
> figures, tables and examples in Pinheiro and Bates (2000) Mixed-Effects
> Models in S and S-PLUS (Springer). If you have not already done so, I
> recommend you get this book and use these script files to facilitate
> your study of it. Doug Bates is one of the world's leading experts in
> mixed-effects modeling, and I have learned a lot from my study of it.
>
> Hope this helps.
> Spencer Graves
>
> Rebecca Sela wrote:
> >I tried to use LME (on a fairly large dataset, so I am not including it),
> >and I got this error message:
> >
> >Error in lme.formula(formula(paste(c(toString(TargetName),
> >"as.factor(nodeInd)"), : nlminb problem, convergence error code = 1
> > message = false convergence (8)
> >
> >Is there any way to get more information or to get the potentially wrong
> >estimates from LME?
> >
> >(Also, the page in the NLMINB documentation,
> >http://netlib.bell-labs.com/cm/cs/cstr/153.pdf, has errors in it, which
> >makes it harder to check on what is happening.)
> >
> >Thank you in advance!
> >
> >Rebecca
> >
> >__
> >R-help@r-project.org mailing list
> >https://stat.ethz.ch/mailman/listinfo/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.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 a data.frame
Hi Joe,
cbind coerces the data to be the same type. Much nicer is:
d <- data.frame(x=x, y=y, z=z)
Cheers
Andrew
On Wed, Feb 13, 2008 at 05:17:32PM -0500, Joe Trubisz wrote:
> OK...newbie question here.
> Either I'm reading the docs wrong, or I'm totally confused.
>
> Given the following:
>
> x<-c("aaa","bbb","ccc")
> y<-rep(0,3)
> z<-rep(0,3)
>
> is.character(x)
> [1] TRUE
>
> is.numeric(y)
> [1] TRUE
>
> Now...I want to create a data frame, but keep the data types.
> In reading the docs, I assume you do it this way:
>
> d<-data.frame(cbind(x=I(x),y=y,z=z)
>
> But, when I do str(d), I get the following:
>
> 'data.frame': 3 obs. of 3 variables:
> $ x: Factor w/ 3 levels "aaa","bbb","ccc": 1 2 3
> $ y: Factor w/ 1 level "0": 1 1 1
> $ z: Factor w/ 1 level "0": 1 1 1
>
> I thought the I() prevents character from becoming factors, right?
> Secondly, how do I force y and z in the data frame to become numeric?
>
> Thanks in advance
> Joe
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 message has been scanned for viruses and
> dangerous content by MailScanner, and is
> believed to be clean.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] skip non-converging nls() in a list
Hi Marc,
use try()
Cheers
Andrew
On Mon, Feb 18, 2008 at 01:45:35PM -0500, Marc Belisle wrote:
> Howdee,
>
> My question appears at #6 below:
>
> 1. I want to model the growth of each of a large number of individuals using
> a 4-parameter logistic growth curve.
>
> 2. nlme does not converge with the random structure that I want to use.
>
> 3. nlsList does not converge for some individuals.
>
> 4. I decided to go around nlsList using:
>
> t(sapply(split(data, list(data$id)),
> function(subd){coef(nls(mass ~ SSfpl(age, A, B, xmid, scal), data =
> subd))}))
>
> 5. This does not converge either:
>
> 'Error in nls(y ~ cbind(1, 1/(1 + exp((xmid - x)/exp(lscal, data = xy,
> :
> singular gradient'
>
> 6. Would anyone know if I can get R to skip non-converging use of nls() so
> that I can at least obtain the parameters for the curves which R can fit?
>
> Thanks for your time,
>
> Marc
>
> ===
> Marc B?lisle
> Professeur adjoint
> Chaire de recherche du Canada en ?cologie spatiale et en ?cologie du paysage
> D?partement de biologie
> Universit? de Sherbrooke
> 2500 Boul. de l'Universit?
> Sherbrooke, Qu?bec
> J1K 2R1 Canada
>
> T?l: +1-819-821-8000 poste 61313
> Fax: +1-819-821-8049
> Courri?l: [EMAIL PROTECTED]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 message has been scanned for viruses and
> dangerous content by MailScanner, and is
> believed to be clean.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 more idiomatic way to write this
How about
x <- x / rep(divs, rep(1000, 1000))
?
Cheers,
Andrew
On Sun, Feb 24, 2008 at 10:36:23PM -0300, Andre Nathan wrote:
> Hello,
>
> I have a vector of 1,000,000 numbers and another vector of 1,000
> divisors. What I'd like to do is to divide the first 1,000 numbers of
> the first vector by the first divisor, then the next 1,000 by the second
> divisor and so on. I came up with this, but I was wondering if there is
> a more idiomatic, R-like way to write it:
>
> x <- ...
> divs <- ...
>
> for (i in seq(from = 1, to = 100, by = 1000)) {
> x[i:(i - 1 + 1000)] <- x[i:(i - 1 + 1000)] / divs[i %/% 1000 + 1]
> }
>
> Any suggestions are welcome.
>
> Thanks in advance,
> Andre
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 message has been scanned for viruses and
> dangerous content by MailScanner, and is
> believed to be clean.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] discrete variable
Pete,
try
x <- c(0, 3, 4, 3, 15, 5, 6, 5)
table(x)
or
length(table(x))
Cheers
Andrew
On Sun, Mar 02, 2008 at 09:27:20PM +0100, Pete Dorothy wrote:
> Hello,
>
> I am sorry for asking such a basic question. I could not find an answer to
> it using google.
>
> I have a discrete variable (a vector x) taking for example the following
> values : 0, 3, 4, 3, 15, 5, 6, 5
>
> Is it possible to know how many different values (modalities) it takes ?
> Here it takes 6 different values but the length of the vector is 8.
>
> I would like to know if there is a way to get the set of the
> modalities {0,3,4,15,5,6} with the number of times each one is taken
> {1,2,1,1,2,1}
>
> Thank you very much
>
> P.S. : is there some useful functions for discrete variables ?
>
> [[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.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] "normal" (as in "Guassian")
On Mon, Mar 03, 2008 at 10:22:41PM +0100, Peter Dalgaard wrote: > Patrick Burns wrote: > > Douglas Bates wrote: > > > > > >> On Mon, Mar 3, 2008 at 8:25 AM, Duncan Murdoch <[EMAIL PROTECTED]> wrote: > >> > >> > >> > >>> On 3/3/2008 9:10 AM, Rogers, James A [PGRD Groton] wrote: > >>> > >>>> As someone of partly French heritage, I would also ask how this > >>>> distribution came to be called "Gaussian". It seems very unfair to de > >>>> Moivre, who discovered the distribution at least half a century earlier. > >>>> :-) > >>>> > >>> Just an example of Stigler's Law. > >>> > >>> > >>> > >> Taking this to a whole new level of "off topic", I wonder if Stigler's > >> Law is self-referential? That is, should Stigler's Law more correctly > >> be attributed to someone else? > >> > > > > No. If Stigler's Law were named after some prior person, > > then it wouldn't be an example of itself. > > > Only if said person actually was first to discover it, surely. I believe that Stigler believes that he was not the first to discover Stigler's Law. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Great difference for piecewise linear function between R and SAS
On Mon, Mar 24, 2008 at 01:09:35PM -0500, Paul Johnson wrote: > On Mon, Mar 24, 2008 at 8:11 AM, Duncan Murdoch <[EMAIL PROTECTED]> wrote: > > On 3/24/2008 9:03 AM, zhijie zhang wrote: > > > Dear Murdoch, > > > > > Thanks very > > > > I would use predict() instead. What you have there doesn't look as > > though it uses the B-spline basis. > > > > The reference given in the ?bs help page is a reasonable starting point, > > but just about any book that covers splines should handle the B-spline > > basis and the linear case. > > > > Duncan Murdoch > > > > Dear Duncan and others: > > Can you please refer us to an understandable book that explains about > b-splines? I've tried a few and the math pretty quickly becomes > unintelligible to a non-math major. I think you should try Simon Wood's "Generalized Additive Models: An introduction in R." 2006 Chapman and Hall. I reviewed it recently and I think it's really very good. Cheers Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] peak finding
Another approach that involves more infrastructure is to fit a smooth line to your points, compute the first derivative, and look for change in sign in the first derivative. eg x <- c(14,15,12,11,12,13,14,15,16,15,14,13,12,11,14,12) smoothed.dx <- predict(smooth.spline(x), deriv=1)$y which(c(smoothed.dx,NA) > 0 & c(NA, smoothed.dx) < 0) My experience is that this approach sometimes requires some fine-tuning, eg in fitting the smooth. I hope that this helps Andrew On Mon, Mar 24, 2008 at 11:23:57PM -0400, Research Scholar wrote: > Hi > Thanks for replying. I meant x[4] is the start of a peak shape and x[14] is > the end of that peak and x[9] is the maxima of the peak. > Thanks, > John > > > > > On Mon, Mar 24, 2008 at 11:09 PM, <[EMAIL PROTECTED]> wrote: > > > It's hard to see how positions 4 and 14 correspond to 'peaks', they look > > like troughs to me. So perhaps this is what you mean: > > > > > x <- c(14,15,12,11,12,13,14,15,16,15,14,13,12,11,14,12) > > > > > y <- which(x == min(x)) > > > y > > [1] 4 14 > > > > as a function: > > > > somefunction <- function(x) which(x == min(x)) > > > > > > Bill Venables > > CSIRO Laboratories > > PO Box 120, Cleveland, 4163 > > AUSTRALIA > > Office Phone (email preferred): +61 7 3826 7251 > > Fax (if absolutely necessary): +61 7 3826 7304 > > Mobile: +61 4 8819 4402 > > Home Phone: +61 7 3286 7700 > > mailto:[EMAIL PROTECTED] > > http://www.cmis.csiro.au/bill.venables/ > > > > -Original Message- > > From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] > > On Behalf Of Research Scholar > > Sent: Tuesday, 25 March 2008 12:54 PM > > To: r-help@r-project.org > > Subject: [R] peak finding > > > > Hi all > > Is there a function that can find the start and end position of peaks > > in a > > set of numbers. > > eg. > > x <- c(14,15,12,11,12,13,14,15,16,15,14,13,12,11,14,12) > > y <- somefunction(x) > > > > y > > 4 14 > > > > > > Thanks > > John > > > >[[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<http://www.r-project.org/posting-guide.html> > > and provide commented, minimal, self-contained, reproducible code. > > > > > > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] [Stat related] Understanding Portmanteau test
On Sat, Nov 08, 2008 at 09:58:14PM -0800, RON70 wrote: > > Still waiting for some input. Did my question void forum rule in any manner? No, no specific rule. It's just not a particularly easy question to answer, because it's not clear in what context you have seen the phrase "portmanteau test". Any given answer might be completely irrelevant and useless to you. > > Sorry to be off-topic. Can somebody please explain me what is Portmanteau > > test? Why it's name is like that? When I would say, a particular test is > > portmanteau test? I did some googling but got no satisfactory answer at > > all. Please anybody help for understanding that? For what it's worth, the Cambridge Dictionary of Statistics defines a "portmanteau test" as a test for assessing the fit of models for time series in the presence of outliers. It provides a citation: Computational Statistics, 1994, 9, 301-310. I hope that helps. Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Spatial ANCOVA in R
Hi Camilo, try gls() in the nlme package. Andrew On Sat, Nov 15, 2008 at 11:34:57PM -0400, Camilo Mora wrote: > Hi: > > Does anyone know if it is possible to run an ANCOVA in R while accounting or > controlling for spatial autocorrelation? I have found usefull information into > how to account for spatial autocorrelaion in regression models but not much > into how to deal with the problem in an ANCOVA. > > Thanks, > > Camilo > > Camilo Mora, Ph.D. > SCRIPPS Institute of Oceanography > University of California San Diego > San Diego, USA > Phone: (858) 822 1642 > http://cmbc.ucsd.edu/People/Faculty_and_Researchers/mora/ > And > Department of Biology > Dalhouisie University > Halifax, Canada > Phone: (902) 494 3910 > http://as01.ucis.dal.ca/fmap/people.php?pid=53 > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Growth rate determination using ANCOVA
Hi David, On Fri, Nov 21, 2008 at 12:01:52PM -0800, dschruth wrote: > I'm a programmer in a biology lab who is starting to use R to automate > some of our statistical analysis of growth rate determination. But I'm > running into some problems as I re-code. > > 1) Hypotheses concerning Slope similarity/difference: >I'm using R's anova(lm()) methods to analyse a model which looks > like this: > growth.metric ~ time * test.tube >I understand that testing the the interaction between time and tube > (time:test.tube) will tell us if the growth rates (for the last three > test tubes) are significantly different from one another (Ho=slopes > are the same). The purpose of doing this test is so that we can be > certain our cultures have fully acclimated to the treatment and aren't > going to change much if we stop measuring. This is an important cost > saving practice too as measurements can go on for years. Yet I'm > worried that our null and alternative hypotheses should be swapped so > that our test is more conservative (Ho=slopes are different ... ie > still acclimating.) Good thinking. > Is there a way to specify my model that flips these hypotheses? > Should I be using a different method? Is this even appropriate? You could think about equivalence tests. See e.g. references in the equivalence package. > 2) Growth Rate is confounded with Variance of Growth Rate >I'm also worried about the fact that rates for cultures with faster > growth are calculated using fewer data points (assuming similar > sampling times between treatments) . The result is that growth ~ var > (growth). Not only does this put a wrinkle in my analysis between > treatments, but it also biases the growth acclimation determining > ANCOVA test above. Faster growing cultures will usually pass the "no > significant difference between slopes test" more easily because there > are fewer points from which to be certain about rejecting Ho. > > Is there a way to control for this? > Perhaps I could include the number of points in my model? Depending on the model that you apply, you might be able to explicitly model the variance to allow for this possibility. I would guess that it's not necessarily only the fewer data points contributing to the greater variation. Faster-growing cultures might also be inherently more variable. > 3) Statistical validity of using subsets of growth.metric measurements > within a test tube >There are some lab members who insist that we can throw out the > beginning and end of our log transformed growth.metric measurements > because they are outliers in determining maximum growth.I've > proposed looping through all possible combinations of 3 or more points > within the growth curve and using the highest or best fitting (best R- > squared) slope. But this idea has been rejected by our PI as not be a > valid thing to do. > > Ideas here? I'm feeling very cautious about giving advice on this question as I don't know enough about the area. Sorry. I hope that this helps, otherwise. Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] confusion with R syntax
Hi Mark,
what's happening here is that R is applying the one-dimensional
subscripting operations sequentially.
Try
x <- seq(1,10)
x[2:4][1]
Cheers
Andrew
On Thu, Oct 11, 2007 at 03:34:22PM -0400, Leeds, Mark (IED) wrote:
> I just noticed something by accident with R syntax that I'm sure is
> correct but I don't understand it. If I have
> a simple numeric vector x and I subscript it, it seems that I can then
> subscript a second time with TRUE
> or FALSE, sort of like a 2 dimensional array in C. Does someone know if
> this is documented somewhere
> Because it's neat but I never knew it existed. To me it seems like a 1
> dimensional vector should
> have only one dimensional indexing ?
>
> x <- seq(1,10)
> > x
> [1] 1 2 3 4 5 6 7 8 9 10
> > x[2:4][c(TRUE,FALSE,TRUE)]
> [1] 2 4
>
> But, it only works for TRUE or FALSE and not numbers so I think it's not
> really 2 dimensional indexing.
>
> x[1][2]
>
> [1] NA
>
> If someone could explain this mechanism or tell me what I should look
> for in the archives, it would
> be appreciated. Thanks.
>
>
> This is not an offer (or solicitation of an offer) to bu...{{dropped:22}}
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Repeated Measures/Linear Mixed Effects function
Hi Keith, it seems like a good starting position. I recommend that you spend some time studying Pinheiro and Bates's book to see where t ogo from here. Cheers Andrew On Mon, Oct 22, 2007 at 02:58:51PM -0800, Keith Cox wrote: > I have three columns of data, Xc, Trt and fish. This was a repeated > measures design with 6 measurements taken from each of 5 fish. Xc is the > actual measurement, Trt is the treatment, and fish is the fish number. Data > can be seen below (hopefully it is in the column format). I would like to > look for differences between treatments in a repeated measures format. I > used the following code > > > > library(nlme) > > > > summary(lme(Xc~trt,data=R.exp,random=~1|fish)) > > > > This seems to work and I would like to know if > > > > 1)this is the right function for my question, if so, then > > 2)in the summary, value is the first column, but what value is this. I > have searched in vain and cannot find the answer. Any help or links to > examples would be appreciated greatly. > > Fixed effects: Xc ~ trt > > Value Std.Error DF t-value p-value > > (Intercept) 103.62 2.619657 20 39.55480 0 > > trtB-33.28 1.991197 20 -16.71356 0 > > trtC-39.38 1.991197 20 -19.77705 0 > > trtD-32.60 1.991197 20 -16.37206 0 > > trtE-47.32 1.991197 20 -23.76460 0 > > trtF-39.58 1.991197 20 -19.87749 0 > > > > > Xc > > Trt > > fish > > > 109.1 > > a > > 1 > > > 73 > > b > > 1 > > > 68.4 > > c > > 1 > > > 74.8 > > d > > 1 > > > 60.3 > > e > > 1 > > > 57 > > f > > 1 > > > 106 > > a > > 2 > > > 72.3 > > b > > 2 > > > 67 > > c > > 2 > > > 70.6 > > d > > 2 > > > 58.2 > > e > > 2 > > > 66.2 > > f > > 2 > > > 102 > > a > > 3 > > > 67.1 > > b > > 3 > > > 61 > > c > > 3 > > > 68.4 > > d > > 3 > > > 50.2 > > e > > 3 > > > 64.7 > > f > > 3 > > > 105 > > a > > 4 > > > 76.6 > > b > > 4 > > > 68.8 > > c > > 4 > > > 77.7 > > d > > 4 > > > 61.8 > > e > > 4 > > > 75.7 > > f > > 4 > > > 96 > > a > > 5 > > > 62.7 > > b > > 5 > > > 56 > > c > > 5 > > > 63.6 > > d > > 5 > > > 51 > > e > > 5 > > > 56.6 > > f > > 5 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > Marlin Keith Cox Ph.D. > > At-Sea Processor Professorship of Fisheries Biology > Science Chair > > Sheldon Jackson College > > Sitka, Alaska 99835 > > 907.747.5296 > > http://www.sheldonjackson.edu > > > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Calculate aggregate differences
Hi Cristian, instead of aggregate, how about something like: n <- dim(my_agg)[1] my_agg$vol.diff <- my_agg$Vol - c(NA, my_agg$Vol[1:(n-1)] my_agg <- my.agg[my.agg$Age > min(my.agg$Age),] (assumes same minimum age for all treatments) (not checked) Cheers, Andrew On Thu, Nov 01, 2007 at 12:09:34PM -0400, [EMAIL PROTECTED] wrote: > Hi everyone > > I am trying to summarize a table with yield estimates of a forest > plantation. For that I have four blocks and four treatments measured over > a period of 10 years (every year). In each plot trees are measured > (diameters and heights). > > With aggregate function I can calculate the average diameter or the total > volume for each plot at any time with something like this: > > my_agg <- aggregate(list(Ht = HTO, Dbh = DBH), list(Age = AGE, Block = > BLOCK, Treat = TREAT), sum) > > my_agg <- aggregate(list(VOL = VOLUME), list(Age = AGE, Block = BLOCK, > Treat = TREAT), sum) > > where HTO is the height of the tree > Ht average height of trees in each plot at time t > DBH is the diameter at 1.3 meter height > Dbh is the average diameter of the plot at time t > VOL is the [i]th tree volume > Vol is the volume of the plot at time t > > To do actual growth analysis I need to calculate the difference in volume > between t and (t-1) for every block/treatment. > > Now the question > > How can I use aggregate to calculate the difference of each time series > and have an output like aggregate does? I know it is imposible to use > diff() within aggregate because it gives me an scalar. > In SAS proc means handles such a thing, and diff() in R is fine for a > single series, but in this case I want to calculate for all the > treatments/blocks at the same time. > > I can do the whole using by and some more less elegant procedures, but I > figure there should be a cleaner way as in PROC MEANS. > > Any suggestions? > > Cristian Montes > NC 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. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Calculate aggregate differences
Hi Cristian, yes, indeed. What I'm not sure about is how doing this by groups will give you a different result, if the object is already sorted according to the group structure. As you are working with differences, each group 'loses' an observation, which will be the observation that is calculated as the spurious difference between the first measurement of the group and the last measurement of the previous group. If that doesn't help, let me suggest that you construct a small worked example that will show us the kind of input and output that you are thinking about. Cheers, Andrew On Thu, Nov 01, 2007 at 12:55:40PM -0400, [EMAIL PROTECTED] wrote: > Thanks Andrew, but what you gave me can actually be done simpler like > > my_agg$vol.diff <- diff(my_agg$Vol) > > or > > my_agg$vol.diff <-c(NA, diff(my_agg$Vol)) # for a list with the same > #length as the aggregated mean list > > However, what I need is to be able to have this done by groups and then > have a table as an output. > > For example if I use by(): > > my_agg$vol <- with(my_agg, by(Vol, list(Block, Treatment), diff)) > > but then I am stucked with a by list that can't be coerced and the problem > becomes how to coerce the by list that has a vector for each group? > > Cristian Montes > NC State University > > > Hi Cristian, > > > > instead of aggregate, how about something like: > > > > > > n <- dim(my_agg)[1] > > my_agg$vol.diff <- my_agg$Vol - c(NA, my_agg$Vol[1:(n-1)] > > my_agg <- my.agg[my.agg$Age > min(my.agg$Age),] > > > > > > (assumes same minimum age for all treatments) > > > > (not checked) > > > > Cheers, > > > > Andrew > > > > On Thu, Nov 01, 2007 at 12:09:34PM -0400, [EMAIL PROTECTED] wrote: > >> Hi everyone > >> > >> I am trying to summarize a table with yield estimates of a forest > >> plantation. For that I have four blocks and four treatments measured > >> over > >> a period of 10 years (every year). In each plot trees are measured > >> (diameters and heights). > >> > >> With aggregate function I can calculate the average diameter or the > >> total > >> volume for each plot at any time with something like this: > >> > >> my_agg <- aggregate(list(Ht = HTO, Dbh = DBH), list(Age = AGE, Block = > >> BLOCK, Treat = TREAT), sum) > >> > >> my_agg <- aggregate(list(VOL = VOLUME), list(Age = AGE, Block = BLOCK, > >> Treat = TREAT), sum) > >> > >> where HTO is the height of the tree > >> Ht average height of trees in each plot at time t > >> DBH is the diameter at 1.3 meter height > >> Dbh is the average diameter of the plot at time t > >> VOL is the [i]th tree volume > >> Vol is the volume of the plot at time t > >> > >> To do actual growth analysis I need to calculate the difference in > >> volume > >> between t and (t-1) for every block/treatment. > >> > >> Now the question > >> > >> How can I use aggregate to calculate the difference of each time series > >> and have an output like aggregate does? I know it is imposible to use > >> diff() within aggregate because it gives me an scalar. > >> In SAS proc means handles such a thing, and diff() in R is fine for a > >> single series, but in this case I want to calculate for all the > >> treatments/blocks at the same time. > >> > >> I can do the whole using by and some more less elegant procedures, but I > >> figure there should be a cleaner way as in PROC MEANS. > >> > >> Any suggestions? > >> > >> Cristian Montes > >> NC 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. > > > > -- > > Andrew Robinson > > Department of Mathematics and StatisticsTel: +61-3-8344-9763 > > University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 > > http://www.ms.unimelb.edu.au/~andrewpr > > http://blogs.mbs.edu/fishing-in-the-bay/ > > -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 in error of mixed models
Hi Bernardo,
it may sound simple-minded, but it looks to me as though you have a
problem with the model. Unfortunately it's difficult for us to
diagnose the problem because you didn't send an example that we can
reproduce. We don't have the "base" data.
Based on my understanding of your code, you might try
lme(logit ~ dis*Modality + non.dis*Modality,
random = ~ dis+non.dis | id,
data = data.bi)
If this doesn't work, try fitting a simpler model, eg
lme(logit ~ dis*Modality + non.dis*Modality,
random = ~ | id,
data = data.bi)
and then using update() to increment towards the full model of
interest. You can also change the model-fitting engine, I have
sometimes had success using Nelder-Mead where nlminb failed. For more
information about this see
?lmeControl
Can I also recommend that you use more spaces in your code?
Good luck
Andrew
On Sun, Nov 04, 2007 at 09:18:09AM -0200, Bernardo Rangel Tura wrote:
>
> Hi R-masters
>
> I read the article: Bivariate analysis of sensitivity and specificity
> produces informative summary measures in diagnostic reviews.
>
> In this paper i proposed a bivariate mixed model and use SAS proc mixed
> to adjust the estimates.
>
>
> I thinks use R to make the same and try with this code:
>
> base<-read.csv("base.csv")
> adj<-.5
> attach(base)
>
> sens<-(VP+adj)/(VP+FN+2*adj)
> log.S<-log(sens/(1-sens))
> var.log.S<-1/(sens*(1-sens)*(VP+FN))
> dis<-rep(1,length(log.S))
> non.dis<-rep(0,length(log.S))
> data.S<-data.frame(id,Modality,log.S,var.log.S,dis,non.dis)
> names(data.S)<-c("id","Modality","logit","var.logit","dis","non.dis")
>
> esp<-(VN+adj)/(VN+FP+2*adj)
> log.E<-log((1-esp)/esp)
> var.log.E<-1/(esp*(1-esp)*(VN+FP))
> dis<-rep(0,length(log.E))
> non.dis<-rep(1,length(log.E))
> data.E<-data.frame(id,Modality,log.E,var.log.E,dis,non.dis)
> names(data.E)<-c("id","Modality","logit","var.logit","dis","non.dis")
>
> data.bi<-rbind(data.S,data.E)
> require(nlme)
> lme(logit~dis*Modality+non.dis*Modality, random=~dis|id+non.dis|
> id,data=data.bi)
>
> but i recive a erro msg :
>
> Error in MEEM(object, conLin, control$niterEM) :
> Singularity in backsolve at level 0, block 1
>
>
> How in solve this problem? Whats is wrong?
>
> Thanks in advance
> --
> Bernardo Rangel Tura, M.D,Ph.D
> National Institute of Cardiology
> Brazil
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] interpreting an LME regression result...
Hi Johan, it's not clear to me that we have enough information to answer your question reliably, but the following thoughts might be useful to you. 1) in your lme() model, all the coefficients are closer to zero than in your lm() model. 2) the change in Z is not very large relative to its standard error under either model. 3) one possible explanation for the change is that the contributions of individual observations to the parameter estimates are weighted differently under lme() than under lm(), and this can be noticeable when the data are unbalanced. Specifically, there might be a large group of individuals who are influencing the parameter estimates strongly in lm(), but who are all from the same family. As individuals, their contribution will lessen in lme(). If your data are unbalanced, then I would look here for the effect. 4) there might be other explanations for the change. I hope that this helps, Andrew On Sat, Nov 10, 2007 at 12:58:21AM -0700, Johan Jackson wrote: > Any help would be most appreciated. (Don't make me get down on my > hands and knees and beg for help, cause I'll do it!!) My boss has me > learning R and doing nested regression with the report due Mon (Friday > night statistics...fun. ). Anyway, here's my problem: > > In a regression equation not accounting for the fact that people are > nested in families, the result for Z variable is VERY strong (beta = > -4511), but this result is much weaker when I use lme and account for > people nested in families (beta = -2613). I'm struggling with a verbal > interpretation of this result. Is it because the effect of Z within > families is not very strong, but between family variation is high? > > > > > - > > If my R syntax would be helpful, it's below: > > BIG EFFECT OF Z (using lm) > summary(lm(Y ~ X + Z + age ,data=pharma.data ,na.action='na.omit')) > > Estimate Std. Error t value Pr(>|t|) > (Intercept) 102.27 24.984.097e-05 *** > X -629.76 267.87 -2.350.020 * > Z-4511.962032.39 -2.220.028 * > age 1.88 1.421.320.188 > > > SMALLER EFFECT OF Z (using lme) > summary(lme(Y ~ X + Z + age, random = ~1|Family.ID, > data=pharma.data, method="ML", na.action='na.omit')) > > Value Std.Error DF t-value p-value > (Intercept) 103.9 20.0 85 5.200 0. > X. -417.3 179.3 85 -2.327 0.0223 > Z -2613.01845.1 85 -1.416 0.1604 > age 1.3 1.2 85 1.126 0.2632 > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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: what are the basis functions in {mgcv}: gam?
Try
Wood S.N. (2006) Generalized Additive Models: An Introduction
with R. Chapman and Hall/CRC Press.
listed in the references in the help file of the function.
It's a great read.
Andrew
On Mon, Mar 23, 2009 at 07:36:44PM -0700, oliviax wrote:
>
> I am writing my thesis with the function gam(), with the package {mgcv}.
>
> My command is: gam(y~s(x1,bs="cr")+s(x2, bs="cr")).
>
> I need help to know what are the default basis funcitons for gam. I have not
> found any detailed reference for this.
>
> Can anyone help me with this??
> --
> View this message in context:
> http://www.nabble.com/help%3A-what-are-the-basis-functions-in-%7Bmgcv%7D%3A-gam--tp22673295p22673295.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.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia (prefer email)
http://www.ms.unimelb.edu.au/~andrewpr Fax: +61-3-8344-4599
http://blogs.mbs.edu/fishing-in-the-bay/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
