Thank you Peter and Petr,
I tried it again eliminating some of the t-tests I don't really use (e.g.
the transformed ones), and it hasn't crashed yet. I know that looping (esp
nested looping) in R isn't efficient, but I was having trouble coming up
with alternatives, if you have any suggestions or could point me towards a
resource that explains how to avoid doing loops in what is typically a
looped programing task I'd appreciate it. Once time provides itself I'll
switch up to the new version of R and try again, this time using try(), and
see what happens.
Thanks again for your time and effort,
---
Russell Pierce
Psychology Department
Graduate Student - Cognitive
(951) 827-2553
University of California, Riverside, 92521
>
>
> So you (Russell, not Petr) need to do some more work. You have the
> problem, and you can do things like printing the data vector(s) that
> t.test complains about. Notice that try() allows you to catch the error
> and do something before exiting.
>
> -p
>
> > Regards
> >
> > Petr
> > [EMAIL PROTECTED]
> > 724008364, 581252140, 581252257
> >
> >
> > [EMAIL PROTECTED] napsal dne 20.06.2008 20:55:08:
> >
> >
> >> Greetings,
> >>
> >> I have stumbled across some unexpected behavior (potential a bug) in,
> >>
> > what I
> >
> >> suspect to be R's (2.6.2 on Ubuntu Linux) t.test function; then again
> >>
> > the
> >
> >> problem may exist in my code. I have shutdown R and started it back up,
> >> re-run the code and re-experienced the error. I have searched on Google
> >>
> > for
> >
> >> the abnormal termination error message "(stderr < 10 *
> >>
> > .Machine$double.eps *
> >
> >> max(abs(mx), abs(my))) stop("data are essentially constant")" but only
> >>
> > found
> >
> >> one instance,
> http://tolstoy.newcastle.edu.au/R/e2/help/07/06/18179.html
> >>
> > ,
> >
> >> but the discussion there did not seem particularly helpful.
> >>
> >> I've included all of my code, amateurish though it may be. I have not
> >> isolated the faulty part, and to me it all looks pretty simple, so I'm
> >>
> > not
> >
> >> sure where I'm going wrong. For background, the goal of this code is to
> >>
> > run
> >
> >> a simulation to explore the problem space of inflation of Type I error
> >>
> > when
> >
> >> decisions to run or not to run more participants are made by preliminary
> >> looks at the data (as in Wagenmakers, 2007). This code is meant to
> >>
> > examine
> >
> >> the problem space given that there is no true difference between the
> >>
> > groups
> >
> >> (as is the case when both a generated from random draws from the normal
> >> distribution). I run an initial number of subjects in two groups (t1N)
> >>
> > then,
> >
> >> if p is < .25 on the t-test I add t2N more subjects to each group. Then
> >>
> > I
> >
> >> perform the t.test again. If the p was > .25 at time 1 I stop. Plainp is
> >> simply storing the p-values from t2 (if it was performed) or from t1 (if
> >>
> > t2
> >
> >> was not performed). In the code I provide t1 starts at 16 since this is
> >> about when the problem becomes more frequent. Please note that it takes
> >> quite a long while to fail, and depending on what the true cause is it
> >>
> > may
> >
> >> not fail at all. On my system it is failing before t1N advances to 17.
> >>
> >> Any suggestions as to how to avoid the error and instructions as to the
> >> cause of it would be appreciated. Thank you for your input and patience.
> >>
> >> logit <- function(p)
> >>
> >> {
> >>
> >> # compute and return logit of p;
> >>
> >> # if p=.5 then logit==0 else sign(logit)==(p>.5)
> >>
> >> return( log(p/(1-p)) )
> >>
> >> }
> >>
> >>
> >> antilogit <- function(x)
> >>
> >> {
> >>
> >> # compute and return antilogit of x;
> >>
> >> # this returns a proportion p for which logit(p)==x;
> >>
> >> return( exp(x)/(1+exp(x)) )
> >>
> >> }
> >>
> >>
> >> plainp <- c() #Clear the plainp value
> >>
> >> t1Nsim <- (100/5) * 1000 * 10 # random chance should provide 10000 cases
> >>
> > at
> >
> >> t1
> >>
> >> contthreshold <- .25 #p value below which we run more subjects
> >>
> >> t1pvals <- rep(NA,t1Nsim) #clear the pvalues
> >>
> >> t2pvals <- rep(NA,t1Nsim) #clear the pvalues
> >>
> >> t1N <- 10 #for debugging
> >>
> >> t2N <- 5 #for debugging
> >>
> >>
> >> for (t1N in 16:50) #Outer loop testing possible values for t1N
> >>
> >> for (t2N in 1:50) #Inner loop testing possible values for t2N
> >>
> >> {
> >>
> >> print(paste("Checking with ",t1N," initial samples and ",t2N," extra
> >> samples",sep="")) #feedback
> >>
> >> for (lcv in 1:t1Nsim) #Run simulation t1Nsim times...
> >>
> >> {
> >>
> >> if (lcv %% 20000 == 0) {print(paste((lcv/t1Nsim)*100,"%",sep=""))}
> >>
> > #feedback
> >
> >> Cgroup <- rnorm(t1N) #Initial random draw for Group1
> >>
> >> Tgroup <- rnorm(t1N) #Initial random draw for Group 2
> >>
> >> currentp <- t.test(Cgroup,Tgroup)[["p.value"]] #Get t1 p value
> >>
> >> t1pvals[lcv] <- currentp #Store t1 p value
> >>
> >> #If p >= .05 or <= continue threshold then run more subjects
> >>
> >> if ((currentp <= contthreshold) & (currentp >= .05)) {
> >>
> >> Cgroup <- c(Cgroup,rnorm(t2N)) #Add t2N subjects to group 1
> >>
> >> Tgroup <- c(Tgroup,rnorm(t2N)) #Add t2N subjects to group 2
> >>
> >> currentp <- t.test(Cgroup,Tgroup)[["p.value"]] #Get t2 p value
> >>
> >> t2pvals[lcv] <- currentp #store t2 p value
> >>
> >> }
> >>
> >> }
> >>
> >> plainp <- ifelse(!is.na(t2pvals),{t2pvals},{t1pvals}) #Make sure we are
> >> looking at the right ps
> >>
> >> table(t1pvals <= .05); round(summary(t1pvals),4) #debugging
> >>
> >> hist(t1pvals) #debugging
> >>
> >> table(plainp <= .05); round(summary(plainp),4) #debugging
> >>
> >> hist(plainp, probability=TRUE,main=paste(t1N,"then",t2N)) #Histogram of
> >> interest
> >>
> >> abline(a=1.00,b=0) #Baseline probability
> >>
> >> dev.copy(jpeg,filename=paste("Sim with ",t1N, "start samples and ",t2N,"
> >> extra samples.jpg",sep=""),height=600,width=800,bg="white") # Create the
> >> image
> >>
> >> dev.off() #Save the image
> >>
> >> chi <- rbind(table(t1pvals <= .05),table(plainp <= .05)) #debugging
> >>
> >> chisq.test(chi) #debugging
> >>
> >> explore <- data.frame(t1=t1pvals,t2=t2pvals,picked=plainp) #debugging
> >>
> >> t.test(explore$picked,explore$t1) #debugging
> >>
> >> t.test(logit(explore$picked),logit(explore$t1)) #debugging
> >>
> >> }
> >>
> >> ---
> >> Russell Pierce
> >> Psychology Department
> >> Graduate Student - Cognitive
> >> (951) 827-2553
> >> University of California, Riverside, 92521
> >>
> >> [[alternative HTML version deleted]]
> >>
> >> ______________________________________________
> >> R-help@r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >>
> > http://www.R-project.org/posting-guide.html
> >
> >> and provide commented, minimal, self-contained, reproducible code.
> >>
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
>
> --
> O__ ---- Peter Dalgaard Øster Farimagsgade 5, Entr.B
> c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
> (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
> ~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
>
>
>
[[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.
