[R] Help writing a Scheffe Contrast function for R

[peter dalgaard]

On Mar 10, 2011, at 05:49 , Tyler Rinker wrote:

> 
> My Question:  What do I need to do to correct the three error codes R gives me and make the function run correctly?
> 
> This is the session, code and R's error message when supplied with data: 
>> rm(list=ls())
>> dat1<-read.table("dat1.csv", header=TRUE, sep=",",na.strings="999")
>> attach(dat1)
>> dat1
>   student  rti score
> 1        1   ns     2
> 2        2   ns     5
> 3        3   ns     2
> 4        4   ns    11
> 5        5 wk2o    10
> 6        6 wk2o    11
> 7        7 wk2o     7
> 8        8 wk2o    12
> 9        9 wk5o    12
> 10      10 wk5o    12
> 11      11 wk2w     5
> 12      12 wk2w     6
> 13      13 wk2w    12
> 14      14 wk5w     5
> 15      15 wk5w     6
> 16      16 wk5w    13
>> anova(lm(score~rti))
> Analysis of Variance Table
> Response: score
>          Df  Sum Sq Mean Sq F value Pr(>F)
> rti        4  83.771  20.943  1.7107 0.2174
> Residuals 11 134.667  12.242               
>> #so MSw is 12.242
>> 
>> #my code for scheffe's post hoc comparison
>> 
>> scheffe<- function(IV,DV,data,group1,group2,MSw,alpha) {
> +  result<-0
> +  J<-length(levels(IV))
> +  d1<-subset(data, IV == "group1") 
> +  d2<-subset(data, IV == "group2")
> +  g1<-d1$DV
> +  g2<-d2$DV


This is your problem. You seem to think that d1$DV extracts d1$score when score is passed in the DV argument, but it goes looking for a DV component. For a quick fix, pass "score" as a character string and use d1[[DV]]. Also beware that your IV argument is passed by value, so it may not be available, and if available, it may not be the object that you are looking for.


> +  y.1<-mean(g1)
> +  y.2<-mean(g2) 
> +  n1<- length(g1)
> +  n2<- length(g2)
> +  N<-length(IV)
> +  psi.hat<-y.1-y.2
> +  se.psi.hat<- sqrt(MSw*((1/n1)+(1/n2)))
> +  t<- psi.hat/se.psi.hat
> +  t.compare<-abs(t)
> +  k<-sqrt((J-1)*( qf((1-alpha),(J-1),(N-J))))
> +   if(t.compare > k) return<-c("reject H0") else result<- c("accept H0")
> +   return(list(cbind(result)), Mean.Group.1 = y.1,n.for.1=n1,   
> +   Mean.Group.2=y.2,n.for.2=n2,N.Data=N,tValue=t,critical.value=k)
> +  }
>> #
>> #doesn't seem to be a problem thus far
>> #let's enter some data
>> 
>> #scheffe(IV,DV,data,group1,group2,MSw,alpha)
>> 
>> scheffe(rti,score,dat1,ns,wk5o,12.242,.05)
> Error in if (t.compare > k) return <- c("reject H0") else result <- c("accept H0") : 
>  missing value where TRUE/FALSE needed
> In addition: Warning messages:
> 1: In mean.default(g1) : argument is not numeric or logical: returning NA
> 2: In mean.default(g2) : argument is not numeric or logical: returning NA
>> 
> 
> 		 	   		  
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com