# Estimating percentages as the dependent variable in regression

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?

-
Regular (OLS) regression does not make assumptions about the distribution of the dependent variable, it makes assumptions about the distribution of the error (as estimated by residuals) from a model. If the number of students per exam varies, you can probably start with OLS regression and check the assumptions. –  Peter Flom Aug 10 '12 at 14:28
That's a good suggestion, @Peter, but I am concerned about subtler yet important violations of OLS assumptions. A student's rank in an exam will depend on the attributes of all the other students taking the exam. This interdependence is not captured by any (routine) application of OLS. –  whuber Aug 10 '12 at 15:12
Another way to appreciate what's going on here is to consider the simplest instance of this problem, where each "exam" involves exactly two of the students. The dependent variable indicates which student did better. If we think of each exam as a contest, this is equivalent to holding a tournament. The question, in effect, wants to come up with a formula for the "strength" of each player in terms of some explanatory variables. –  whuber Aug 10 '12 at 15:22
You are right @whuber. How to address this? I think an IRT model would do it, but it's been a long time since I studied those. –  Peter Flom Aug 10 '12 at 16:47
@whuber I think that IRT could be a possibility, a ordered probit (or logit) for example. But it's also possible to consider the percentage rank as a bounded continous variable(like the logistic regression suggested by Mike), this would be a valid approximation, whereas the ordered probit would be an straight model to the problem. The advantage of the logistic is parsimony; to account for students correlation, a White or sandwich variance estimator could be implemented. I understood that the sample has 38 exams for several students(longitudinal) so this should be possible. –  JDav Aug 10 '12 at 17:09
show 3 more comments

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."

-
That's nice, and in other contexts would be great advice, but using a glm is not going to fix the problems identified in my answer. –  whuber Aug 13 '12 at 1:51

### Synopsis

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

### Introduction

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.

### Analysis

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.

### Simulation

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

-
I agree with your claims regarding $\varepsilon_i$ (which in answer corresponds to $v_i$) and I hope I got the complete picture of your analysis. Outcomes on $p_{ij}$ would depend on unobservables characteristics, which may be correlated with the observables, like your SAT scores (and others). But isn't trivial to remove the unobservable individual effect by applying a "within" transformation ($y_{ij}$-mean($y_{ij}$)) or any other (difference) ? (+1) –  JDav Aug 12 '12 at 23:07
If your simulation preserves the positive correlation between unobserved abilities and SAT and other scores then the OLS paramters may be upward biased right ?(the asymp. bias is proportional to the controls-random term correlation) so this could explain your "negtive" parameters are atracted towards zero ? –  JDav Aug 12 '12 at 23:28
+1, I'm a big fan of using simulations to help work through statistical ideas. –  gung Aug 13 '12 at 15:23

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%.

-
Logistic regression applies only to binary response variables, which does not appear to be the case here. Are you perhaps suggesting ordinary regression of the logit of the rank percentage? –  whuber Aug 10 '12 at 15:09
Neither logistic regression nor ordinary linear regression is appropriate here for the reasons whuber gave. Also if the model is to have Gaussian residuals I don't see how that is going to translate into a near uniform distribution for the response. The dependence between individual students in the rankings is what I think is most crucial to the analysis. –  Michael Chernick Aug 10 '12 at 15:34
I agree with Mike. Could this be another nomenclature issue ? logistic regression is also refered to the following model: $\ln(p/1-p)=\beta'x_i + u_i$ where $p_i$ is bounded between 0 and 1. –  JDav Aug 10 '12 at 16:44
@J Dav I am unable to find any reference that uses "logistic regression" in the sense you describe. –  whuber Aug 10 '12 at 17:08
I think he mean to do a Generalized Linear Model, with a conditional Gaussian distribution, with a logit-link to the mean. (At least I'm going to assume that's what he meant.) –  Shea Parkes Aug 13 '12 at 2:45
show 2 more comments

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.

-
The idea of using a logistic transformation is a good one, but this kind of regression model has problems (as explained in my answer). –  whuber Aug 12 '12 at 21:17
Just complemented the answer with the estimation method under correlation between the unobserved abilities and observables within $\mu_{ij}$ following @whuber remarks –  JDav Aug 12 '12 at 23:20

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."

-