## all-DIF.R -- detects uniform and non-uniform DIF
##  NOTE: you must have run this with RFILE=filename.rf

##  This version has been changed for 2009 data
##     if you need the 2007 survey version, use
##     all-DIF-orig.R

## Set a few parameters

## This is the base path where the program will look for files
##   and save graphs. The full path will be mypath/[st]meas/<measname>/
mypath <- "/home/sl70/surveys/surv2009/"
## Enter measure name here
measname <- "trpa"
## Enter s for student, t for teacher
tors <- "t"

## Enter the codes for the groups to compare here
##  these are take from the values of the vector a below
group1 <- "e"
group2 <- "h"

## If you want to print probability curves, set the following to T; otherwise, F
print.probability.curves <- T

## How many is good here? Will different numbers change our inferences?

number.of.strata <- 12

## Some constants -- comparison codes
a <- c("e", "h",  "s", "c", "u")
## Comparison group titles
b <- c("Elementary", "High School", "Senior", "Combined 6-10", "Upper Grades")
## Student or teacher measure
c <- c("s", "t")
## Student or teacher title
d <- c("Student", "Teacher")



## First we need to read in person, item and step files

## first items
## You might need to change these if you write out the item file in a different format
##   Also, I use HLINES=N in my control file to remove the headers from the output files
##   If you don't do this, you need to add skip=4 to the read.fwf line to skip the headers
if.widths <- c(1, 5, 8, 3, 6, 7, 7, 7, 7, 7, 7, 7, 7, 2, 2, 85)

fe <- paste(mypath, tors, "meas/testing-", measname, "/", measname,  "2.if", sep="")
bar <- read.fwf(file=fe,widths=if.widths)
names(bar) <- c("semi", "entry", "measure", "st", "count", "score",
                "error", "inmsq", "instd", "outms", "outstd", "displc", "ptbis", "group",
                "r", "name")
okrow <- (bar$st >= -1)
items <- bar[okrow,]

## need the number of items later when we read in the response file
number.items <- length(items$entry)

rm("bar")

## then steps
## item.group table has to have:
##   item.group
##   group step.num step.diff

steps <- read.table(paste(mypath, tors, "meas/testing-", measname, "/", measname, ".sf", sep=''), header=F)
if (length(steps[1,]) > 2) {
  names(steps) <- c("item.num", "step.num", "step.diff")
  ## Make item-group correspondence
   item.group <- unique(merge(items[c(2, 14)], steps, by.x="entry", by.y="item.num")[2:4])
} else {
  names(steps) <- c("step.num", "step.diff")
  item.group <- cbind(items[1:length(steps$step.diff), 14], steps)
  names(item.group) <- c("group", "step.num", "step.diff")
} 

## read in people

pf.widths <- c(1, 5, 8, 3, 6, 7, 7, 7, 7, 7, 7, 7, 7, 5, 12, 12)

foo <- read.fwf(file=paste(mypath, tors, "meas/testing-", measname, "/", measname, ".pf", sep=""),
                widths=pf.widths)
names(foo) <- c("semi", "entry", "measure", "st", "count", "score",
                "error", "inmsq", "instd", "outms", "outstd", "displc", "ptbis", "unit",
                "id", "str")
person.order <- order(foo$measure)
all.people <- foo[person.order,]

rm("foo")

min.meas <- floor(all.people$measure[1])
max.meas <- ceiling(all.people$measure[length(all.people$measure)])


## These are the conditions for selecting the groups to compare
##   I select records by order in the data set. If you use different
##   selection criteria, use them here. The conditions should be
##   in the same order as the codes in vector a above.


group.cond <- c(expression(all.people$entry<=2000),                          ### elementary
                expression(all.people$entry<=4000 & all.people$entry>2000),  ### hs
                expression(all.people$entry<=6000 & all.people$entry>4000),  ### seniors
                expression(all.people$entry<=4000) ,                         ### grades 6-10
                expression(all.people$entry>4000)                            ### upper grades
                )

group1.num <- which(a==group1)
group2.num <- which(a==group2)

compare.1 <-  eval(group.cond [ group1.num ])
compare.2 <-  eval(group.cond [ group2.num ])


## Read in response file -- do this instead of recoding

## if there are more than 14 items then the id field in the response file grows,
##   so we have to adjust the field widths
id.length <- ifelse(number.items<15, 30, 16 + number.items)

rf.widths <- c(id.length, number.items)
foor <- read.fwf(file=paste(mypath, tors, "meas/testing-", measname, "/", measname, ".rf", sep=''), widths=rf.widths)
responses <- foor[person.order,]
names(responses) <- c("id", "responses")

rm("foor")

item.string <- as.character(responses$responses)

## let's read these all into a big matrix -- make it easier to manipulate
##   ncol should be number of items; nrow is number of people

item.matrix <- matrix(0, ncol=length(items$entry), nrow=length(all.people$measure))

## Note to self: the item string is length(bar$entry) long,
##   but the item matrix length has length(items$entry) columns.
##   Be careful with this difference

for (ii in 1:length(items$entry)) {
  item.matrix[ ,ii] <-  as.numeric(substring(item.string, items$entry[ii], items$entry[ii])) 
}

item.matrix.compare.1 <- item.matrix[compare.1,]
item.matrix.compare.2 <- item.matrix[compare.2,]




## for all items uncomment out this line -- this will
##   make one graph for each item in the measure
for (k in 1:length(items$measure)) {
## Otherwise, for a single item, uncomment the next line and insert
##  the appropriate number instead of 6  
## for (k in 6) {

  ## so for item k the thresholds will be item.group$step.diff[item.group$group==items[k, 14] ]

   thresh <- item.group$step.diff[ item.group$group==items[k, 14] ]

   ## This calculates the expected scores
   ## This method make one j x n matrix for each of the k items

   alldenom <- matrix(0, ncol=length(thresh), nrow=length(all.people$measure))
   for(i in 1:length(thresh)) {
     alldenom[,i] <- exp(i * (all.people$measure-items$measure[k]) - sum(thresh[1:i]))
   }

   denom <- apply(alldenom, 1, sum)

   numer <- matrix(0, ncol=length(thresh), nrow=length(all.people$measure))

   for(the.cat in 1:length(thresh)) {
     numer[, the.cat] <-  the.cat * exp(the.cat * (all.people$measure-items$measure[k]) - sum(thresh[1:the.cat]))
   }

   expected.score.curve <- apply(numer, 1, sum)/denom

   ## postscript(paste("expected-score-ogive-", measname, ".ps", sep=''), horiz=F, height=8, width=7)
   ## 
   ## plot(all.people$measure, expected.score.curve, type="l",
   ##      ylim=c(1, 4), xlim=c(min.meas, max.meas),
   ##      xlab="Person minus item measure",
   ##      ylab="Expected Score",
   ##      main="Expected Score Ogive: Means")
   ## 
   ## 
   ## dev.off()

   ## this calculates category probabilities -- need this for calculating the score variance later

   for(the.cat in 1:length(thresh)) {
     numer[, the.cat] <-   exp(the.cat * (all.people$measure-items$measure[k]) - sum(thresh[1:the.cat]))
   }

   probability.score.curves <- numer/denom

   ## want to do this only once for each group
   
   if (print.probability.curves) {

     cat.colors <- hcl(h=seq(300, 30, length=length(thresh)), l=70, c=40)



     postscript(paste(mypath, tors, "meas/testing-", measname, "/","category-probabilities-", measname,
                      item.group$group==items[k, 14], ".ps", sep=''), horiz=F, height=8, width=7)

     plot(all.people$measure, probability.score.curves[,1], type="l", col=cat.colors[1],
          ylim=c(0, 1), xlim=c(min.meas, max.meas),
          xlab="Person minus item measure",
          ylab="Probability of response",
          main="Category Probabilities: Modes"
          )
     for (cc in 2:length(thresh)) {
       lines(all.people$measure, probability.score.curves[,cc], type="l", col=cat.colors[cc])
     }

     dev.off()
   }

   ## get vectors of individual item responses; also subsets by level

   this.item <- item.matrix[,k]
   this.item.compare.1 <- item.matrix.compare.1[,k]
   this.item.compare.2 <- item.matrix.compare.2[,k]

   ## Make factors for the strata
   ##  note we use rank() to get strata with approximately the same numbers of people in each.
   ##  Otherwise, the anova might not come out right because of unequal cell sizes
   measure.strata <- ceiling(rank(all.people$measure, ties.method="random")/(length(all.people$measure)/number.of.strata))

   print("Numbers of people in strata")
   print(table(measure.strata))

   ## Calculate mean scores and mean measures for each stratum and
   ##   separately by school level
   mean.score <- tapply(this.item, measure.strata, mean, na.rm=T)
   mean.meas <-  tapply(all.people$measure, measure.strata, mean, na.rm=T)

   mean.meas.compare.1 <-  tapply(all.people$measure[compare.1], measure.strata[compare.1], mean, na.rm=T)
   mean.meas.compare.2 <-  tapply(all.people$measure[compare.2], measure.strata[compare.2], mean, na.rm=T)

   mean.score.compare.1 <- tapply(this.item.compare.1, measure.strata[compare.1], mean, na.rm=T)
   mean.score.compare.2 <- tapply(this.item.compare.2, measure.strata[compare.2], mean, na.rm=T)

   print(paste("Mean Measures and Scores by Strata for", items$name[k], "for", b[group1.num]))
   print(mean.meas.compare.1)
   print(mean.score.compare.1)

   print(paste("Mean Measures and Scores by Strata for", items$name[k], "for", b[group2.num]))
   print(mean.meas.compare.2)
   print(mean.score.compare.2)

   
   ## Calculate score variance, used to calculate the standardized residuals
   ##  Formula is in RSA p. 108, equation 5.7.2:
   ##    \[ W_{ni} = \sum_{k=0}^m (k - E_{ni})^2\pi_{nik} \]

   score.variance <- rep(0, 4000)

   for(cat in 1:length(thresh)) {
     score.variance <- score.variance + (cat - expected.score.curve)^2 * probability.score.curves[,cat]
   }

   std.resid <- (this.item - expected.score.curve)/sqrt(score.variance)

   ## Do anova on residual by measure strata and school level

   aov.results <- aov(std.resid ~ measure.strata * compare.2)

   print(summary(aov.results))

   ## get significance of effects from the output
   compare.2.compare.1.effect.sig <- 1-pf(summary(aov.results)[[1]]$"F value"[2], ### F statistic for the school level effect
                                          summary(aov.results)[[1]]$Df[2],        ### Number of degrees of freedom in numerator
                                          summary(aov.results)[[1]]$Df[4])        ### Number of degrees of freedom in denominator

   interaction.effect.sig <- 1-pf(summary(aov.results)[[1]]$"F value"[3],
                                  summary(aov.results)[[1]]$Df[3],
                                  summary(aov.results)[[1]]$Df[4])                   

   ##  Significance marks: <0.001: ***, <0.01: **, <0.05: *, >=0.05, .
   compare.2.compare.1.effect.sig.mark <- ifelse(compare.2.compare.1.effect.sig<0.001, "***",
                                                 ifelse(compare.2.compare.1.effect.sig<0.01, "**",
                                                        ifelse(compare.2.compare.1.effect.sig<0.05, "*", ".")))
   interaction.effect.sig.mark <- ifelse(interaction.effect.sig<0.001, "***",
                                         ifelse(interaction.effect.sig<0.01, "**",
                                                ifelse(interaction.effect.sig<0.05, "*", ".")))

   postscript(paste(mypath, tors, "meas/testing-", measname, "/", measname, "DIF", k, ".ps", sep=''), horiz=F, height=8, width=7)

   plot(all.people$measure, expected.score.curve, type="l",
        ylim=c(1, length(thresh)), xlim=c(min.meas, max.meas),
        xlab="Person minus item measure",
        ylab="Expected Score",
        main=paste("Expected Score Ogive\n", items$name[k]))

   lines(mean.meas.compare.1, mean.score.compare.1, type="l", col="red")
   lines(mean.meas.compare.2, mean.score.compare.2, type="l", col="blue")
   legend(max.meas-5, 1.7, legend=c(b[group1.num],
                             paste(b[group2.num], compare.2.compare.1.effect.sig.mark),
                             paste("Interaction", interaction.effect.sig.mark)),
          lty=c(1, 1, NA), col=c("red", "blue", NA))

   dev.off()


 }