Berurusan dengan regresi dari variabel respons yang dibatasi secara tidak biasa

Saya mencoba untuk memodelkan variabel respons yang secara teoritis dibatasi antara -225 dan +225. Variabelnya adalah skor total yang diperoleh subjek saat bermain game. Meskipun secara teoritis dimungkinkan untuk mata pelajaran untuk mencetak +225. Meskipun demikian karena skor tidak hanya bergantung pada aksi subjek, tetapi juga aksi aksi lainnya, skor maksimum siapa pun yang mencapai 125 (ini adalah 2 pemain tertinggi yang bermain satu sama lain dapat mencetak keduanya) ini terjadi dengan frekuensi yang sangat tinggi. Skor terendah adalah +35.

Batas 125 ini menyebabkan kesulitan dengan regresi linier. Satu-satunya hal yang bisa saya pikirkan adalah mengubah skala respons menjadi antara 0 dan 1 dan menggunakan regresi beta. Jika saya melakukan ini meskipun saya tidak yakin saya benar-benar dapat membenarkan mengatakan 125 adalah batas atas (atau 1 setelah transformasi) karena mungkin untuk mencetak +225. Lebih jauh, jika saya melakukan ini, apa yang menjadi batas bawah saya, 35?

Terima kasih,

Jonathan

Meskipun saya tidak sepenuhnya yakin apa masalah Anda dengan regresi linier adalah saya sekarang sedang menyelesaikan sebuah artikel tentang bagaimana menganalisis hasil yang terbatas. Karena saya tidak terbiasa dengan regresi Beta, mungkin orang lain akan menjawab pilihan itu.

Dengan pertanyaan Anda, saya mengerti bahwa Anda mendapatkan prediksi di luar batas. Dalam hal ini saya akan pergi untuk regresi kuantil logistik . Regresi kuantitatif adalah alternatif yang sangat rapi untuk regresi linier reguler. Anda dapat melihat berbagai kuantil dan mendapatkan gambaran yang lebih baik dari data Anda daripada yang mungkin dengan regresi linier biasa. Juga tidak memiliki asumsi mengenai distribusi ¹ .

Transformasi variabel sering dapat menyebabkan efek lucu pada regresi linier, misalnya Anda memiliki signifikansi dalam transformasi logistik tetapi itu tidak diterjemahkan ke dalam nilai reguler. Ini tidak terjadi dengan kuantil, median selalu median terlepas dari fungsi transformasi. Ini memungkinkan Anda untuk mengubah bolak-balik tanpa mendistorsi apa pun. Bottai menyarankan pendekatan ini untuk membatasi hasil ² , ini merupakan metode yang sangat baik jika Anda ingin melakukan prediksi individu tetapi memiliki beberapa masalah ketika Anda tidak ingin melihat beta dan menafsirkannya dengan cara non-logistik. Rumusnya sederhana:

$logit(y) = log(\frac{y + \epsilon}{max(y) - y + \epsilon})$

Di mana adalah skor Anda dan adalah angka kecil yang berubah - ubah .$y$$\epsilon$

Berikut adalah contoh yang saya lakukan beberapa waktu lalu ketika saya ingin bereksperimen dengan itu di R:

library(rms)
library(lattice)
library(cairoDevice)
library(ggplot2)

# Simulate some data
set.seed(10)
intercept <- 0
beta1 <- 0.5
beta2 <- 1
n = 1000
xtest <- rnorm(n,1,1)
gender <- factor(rbinom(n, 1, .4), labels=c("Male", "Female"))
random_noise  <- runif(n, -1,1)

# Add a ceiling and a floor to simulate a bound score
fake_ceiling <- 4
fake_floor <- -1

# Simulate the predictor
linpred <- intercept + beta1*xtest^3 + beta2*(gender == "Female") + random_noise

# Remove some extremes
extreme_roof <- fake_ceiling + abs(diff(range(linpred)))/2
extreme_floor <- fake_floor - abs(diff(range(linpred)))/2
linpred[ linpred > extreme_roof|
    linpred < extreme_floor ] <- NA

#limit the interval and give a ceiling and a floor effect similar to scores
linpred[linpred > fake_ceiling] <- fake_ceiling
linpred[linpred < fake_floor] <- fake_floor

# Just to give the graphs the same look
my_ylim <- c(fake_floor - abs(fake_floor)*.25, 
             fake_ceiling + abs(fake_ceiling)*.25)
my_xlim <- c(-1.5, 3.5)

# Plot
df <- data.frame(Outcome = linpred, xtest, gender)
ggplot(df, aes(xtest, Outcome, colour = gender)) + geom_point()

Ini memberikan pencar data berikut, seperti yang Anda lihat jelas dibatasi dan tidak nyaman :

###################################
# Calculate & plot the true lines #
###################################
x <- seq(min(xtest), max(xtest), by=.1)
y <- beta1*x^3+intercept
y_female <- y + beta2
y[y > fake_ceiling] <- fake_ceiling
y[y < fake_floor] <- fake_floor
y_female[y_female > fake_ceiling] <- fake_ceiling
y_female[y_female < fake_floor] <- fake_floor

tr_df <- data.frame(x=x, y=y, y_female=y_female)
true_line_plot <- xyplot(y  + y_female ~ x, 
                         data=tr_df,
                         type="l", 
                         xlim=my_xlim, 
                         ylim=my_ylim, 
                         ylab="Outcome", 
                         auto.key = list(
                           text = c("Male"," Female"),
                           columns=2))

##########################
# Test regression models #
##########################

# Regular linear regression
fit_lm <- Glm(linpred~rcs(xtest, 5)+gender, x=T, y=T)
boot_fit_lm <- bootcov(fit_lm, B=500)
p <- Predict(boot_fit_lm, xtest=seq(-2.5, 3.5, by=.001), gender=c("Male", "Female"))
lm_plot <- plot(p, 
             se=T, 
             col.fill=c("#9999FF", "#BBBBFF"), 
             xlim=my_xlim, ylim=my_ylim)

Ini menghasilkan gambar berikut di mana perempuan jelas di atas batas atas:

# Quantile regression - regular
fit_rq <- Rq(formula(fit_lm), x=T, y=T)
boot_rq <- bootcov(fit_rq, B=500)
# A little disturbing warning:
# In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique

p <- Predict(boot_rq, xtest=seq(-2.5, 3.5, by=.001), gender=c("Male", "Female"))
rq_plot <- plot(p, 
             se=T, 
             col.fill=c("#9999FF", "#BBBBFF"), 
             xlim=my_xlim, ylim=my_ylim)

Ini memberikan plot berikut dengan masalah yang sama:

# The logit transformations
logit_fn <- function(y, y_min, y_max, epsilon)
    log((y-(y_min-epsilon))/(y_max+epsilon-y))


antilogit_fn <- function(antiy, y_min, y_max, epsilon)
    (exp(antiy)*(y_max+epsilon)+y_min-epsilon)/
        (1+exp(antiy))

epsilon <- .0001
y_min <- min(linpred, na.rm=T)
y_max <- max(linpred, na.rm=T)

logit_linpred <- logit_fn(linpred, 
                            y_min=y_min,
                            y_max=y_max,
                            epsilon=epsilon)

fit_rq_logit <- update(fit_rq, logit_linpred ~ .)
boot_rq_logit <- bootcov(fit_rq_logit, B=500)

p <- Predict(boot_rq_logit, 
             xtest=seq(-2.5, 3.5, by=.001), 
             gender=c("Male", "Female"))

# Change back to org. scale
# otherwise the plot will be
# on the logit scale
transformed_p <- p
transformed_p$yhat <- antilogit_fn(p$yhat,
                                    y_min=y_min,
                                    y_max=y_max,
                                    epsilon=epsilon)
transformed_p$lower <- antilogit_fn(p$lower, 
                                     y_min=y_min,
                                     y_max=y_max,
                                     epsilon=epsilon)
transformed_p$upper <- antilogit_fn(p$upper, 
                                     y_min=y_min,
                                     y_max=y_max,
                                     epsilon=epsilon)

logit_rq_plot <- plot(transformed_p, 
             se=T, 
             col.fill=c("#9999FF", "#BBBBFF"), 
             xlim=my_xlim)

Regresi kuantil logistik yang memiliki prediksi terikat sangat bagus:

Di sini Anda dapat melihat masalah dengan Beta yang dalam mode mentransformasi ulang berbeda di berbagai wilayah (seperti yang diharapkan):

# Some issues trying to display the gender factor
contrast(boot_rq_logit, list(gender=levels(gender), 
                             xtest=c(-1:1)), 
         FUN=function(x)antilogit_fn(x, epsilon))

   gender xtest Contrast   S.E.       Lower      Upper       Z      Pr(>|z|)
   Male   -1    -2.5001505 0.33677523 -3.1602179 -1.84008320  -7.42 0.0000  
   Female -1    -1.3020162 0.29623080 -1.8826179 -0.72141450  -4.40 0.0000  
   Male    0    -1.3384751 0.09748767 -1.5295474 -1.14740279 -13.73 0.0000  
*  Female  0    -0.1403408 0.09887240 -0.3341271  0.05344555  -1.42 0.1558  
   Male    1    -1.3308691 0.10810012 -1.5427414 -1.11899674 -12.31 0.0000  
*  Female  1    -0.1327348 0.07605115 -0.2817923  0.01632277  -1.75 0.0809  

Redundant contrasts are denoted by *

Confidence intervals are 0.95 individual intervals

Referensi

1. R. Koenker dan G. Bassett Jr, "Kuantitas Regresi," Econometrica: jurnal Masyarakat Ekonometrik, hlm. 33-50, 1978.
2. M. Bottai, B. Cai, dan RE McKeown, "Regresi kuantil logistik untuk hasil yang terbatas," Statistics in Medicine, vol. 29, tidak. 2, hlm. 309–317, 2010.

Bagi yang penasaran plot dibuat menggunakan kode ini:

# Just for making pretty graphs with the comparison plot
compareplot <- function(regr_plot, regr_title, true_plot){
  print(regr_plot, position=c(0,0.5,1,1), more=T)
  trellis.focus("toplevel")
  panel.text(0.3, .8, regr_title, cex = 1.2, font = 2)
  trellis.unfocus()
  print(true_plot, position=c(0,0,1,.5), more=F)
  trellis.focus("toplevel")
  panel.text(0.3, .65, "True line", cex = 1.2, font = 2)
  trellis.unfocus()
}

Cairo_png("Comp_plot_lm.png", width=10, height=14, pointsize=12)
compareplot(lm_plot, "Linear regression", true_line_plot)
dev.off()

Cairo_png("Comp_plot_rq.png", width=10, height=14, pointsize=12)
compareplot(rq_plot, "Quantile regression", true_line_plot)
dev.off()

Cairo_png("Comp_plot_logit_rq.png", width=10, height=14, pointsize=12)
compareplot(logit_rq_plot, "Logit - Quantile regression", true_line_plot)
dev.off()

Cairo_png("Scat. plot.png")
qplot(y=linpred, x=xtest, col=gender, ylab="Outcome")
dev.off()