I have the rank percentages of students in 38 exams as the dependent variable in my study. A rank percentage is calculated by (a student's rank / number of students in an exam). This dependent variable has nearly uniform distribution and I want to estimate the effects of some variables on the dependent variable.
Which regression approach do I use?
If you are working with Stata have a look at the following example: http://www.ats.ucla.edu/stat/stata/faq/proportion.htm
Here is a quote from this webpage:
"How does one do regression when the dependent variable is a proportion?
Proportion data has values that fall between zero and one. Naturally, it would be nice to have the predicted values also fall between zero and one. One way to accomplish this is to use a generalized linear model (glm) with a logit link and the binomial family. We will include the robust option in the glm model to obtain robust standard errors which will be particularly useful if we have misspecified the distribution family."
Regression results may have some limited value when carefully interpreted. Unavoidable forms of variation will cause coefficient estimates to shrink substantially towards zero. A better model is needed that handles variation in a more appropriate way.
(A maximum likelihood model can be constructed but may be impracticable due to the computation needed, which involves numerical evaluation of multidimensional integrals. The numbers of dimensions are equal to the numbers of students enrolled in the classes.)
As a narrative to inform our intuition, imagine that these 38 exams were given in 38 separate courses during one semester at a small school with enrollment of 200 college students. In a realistic situation those students will have varying abilities and experiences. As surrogate measures of these abilities and experiences we might take, say, scores on the SAT math and verbal tests and year in college (1 through 4).
Typically, students will enroll in courses according to their abilities and interests. Freshmen take introductory courses and introductory courses are populated primarily by freshmen. Upperclassmen and talented freshmen and sophomores take the advanced and graduate-level courses. This selection partially stratifies the students so that the innate abilities of students within any class are typically more homogeneous than the spread of abilities throughout the school.
Thus, the most capable students may find themselves scoring near the bottom of the difficult, advanced classes in which they enroll, while the least capable students may score near the top of the easy introductory classes they take. This may confound a direct attempt to relate the exam ranks directly to attributes of students and the classes.
Index the students with $i$ and let the attributes of student $i$ be given by the vector $\mathbf{x}_i$. Index the classes with $j$ and let the attributes of class $j$ be given by the vector $\mathbf{z}_j$. The set of students enrolled in class $j$ is $A_j$.
Assume the "strength" of each student $s_i$ is a function of their attributes plus some random value, which may as well have zero mean:
$$s_i = f(\mathbf{x}_i, \beta) + \varepsilon_i.$$
We model the exam in class $j$ by adding independent random values to the strength of each student enrolled in the class and converting those to ranks. Whence, if student $i$ is enrolled in class $j$, their relative rank $r_{i,j}$ is determined by their position in the sorted array of values
$$\left(s_k + \delta_{k,j}, k \in A_j\right).$$
This position $r_{i,j}$ is divided by one more than the total class enrolment to give the dependent variable, the percentage rank:
$$p_{i,j} = \frac{r_{i,j}}{1 + |A_j|}.$$
I claim that the regression results depend (quite a bit) on the sizes and structure of the random (unobserved) values $\varepsilon_i$ and $\delta_{i,j}$. The results also depend on precisely how students are enrolled in classes. This should be intuitively obvious, but what is not so obvious--and appears difficult to analyze theoretically--is how and how much the unobserved values and the class structures affect the regression.
Without too much effort we can simulate this situation to create and analyze some sample data. One advantage of the simulation is that it can incorporate the true strengths of the students, which in reality are not observable. Another is that we can vary the typical sizes of the unobserved values as well as the class assignments. This provides a "sandbox" for assessing proposed analytical methods such as regression.
To get started, let's set the random number generator for reproducible results and specify the size of the problem. I use R because it's available to anyone.
set.seed(17)
n.pop <- 200 # Number of students
n.classes <- 38 # Number of classes
courseload <- 4.5 # Expected number of classes per student
To provide realism, create n.classes classes of varying difficulties on two scales (mathematical and verbal, with a negative correlation), conducted at varying academic levels (ranging from 1=introductory to 7=research), and with variable ease. (In an "easy" class, differences among the amounts of student learning may be large and/or the exam may provide little discrimination among the students. This is modeled by random terms $\delta_{i,j}$ that, for class $j$ tend to be large. The exam results will then be almost unpredictable from the student strength data. When the class is not "easy," these random terms are negligibly small and the student strengths can perfectly determine the exam rankings.)
classes <- data.frame(cbind(
math <- runif(n.classes),
rbeta(n.classes, shape1=(verbal <- (1-math)*5), shape2=5-verbal),
runif(n.classes, min=0, max=7),
rgamma(n.classes, 10, 10)))
rm(math, verbal)
colnames(classes) <- c("math.dif", "verbal.dif", "level", "ease")
classes <- classes[order(classes$math.dif + classes$verbal.dif + classes$level), ]
row.names(classes) <- 1:n.classes
plot(classes, main="Classes")
The students are spread among the four years and endowed with random values of their attributes. There are no correlations among any of these attributes:
students <- data.frame(cbind(
as.factor(ceiling(runif(n.pop, max=4))),
sapply(rnorm(n.pop, mean=60, sd=10), function(x) 10*median(c(20, 80, floor(x)))),
sapply(rnorm(n.pop, mean=55, sd=10), function(x) 10*median(c(00, 80, floor(x)))),
rnorm(n.pop)
))
colnames(students) <- c("year", "math", "verbal", "ability")
plot(students, main="Students")
The model is that each student has an inherent "strength" determined partly by their attributes and partly by their "ability," which is the $\varepsilon_i$ value. The strength coefficients beta, which determine the strength in terms of other attributes, are what the subsequent data analysis will seek to estimate. If you want to play with this simulation, do so by changing beta. The following is an interesting and realistic set of coefficients reflecting continued student learning throughout college (with a large amount between years 2 and 3); where 100 points on each part of the SAT are worth about one year of school; and where about half the variation is due to the "ability" values not captured by SAT scores or year in school.
beta <- list(year.1=0, year.2=1, year.3=3, year.4=4, math=1/100, verbal=1/100, ability=2, sigma=0.01)
students$strength <- (students$year==1)*beta$year.1 +
(students$year==2)*beta$year.2 +
(students$year==3)*beta$year.3 +
(students$year==4)*beta$year.4 +
students$math*beta$math +
students$verbal*beta$verbal +
students$ability*beta$ability
students <- students[order(students$strength), ]
row.names(students) <- 1:n.pop
(Bear in mind that students$ability is unobservable: it is an apparently random deviation between the strength predicted from the other observable attributes and the actual strength on exams. To remove this random effect, set beta$ability to zero. beta$sigma will multiply the ease values: it's basically the standard deviation of the $\delta_{i,j}$ relative to the range of strengths of students in a given course. Values around $.01$ to $.2$ or so seem reasonable to me.)
Let the students pick courses to match their abilities. Once they do that, we can compute the class sizes and stash those with the classes dataframe for later use. The value of spread in the assignments <-... line determines how closely the students are sectioned into classes by ability. A value close to $0$ essentially pairs the weakest students with the easiest courses. A value close to the number of classes spreads the students out a little more. Much larger values than that start to get unrealistic, because they tend to put weaker students into the most difficult courses.
pick.classes <- function(i, k, spread) {
# i is student strength rank
# k is number to pick
p <- pmin(0.05, diff(pbeta(0:n.classes/n.classes, i/spread, (1+n.pop-i)/spread)))
sample(1:n.classes, k, prob=p)
}
students$n.classes <- floor(1/2 + 2 * rbeta(n.pop,10,10) * courseload)
assignments <- lapply(1:n.pop, function(i) pick.classes(i, students$n.classes[i], spread=1))
enrolment <- function(k) length(seq(1, n.pop)[sapply(assignments, function(x) !is.na(match(k, x)))])
classes$size <- sapply(1:n.classes, enrolment)
classes$variation <- by(data, data$Class, function(x) diff(range(x$strength)))
(As an example of what this step has accomplished, see the figure further below.)
Now apply the model: the abilities of the students in each class are independently varied--more for easy exams, less for hard (discriminating) exams--to determine their exam scores. These are summarized as ranks and "pranks", which are rank percents. The pranks for a class of $n$ students range from $1/(n+1)$ through $n/(n+1)$ in increments of $1/(n+1)$. This will later make it possible to apply transformations such as the logistic function (which is undefined when applied to values of $0$ or $1$).
exam.do <- function(k) {
s <- seq(1, n.pop)[sapply(assignments, function(x) !is.na(match(k, x)))]
e <- classes$ease[k]
rv <- cbind(rep(k, length(s)), s, order(rnorm(length(s), students$strength[s], sd=e*beta$sigma*classes$variation[k])))
rv <- cbind(rv, rv[,3] / (length(s)+1))
dimnames(rv) <- list(NULL, c("Class", "Student", "Rank", "Prank"))
rv
}
data.raw <- do.call(rbind, sapply(1:n.classes, exam.do))
To these raw data we attach the student and class attributes to create a dataset suitable for analysis:
data <- merge(data.raw, classes, by.x="Class", by.y="row.names")
data <- merge(data, students, by.x="Student", by.y="row.names")
Let's orient ourselves by inspecting a random sample of the data:
> data[sort(sample(1:dim(data)[1], 5)),]
Row Student Class Rank Prank math.dif verbal.dif level ease Size year math verbal ability strength n.classes
118 28 1 22 0.957 0.77997 6.95e-02 0.0523 1.032 22 2 590 380 0.576 16.9 4
248 55 5 24 0.889 0.96838 1.32e-07 0.5217 0.956 26 3 460 520 -2.163 19.0 5
278 62 6 22 0.917 0.15505 9.54e-01 0.4112 0.497 23 2 640 510 -0.673 19.7 4
400 89 10 16 0.800 0.00227 1.00e+00 1.3880 0.579 19 1 800 350 0.598 21.6 5
806 182 35 18 0.692 0.88116 5.44e-02 6.1747 0.800 25 4 610 580 0.776 30.7 4
Record 118, for example, says that student #28 enrolled in class #1 and scored 22nd (from the bottom) on the exam for a percentage rank of 0.957. This class's overall level of difficulty was 0.0523 (very easy). A total of 22 students were enrolled. This student is a sophomore (year 2) with 590 math, 380 verbal SAT scores. Their overall inherent academic strength is 16.9. They were enrolled in four classes at the time.
This dataset comports with the description in the question. For instance, the percentage ranks indeed are almost uniform (as they must be for any complete dataset, because the percentage ranks for a single class have a discrete uniform distribution).
Remember, by virtue of the coefficients in beta, this model has assumed a strong connection between examination scores and the variables shown in this dataset. But what does regression show? Let's regress the logistic of the percentage rank against all the observable student characteristics that might be related to their abilities, as well as the indicators of class difficulty:
logistic <- function(p) log(p / (1-p))
fit <- lm(logistic(Prank) ~ as.factor(year) + math + verbal + level, data=data)
summary(fit)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.577788 0.421579 -6.11 1.5e-09 ***
as.factor(year)2 0.467846 0.150670 3.11 0.0020 **
as.factor(year)3 0.984671 0.164614 5.98 3.2e-09 ***
as.factor(year)4 1.109897 0.171704 6.46 1.7e-10 ***
math 0.002599 0.000538 4.83 1.6e-06 ***
verbal 0.002130 0.000514 4.14 3.8e-05 ***
level -0.208495 0.036365 -5.73 1.4e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.48 on 883 degrees of freedom
Multiple R-squared: 0.0661, Adjusted R-squared: 0.0598
F-statistic: 10.4 on 6 and 883 DF, p-value: 3.51e-11
Diagnostic plots (plot(fit)) look fastastic: the residuals are homoscedastic and beautifully normal (albeit slightly short tailed, which is no problem); no outliers; and no untoward influence in any observation.
As you can see, everything is highly significant, although the small R-squared might be disappointing. The coefficients all have the roughly the correct signs and relative sizes. If we were to multiply them by $3.5$, they would equal $(-9, 1.6, 3.4, 3.9, 0.009, 0.007, -0.7)$. The original betas were $(*, 1, 3, 4, 0.010, 0.010, *)$ (where $*$ stands for a coefficient that was not explicitly specified).
Notice the high significance of level, which is an attribute of the classes, not of the students. Its size is pretty large: the class levels range from near $0$ to near $7$, so multiplying this range by the estimated coefficient of level show it has the same size of effect as any of the other terms. Its negative sign reflects a tendency for students to do a little bit worse in the more challenging classes. It is very interesting to see this behavior emerge from the model, because the level was never explicitly involved in determining the examination outcomes: it only affected how the students chose their classes.
(By the way, using the percentage ranks untransformed in the regression does not qualitatively change the results reported below.)
Let's vary things a bit. Instead of setting spread to $1$, we were to use $38$, thereby causing a wider (more realistic) distribution of students throughout the classes. Rerunning everything from the top gives these results:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.902006 0.349924 -14.01 < 2e-16 ***
as.factor(year)2 0.605444 0.130355 4.64 3.9e-06 ***
as.factor(year)3 1.707590 0.134649 12.68 < 2e-16 ***
as.factor(year)4 1.926272 0.136595 14.10 < 2e-16 ***
math 0.004667 0.000448 10.41 < 2e-16 ***
verbal 0.004019 0.000434 9.25 < 2e-16 ***
level -0.299475 0.026415 -11.34 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.3 on 883 degrees of freedom
Multiple R-squared: 0.282, Adjusted R-squared: 0.277
F-statistic: 57.9 on 6 and 883 DF, p-value: <2e-16
(In this scatterplot of class assignments, with spread set to $38$, students are sorted by increasing strength and classes are sorted by increasing level. When spread originally was set to 1, the assignment plot fell in a tight diagonal band. Weaker students tend to take easier classes and stronger students take harder classes, but there are plenty of exceptions.)
This time the R-squared is much improved (although still not great). However, all the coefficients have increased by 20 - 100%. This table compares them along with some additional simulations:
Simulation Intercept Year.2 Year.3 Year.4 Math Verbal Level R^2
Beta * 1.0 3.0 4.0 .010 .010 * *
Spread=1 -2.6 0.5 1.0 1.1 .003 .002 -0.21 7%
Spread=38 -4.9 0.6 1.7 1.9 .005 .004 -0.30 25%
Ability=1 -8.3 0.9 2.6 3.3 .008 .008 -0.63 58%
No error -11.2 1.1 3.3 4.4 .011 .011 -0.09 88%
Keeping spread at $38$ and changing ability from $2$ to $1$ (which is a very optimistic assessment of how predictable the student strengths are) yielded the penultimate line. Now the estimates (for student year and student SAT scores) are getting reasonably close to the true values. Finally, setting both ability and sigma to $0$, to remove the error terms $\varepsilon_i$ and $\delta_{i,j}$ altogether, gives a high R squared and produces estimates close to the correct values. (It is noteworthy that the coefficient for level then decreases by an order of magnitude.)
This quick analysis shows that regression, at least as performed here, is going to confound unavoidable forms of variation with the coefficients. Furthermore, the coefficients also depend (to some extent) on how students are distributed among classes. This can partially be accommodated by including class attributes among the independent variables in the regression, as done here, but even so the effect of student distribution does not disappear.
Any lack of predictability of true student performance, and any variation in student learning and actual performance on examinations, apparently cause the coefficient estimates to shrink towards zero. They appear to do so uniformly, suggesting that the relative coefficients may still be meaningful.
The measure @user13203 proposes may be considered as a continous bounded underperformance score, the lower the better the performance: $y_{ij} $ i-th student underperformance at j-th exam.
Using a logit linearizing transformation where $\mu_{ij}$ may depend on observable student or exams characteristics :
$\ln(y_{ij}/(1-y_{ij})) = \mu_{ij} + e_{ij} + v_i $
student's unobserved skills are modeled through the random component $v_i$ while $e_{ij}$ models other non systematic unobservables. Correlation between responses (examinations) may be addresed by assming a general covariance structure for $e_{ij}$. Why not a White (or sandwich/robust) variance structure ? Moreover, some of the responses correlation can be accounted within the $\mu_{ij}$ (conditional dependence).
(This is just an idea from my biased experience, comments and critics are more than welcome.)
Unobservable abilities are likely to be correlated with students or exams observables attributes within $\mu_{ij}$. This assumptions makes this model a RE with correlated error components, that can be estimated by ML or a two stage estimator: first stage: a within (or analog) transformation that eliminates $v_i$. Second stage: OLS on the transformed model.
You might want to try logistic regression. The logit transform $\ln(\frac{p}{1-p})$ will spread your response variable out over the real line so you won't get absurd predicted rank percentages like -3% or +110%.
A perfect model in this case will map the inputs (whatever covariates you have) to the outputs (the rank of the student in the class). Another way to think of this is by mapping first to the scores, and then mapping those scores to the rank. I'm going to ignore error for now.
test score: $y = \sum \beta x$
rank: $ r = R(y)$
In which $R$ is the ranking function. The problem is that $R$ is a non-linear function that depends entirely on the data itself. If we assume that we have an infinite amount of data, then we know the complete distribution of $y$, and $R(y)$ is essentially the cumulative density function. It tells you what percent of people scored worse than you on the test, the area to the left of your score.
This appears to be quite similar to the functional form of the generalized linear model. I think this is why the logistic regression approach was proposed by @Mike Anderson. If your exam scores were logistically distributed, then the link function to use would be the logit (its inverse is the cumulative density function we care about). Similarly, if the scores were normally distributed, the probit function would be the link function.
For your regression, the only way to estimate ranks is to say "given that my data are distributed as X, this point is in the 34th percentile". Otherwise, how do you know what a two point increase in your test score translates to in terms of rank? The caveat is that you have to estimate that distribution in order to choose your link function (certain functional forms will make your life a lot easier). Furthermore, this model isn't going to say "you were the 6th best out of a class of 38", rather "if the test scores were distributed how we think they are, your score would put you in the 15th percentile."
