Unexpected Behavior (potentially) in t.test

Greetings,

I have stumbled across some unexpected behavior (potential a bug) in, 

suspect to be R's (2.6.2 on Ubuntu Linux) t.test function; then again 

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 

the abnormal termination error message "(stderr < 10 * 

max(abs(mx), abs(my))) stop("data are essentially constant")" but only 

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 

sure where I'm going wrong. For background, the goal of this code is to 

a simulation to explore the problem space of inflation of Type I error 

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 

the problem space given that there is no true difference between the 

(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) 

if p is < .25 on the t-test I add t2N more subjects to each group. Then 

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 

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 

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 

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=""))} 

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]]

Hi

I tried your code with R 2.8.0 devel,  gave up after about 4 hours (more 
or less t1N was 18) but until i interrupted it manually there was no 
problem neither error. Maybe you could try it with new R version. Besides 
I did not debug or profile your code so maybe you could try to debug it 
yourself. Especially 3 nested loops seems to me quite ineffective together 
with incremental Cgroup and Tgroup increasing.

  

Regards

Petr
petr.pikal at precheza.cz
724008364, 581252140, 581252257


r-help-bounces at r-project.org 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]]