Szacowanie odsetków jako zmiennej zależnej w regresji

Mam procentowe stopnie studentów na 38 egzaminach jako zmienną zależną w moim badaniu. Procent rangi jest obliczany na podstawie (rangi studenta / liczby studentów na egzaminie). Ta zmienna zależna ma prawie jednolity rozkład i chcę oszacować wpływ niektórych zmiennych na zmienną zależną.

Jakiego podejścia regresji używam?

Jeśli pracujesz ze Statą, spójrz na następujący przykład: http://www.ats.ucla.edu/stat/stata/faq/proportion.htm

Oto cytat z tej strony:

„Jak wykonuje się regresję, gdy zmienna zależna jest proporcją?

Dane dotyczące proporcji zawierają wartości od zera do jednego. Oczywiście byłoby miło, gdyby przewidywane wartości również zawierały się w przedziale od zera do jednego. Jednym ze sposobów osiągnięcia tego jest użycie uogólnionego modelu liniowego (glm) z łączem logit i rodziną dwumianową. Uwzględnimy solidną opcję w modelu glm, aby uzyskać solidne standardowe błędy, które będą szczególnie przydatne, jeśli źle sprecyzowaliśmy rodzinę dystrybucji ”.

Streszczenie

Wyniki regresji mogą mieć pewną ograniczoną wartość, jeśli zostaną dokładnie zinterpretowane. Nieuniknione formy zmienności spowodują, że oszacowania współczynników znacznie się skurczą w kierunku zera. Potrzebny jest lepszy model, który obsługuje zmiany w bardziej odpowiedni sposób.

(Można zbudować model maksymalnego prawdopodobieństwa, ale może on być niewykonalny ze względu na wymagane obliczenia, które obejmują numeryczną ocenę całek wielowymiarowych. Liczby wymiarów są równe liczbie uczniów zapisanych na zajęcia.)

Wprowadzenie

Jako narrację informującą naszą intuicję, wyobraź sobie, że tych 38 egzaminów odbyło się w 38 osobnych kursach w ciągu jednego semestru w małej szkole z zapisaniem 200 studentów. W realistycznej sytuacji uczniowie będą mieli różne umiejętności i doświadczenia. Jako zastępcze miary tych umiejętności i doświadczeń możemy wziąć, powiedzmy, wyniki z matematyki SAT i testów werbalnych oraz rok na studiach (od 1 do 4).

Zazwyczaj uczniowie zapisują się na kursy zgodnie z ich umiejętnościami i zainteresowaniami. Początkujący rozpoczynają kursy wprowadzające, a kursy wprowadzające są wypełniane głównie przez pierwszoklasistów. Uczniowie wyższych klas i utalentowani pierwszoklasanci i absolwenci uczęszczają na kursy dla zaawansowanych i absolwentów. Ten wybór częściowo rozwarstwia uczniów, tak że wrodzone zdolności uczniów w każdej klasie są zazwyczaj bardziej jednorodne niż rozprzestrzenianie się umiejętności w szkole.

Zatem najzdolniejsi uczniowie mogą zdobyć punkty w dolnej części trudnych, zaawansowanych klas, w których się zapisują, podczas gdy najsłabsi uczniowie mogą zdobyć punkty w górnej części łatwych zajęć wprowadzających, które biorą. Może to mylić bezpośrednią próbę powiązania stopni egzaminu bezpośrednio z atrybutami uczniów i klas.

Analiza

Zindeksuj uczniów za pomocą i pozwól, aby atrybuty ucznia i zostały podane przez wektor x i . Indeksuj klasy za pomocą j i niech atrybuty klasy j będą podane przez wektor z j . Zbiór studentów uczestniczących w klasie j jest j .$i$$i$$\mathbf{x}_i$$j$$j$$\mathbf{z}_j$$j$$A_j$

Przyjmuje się, że „siły” każdy student jest funkcją ich atrybutów oraz niektóre wartość losową, która równie dobrze może mieć zerową średnią:$s_i$

Egzamin modelujemy w klasie , dodając niezależne losowe wartości do siły każdego ucznia zapisanego do klasy i przekształcając je w szeregi. Skąd, jeśli uczeń I nie zapisze się do klasy j , ich względna pozycja R i , j jest określona przez ich pozycję w posortowanej tablicy wartości$j$$i$$j$$r_{i,j}$

Ta pozycja jest podzielona przez jeden więcej niż całkowita rejestracja klas, aby uzyskać zmienną zależną, stopień procentowy:$r_{i,j}$

Twierdzę, że wyniki regresji zależą (całkiem sporo) od wielkości i struktury losowych (nieobserwowanych) wartości oraz δ i , j . $\varepsilon_i$$\delta_{i,j}$ Wyniki zależą również od tego, w jaki sposób uczniowie są zapisani na zajęcia. Powinno to być intuicyjnie oczywiste, ale tym, co nie jest tak oczywiste - i wydaje się trudne do analizy teoretycznej - jest to, w jaki sposób i jak bardzo nieobserwowane wartości i struktury klas wpływają na regresję.

Symulacja

Bez większego wysiłku możemy symulować tę sytuację, aby tworzyć i analizować niektóre przykładowe dane. Jedną z zalet symulacji jest to, że może ona obejmować prawdziwe mocne strony uczniów, których w rzeczywistości nie da się zaobserwować. Innym jest to, że możemy zmieniać typowe rozmiary nieobserwowanych wartości, a także przypisania klas. Zapewnia to „piaskownicę” do oceny proponowanych metod analitycznych, takich jak regresja.

Na początek ustawmy generator liczb losowych dla powtarzalnych wyników i określmy rozmiar problemu. Używam, Rponieważ jest dostępny dla każdego.

set.seed(17)
n.pop <- 200      # Number of students
n.classes <- 38   # Number of classes
courseload <- 4.5 # Expected number of classes per student

Aby zapewnić realizm, twórz n.classesklasy o różnych trudnościach w dwóch skalach (matematycznej i werbalnej, z ujemną korelacją), prowadzonych na różnych poziomach akademickich (od 1 = wprowadzenie do 7 = badania) i ze zmienną łatwością. (W „łatwej” klasie różnice między poziomem uczenia się uczniów mogą być duże i / lub egzamin może w niewielkim stopniu dyskryminować uczniów. Jest to modelowane losowymi terminami które dla klasy j są zwykle duże Wyniki egzaminu będą wówczas prawie nieprzewidywalne na podstawie danych dotyczących siły ucznia. Gdy klasa nie jest „łatwa”, te losowe terminy są pomijalnie małe, a siła ucznia może idealnie określić ranking egzaminu.)$\delta_{i,j}$$j$

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

Studenci są rozłożeni na cztery lata i obdarzeni losowymi wartościami swoich atrybutów. Nie ma korelacji między żadnym z tych atrybutów:

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

$\varepsilon_i$betabeta

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

students$abilitybeta$abilitybeta$sigmaease$\delta_{i,j}$$.01$$.2$

classesspreadassignments <-...$0$

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

(Jako przykład tego, co osiągnął ten krok, patrz rysunek poniżej).

$n$$1/(n+1)$$n/(n+1)$$1/(n+1)$$0$$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))

Do tych surowych danych dołączamy atrybuty ucznia i klasy, aby utworzyć zestaw danych odpowiedni do analizy:

data <- merge(data.raw, classes, by.x="Class", by.y="row.names")
data <- merge(data, students, by.x="Student", by.y="row.names")

Zorientujmy się, sprawdzając losową próbkę danych:

> 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

Na przykład rekord 118 mówi, że uczeń nr 28 zapisał się do klasy nr 1 i uzyskał 22 miejsce (od dołu) na egzaminie, uzyskując stopień procentowy 0,957. Ogólny poziom trudności tej klasy wynosił 0,0523 (bardzo łatwo). Łącznie zapisano 22 studentów. Ten student jest studentem drugiego roku (rok 2) z 590 matematyki, 380 ustnych wyników SAT. Ich ogólna wrodzona siła akademicka wynosi 16,9. W tym czasie zapisano ich na cztery klasy.

Ten zestaw danych jest zgodny z opisem w pytaniu. Na przykład rangi procentowe rzeczywiście są prawie jednolite (tak jak muszą być w przypadku każdego kompletnego zestawu danych, ponieważ rangi procentowe dla pojedynczej klasy mają dyskretny rozkład równomierny).

Pamiętaj, ze względu na współczynniki w betatym modelu przyjęto silny związek między wynikami badań a zmiennymi przedstawionymi w tym zbiorze danych. Ale co pokazuje regresja? Zresetujmy logistykę rangi procentowej względem wszystkich obserwowalnych cech uczniów, które mogą być związane z ich umiejętnościami, a także wskaźników trudności klasowych:

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 

Wykresy diagnostyczne ( plot(fit)) wyglądają szybko: reszty są homoscedastyczne i pięknie normalne (aczkolwiek nieco krótkie, co nie stanowi problemu); brak wartości odstających; i brak niepożądanego wpływu na jakąkolwiek obserwację.

$3.5$$(-9, 1.6, 3.4, 3.9, 0.009, 0.007, -0.7)$$(*, 1, 3, 4, 0.010, 0.010, *)$$*$ oznacza współczynnik, który nie został wyraźnie określony).

level$0$$7$level

(Nawiasem mówiąc, użycie rang procentowych nieprzetworzonych w regresji nie zmienia jakościowo wyników przedstawionych poniżej.)

spread$1$$38$

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

spread$38$spread1

Tym razem R-kwadrat został znacznie poprawiony (choć nadal nie jest świetny). Jednak wszystkie współczynniki wzrosły o 20 - 100%. Ta tabela porównuje je wraz z kilkoma dodatkowymi symulacjami:

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%

spread$38$ability$2$$1$abilitysigma$0$$\varepsilon_i$$\delta_{i,j}$level

Ta szybka analiza pokazuje, że regresja, przynajmniej tak jak tutaj przeprowadzona, będzie mylić nieuniknione formy zmienności ze współczynnikami. Ponadto współczynniki zależą również (do pewnego stopnia) od tego, w jaki sposób uczniowie są podzieleni między klasy. Można to częściowo uwzględnić poprzez uwzględnienie atrybutów klasy wśród zmiennych niezależnych w regresji, jak to tutaj zrobiono, ale mimo to efekt rozkładu uczniów nie znika.

Jakikolwiek brak przewidywalności prawdziwych wyników uczniów oraz wszelkie różnice w uczeniu się uczniów i rzeczywistych wynikach egzaminów, najwyraźniej powodują, że szacunki współczynników kurczą się do zera. Wydaje się, że robią to jednolicie, co sugeruje, że względne współczynniki mogą być nadal znaczące.

$y_{ij}$

$\mu_{ij}$

$\ln(y_{ij}/(1-y_{ij})) = \mu_{ij} + e_{ij} + v_i$

nieobserwowane umiejętności ucznia są modelowane za pomocą losowego komponentu podczas gdy modeluje inne niesystematyczne nieobserwowalne. Korelację między odpowiedziami (badaniami) można zwiększyć, przyjmując ogólną strukturę kowariancji dla . Dlaczego nie ma białej (lub warstwowej / solidnej) struktury wariancji? Co więcej, niektóre korelacje odpowiedzi można uwzględnić w (zależność warunkowa).e i j e i j μ i $v_i$$e_{ij}$$e_{ij}$$\mu_{ij}$

(To tylko pomysł z mojego stronniczego doświadczenia, komentarze i krytycy są mile widziane.)

Nieobserwowalne zdolności prawdopodobnie będą skorelowane z atrybutami obserwowalnych uczniów lub egzaminów w . Dzięki tym założeniom model ten jest RE z skorelowanymi składnikami błędu, które można oszacować za pomocą ML lub estymatora dwustopniowego: pierwszy etap: transformacja wewnątrz (lub analogowa), która eliminuje . Drugi etap: OLS na przekształconym modelu. v $\mu_{ij}$$v_i$

Możesz spróbować regresji logistycznej. Logit transform rozłoży twoją zmienną odpowiedzi na rzeczywistą linię, więc nie dostaniesz absurdalnie przewidywanych procentów rangi, takich jak -3% lub + 110%.$\ln(\frac{p}{1-p})$

Idealny model w tym przypadku odwzoruje dane wejściowe (bez względu na towarzyszące im zmienne) na dane wyjściowe (ranga ucznia w klasie). Innym sposobem myślenia o tym jest mapowanie najpierw na wyniki, a następnie mapowanie tych wyników na rangę. Na razie zamierzam zignorować błąd.

$y = \sum \beta x$

$r = R(y)$

$R$$R$$y$$R(y)$

Wydaje się to być dość podobne do formy funkcjonalnej uogólnionego modelu liniowego. Myślę, że właśnie dlatego regresję logistyczną zaproponował @Mike Anderson. Jeśli wyniki twojego egzaminu byłyby rozłożone logistycznie, wówczas funkcją linku byłaby logit (odwrotnie jest to funkcja gęstości skumulowanej, na której nam zależy). Podobnie, jeśli wyniki byłyby normalnie rozłożone, funkcją probit byłaby funkcja łączenia.

Dla twojej regresji jedynym sposobem oszacowania rang jest powiedzenie „biorąc pod uwagę, że moje dane są dystrybuowane jako X, ten punkt jest w 34. percentylu”. W przeciwnym razie, skąd wiesz, co oznacza wzrost wyniku testu o dwa punkty w kategoriach rangi? Zastrzeżenie polega na tym, że musisz oszacować ten rozkład, aby wybrać funkcję linku (niektóre formy funkcjonalne znacznie ułatwią ci życie). Co więcej, ten model nie powie „byłeś szóstym najlepszym z klasy 38”, a raczej „gdyby wyniki testu były podzielone tak, jak nam się wydaje, wynik byłby na poziomie 15 percentyla”.