Alternative Modelling of Binary Response Data

Chapter 7

Her ne kadar link fonksiyon olarak logit fonksiyon sık kullanılsa da, başka link fonksiyonlar da kullanılabilir.

Chapter 7.1 Probit ve Log-Log Modeller

Binary responselu logistic modellere alternatiflerden iki tanesinden bahsedelim. Bu \(\pi(x) = P(Y =1)\) şeklindeki modellerin formu \(g[\pi(x)]=\alpha+\beta_1x_1+\beta_2x_2+\cdots+\beta_px_p\) şeklindedir. g ile gösterien fonksiyon da logit dışındaki diğer link fonksiyonlardır.

7.1.1 Probit Modeller ve Latent Değişkenler

Probit modeller de logistic regresyon gibi S şeklinde ilerlerler.

\(probit[P(Y = 1)]=\alpha+\beta_1x_1+\beta_2x_2+\cdots+\beta_px_p\) Buradaki link fonksiyonu probit fonksiyonudur ve \(P(Y=1)\) ifadesini standard normal z-skoruna çevirir. Eğrinin sol tarafındaki olasılık \(P(Y=1)\)i gösterir. Örneğin, \(probit(0.05)= -1.645\)e eşittir. Bunun nedeni normal dağılımın %5’i -1.645 ve altını gösterir.

qnorm(0.05)
[1] -1.644854
qnorm(0.95)
[1] 1.644854

İstatistiksel terminolojide gözlemlenen sonucun altında yattığı düşünülen ve gözlemlenemeyen parametrelere latent variable adı verilir.

Örnek olarak \(y^*\) latent variable için gözlemlenen binary response:

\[ y = \left\{ \begin{array}{lr} 0 & y^* \leq \tau \\ 1 & y^* > \tau \end{array} \right. \]

Ek olarak bu latent değişkenin lineer model ile modellenebileceğini varsayalım.

\[ y^* = \alpha + \beta_1x1 +\cdots + \beta_px_p + \epsilon \] Buradaki \(\epsilon\) normal dağılıyorsa ve eşit varyansa sahipse probit model binary sonuç için mümkündür. Ek olarak \(var(\epsilon)=1\) olacak şekilde latent yanıt ölçeklendirilirse probit modeldeki etkiler latent variable model ile aynıdır. Böylece probit modeldeki \(\hat\beta_j\), \(x_j\)deki 1 birimlik değişiklik için \(E(y^*)\)deki değişikliği göstermiş olur(diğer parametreler sabit kaldığında). Başka bir deyişle, \(\hat\beta_j\); \(y^*\)daki standard deviasyon kaymasını gösterir.

Buradaki link fonksiyonu \(\phi^{-1}\)dir. Standart kümülatif dağılım fonksiyonunun tersidir. \(\phi^{-1}[\pi(x)]=\alpha + \beta x\).

7.1.2 Probit Modelin Değerlendirilmesi

Tek bağımsız değişkenli probit modeller için \(\pi(x)\) olasılık eğrisi normal bir dağılım fonksiyonunu andırır ve \(\mu = -\alpha/\beta\) ile \(\sigma = 1/|\beta|\)dır. Normal dağılımın %68’i ortalamanın 1 SD etrafında kaldığına göre, \(1/|\beta|\) değeri \(\pi(x) = 0.16 - 0.84\)intervaline karşılık gelir. \(\pi(x)\)’in kısmi türevi değişiklik hızını gösterir.

\[ \beta\phi(\alpha + \beta_x) \qquad \phi(\centerdot)= \text st.norm.dens.func \]

En yüksek değişiklik hızı parantezin içi sıfıra eşitken gerçekleşir. (\(x = (-\alpha/\beta)\))

Logistic regresyonda ise \(\pi(x)\) eğrisi logistic cdf eğrisidir ve standard deviasyonu \(\pi/|\beta|\sqrt{3}\)tür. \(x = (-\alpha/\beta)\) iken değişiklik hızı \(0.25\beta\) kadardır.

Her iki modelin de uyumunun yüksek olduğu durumlarda, logistic regresyondaki tahminler probit modellerin 1.6 ile 1.8 katı kadardır.

Probit modeller Fisher skorlaması ile ve Maximum Likelihood uygulanarak çözülür.

Bu örnekte, çeşitli konsantrasyonlarda karbon disülfit gazına 5 saat maruz bırakılan böceklerin ölümünü inceleyeceğiz.

prob_df <-data.frame(log.dose=c(1.691,1.724,1.755,1.784,1.811,1.837,1.861,1.884),
n=c(59,60,62,56,63,59,62,60),
y=c(6,13,18,28,52,53,61,60))
prob_df
  log.dose  n  y
1    1.691 59  6
2    1.724 60 13
3    1.755 62 18
4    1.784 56 28
5    1.811 63 52
6    1.837 59 53
7    1.861 62 61
8    1.884 60 60
options(contrasts = c("contr.treatment","contr.treatment"))
prob_1 <- glm(y/n ~ log.dose, weights = n, family = binomial(link = "probit"), data = prob_df)
summary(prob_1)

Call:
glm(formula = y/n ~ log.dose, family = binomial(link = "probit"), 
    data = prob_df, weights = n)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -34.956      2.649  -13.20   <2e-16 ***
log.dose      19.741      1.488   13.27   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 284.202  on 7  degrees of freedom
Residual deviance:   9.987  on 6  degrees of freedom
AIC: 40.185

Number of Fisher Scoring iterations: 4
library(ggeffects)
library(ggplot2)
plot(ggpredict(prob_1))

Yukarıdaki formüllerden ölüm hızının en hızlı arttığı nokta olarak \(34.956/19.741=1.77\). SD ise \(1/19.73=0.051\) şeklinde bulunur.

cbind(fitted(prob_1),prob_df$log.dose)
       [,1]  [,2]
1 0.0577367 1.691
2 0.1781060 1.724
3 0.3780390 1.755
4 0.6032833 1.784
5 0.7866532 1.811
6 0.9045852 1.837
7 0.9626183 1.861
8 0.9873227 1.884
probseq <- seq(1.67,1.87,by = 0.01)
probcdf <- data.frame( y = pnorm(probseq,1.77,0.05), x = 1:21)
ggplot(probcdf, aes(x = x, y = y)) + geom_point() +
  geom_hline(yintercept = c(0.16,0.5,0.84), color = 1:3, linetype = 3, linewidth = 0.8) +
  geom_vline(xintercept = c(6,11,16),color = 1:3, linetype = 2, linewidth = 0.8) + theme_bw()

Görüldüğü gibi toksik gaz dozu arttıkça ölüm ihtimali 0.0577’den 0.987’ye yükselmektedir. X’teki her 0.1 birimlik artışa karşılık(ör: 1.7 den 1.8’e) latent değişkenin \(y^*\) nin yaklaşık \(0.1*19.73 = 2SD\) kadar yükseleceğini söyleyebiliriz.

7.2 Bayesçi Yaklaşım

Modellerin birden çok parametresi olması doğaldır. Araştırmacı bu parametrelerin bir kısmı hakkında önceden bilgiye sahip olabilir. Basitçe, \(\beta\)değerlerinin prior dağılımlarının olası tüm parametre değerleri boyunca sabit olduğu düşünülebilir. O zaman posterior dağılım kısaca likelihood fonksiyonunun bir sabitle çarpımıdır. Posterior dağılımın mod değeri de ML tahminine eşittir. Ancak verisetine göre posterior dağılımlar oldukça çarpık da olabilir. Bu gibi durumlarda posterior mean ve mode arasında ciddi farklılıklar görünür.

\(\beta\)nın prior dağılımında çok değişkenli normal dağılım oldukça esneklik sağlar. Eğer tüm parametreler için \(N (0,\sigma^2)\) prior olarak kullanılacaksa, bağımsız değişkenler standartlaştırılmalıdır. Bu şekilde etkilerin karşılaştırılabilmesi mümkündür. Diğer türlü ölçeklendirme göz önünde bulundurulmalıdır.

Normal prior değerlerinde yüksek SD değerleri, prior olasılık dağılımının U şeklinde olmasına yol açar. Olasılıkların yarısı 0’a çok yakınken diğer yarısı 1’e çok yakındır. Her ne kadar bu durum işe yarar gibi gözükse de posterior dağılım üzerine etkisi oldukça sınırlıdır ve ML tahmincisine benzer sonuçlar verir. Bunun yerine uniform SD kullanılması daha faydalı olacaktır.

Binomial regression modellerinde Jerreys prior dağılımı genellikle uygundur. Probit ve logit linklerde, bu prior simetriktir, ve 0 noktasında unimodaldir. Her iki ucu da çok değişkenli t dağılımına kıyasla daha incedir.

Şimdi bir örnek görelim. Bu örnekte endometrium kanserini etkileyen faktörler araştırılmış. HG = 1:high-grade, 0: low-grade; risk faktörleri olarak da NV(x1): neovaskülarizasyon(0,1), PI(x2) : pulsatility index(range 0 - 49), EH(x3): endometrium doku kalınlığı (range 0.27-3.61)

library(vcd)
Loading required package: grid
library(vcdExtra)
Loading required package: gnm
library(ggplot2)
library(bayesplot)
This is bayesplot version 1.10.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
   * Does _not_ affect other ggplot2 plots
   * See ?bayesplot_theme_set for details on theme setting
library(rstanarm)
Loading required package: Rcpp
This is rstanarm version 2.26.1
- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
  options(mc.cores = parallel::detectCores())

Attaching package: 'rstanarm'
The following object is masked from 'package:gnm':

    se
theme_set(bayesplot::theme_default())

Endo <- read.table("http://www.stat.ufl.edu/~aa/cat/data/Endometrial.dat",header=TRUE)
head(Endo)
  NV PI   EH HG
1  0 13 1.64  0
2  0 16 2.26  0
3  0  8 3.14  0
4  0 34 2.68  0
5  0 20 1.28  0
6  0  5 2.31  0
str(Endo)
'data.frame':   79 obs. of  4 variables:
 $ NV: int  0 0 0 0 0 0 0 0 0 0 ...
 $ PI: int  13 16 8 34 20 5 17 10 26 17 ...
 $ EH: num  1.64 2.26 3.14 2.68 1.28 2.31 1.8 1.68 1.56 2.31 ...
 $ HG: int  0 0 0 0 0 0 0 0 0 0 ...

Main-effects model olarak:

\[ logit[P(Y =1)] = \alpha + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 \] Burada \(x_2\) ve \(x_3\) standardize edeceğiz.

Endo$PI2 <- scale(Endo$PI); Endo$EH2 <- scale(Endo$EH)
endo_ml <- glm(HG ~ NV + PI2 + EH2, family=binomial, data=Endo)
summary(endo_ml)

Call:
glm(formula = HG ~ NV + PI2 + EH2, family = binomial, data = Endo)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -1.2517     0.3688  -3.394 0.000688 ***
NV            18.1856  1715.7509   0.011 0.991543    
PI2           -0.4217     0.4432  -0.952 0.341333    
EH2           -1.9219     0.5599  -3.433 0.000597 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 104.903  on 78  degrees of freedom
Residual deviance:  55.393  on 75  degrees of freedom
AIC: 63.393

Number of Fisher Scoring iterations: 17

Görüldüğü gibi tama yakın ayrışma nedeniyle \(\beta_1\)’in ML tahmini \(\infty\) olarak bulundu. Neovaskülarizasyonun önemli bir kanser belirteci olduğu göz önünde bulundurulduğunda bu değişkeni modelden çıkartmak çok mantıklı bir yaklaşım değildir.

Bayesçi yaklaşımda \({\beta_j}\) için \(N (\mu,\sigma^2)\) prior dağılımlar alınır. Prior dağılımlar \(\mu = 0\) ve (0,1) şeklinde kodlandığında prior etkilerin işaretinin + ya da - olduğunun tam bilinemediği varsayılmış olur. Bunun yerine, \(x1\) değişkenini (-0.5,0.5) şeklinde kodlayalım. Böylece prior dağılım simetrik olur ve her grubun logitleri aynı prior dağılıma ve ortalamaya sahip olur. Bununla birlikte \(\beta_1\)in şartlı log odds ratio olarak yorumlanması değişmez

Endo$NV2 <- Endo$NV - 0.5
bstan1 <- stan_glm(HG ~ NV2 + PI2 + EH, data = Endo,
                   family = binomial(link = "logit"),
                   prior = normal(0,1),
                   prior_intercept = normal(0,1),
                   cores =8, 
                   iter = 1000,
                   seed = 123)

bstan2 <- stan_glm(HG ~ NV2 + PI2 + EH, data = Endo,
                   family = binomial(link = "logit"),
                   prior = normal(0,10),
                   prior_intercept = normal(0,1),
                   cores =8, 
                   iter = 1000,
                   seed = 123)
coef(bstan1) |> round(3)
(Intercept)         NV2         PI2          EH 
      3.504       1.712      -0.173      -2.160 
coef(bstan2) |> round(3)
(Intercept)         NV2         PI2          EH 
      8.051       7.984      -0.432      -3.129 
round(posterior_interval(bstan1, prob = 0.95), 2) 
             2.5% 97.5%
(Intercept)  1.91  5.40
NV2          0.41  3.17
PI2         -0.77  0.41
EH          -3.29 -1.19
round(posterior_interval(bstan2, prob = 0.95), 2)
             2.5% 97.5%
(Intercept)  4.02 13.22
NV2          2.23 17.31
PI2         -1.38  0.37
EH          -5.04 -1.65
coef(endo_ml) |> round(3)
(Intercept)          NV         PI2         EH2 
     -1.252      18.186      -0.422      -1.922 
confint(endo_ml) |> round(3)
Waiting for profiling to be done...
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
collapsing to unique 'x' values
              2.5 %  97.5 %
(Intercept)  -2.057  -0.583
NV          -69.228 506.194
PI2          -1.371   0.382
EH2          -3.169  -0.951
#launch_shinystan(bstan1)
plot(bstan1)

plot(bstan2)

Görüldüğü üzere bayesçi yaklaşımla neovaskülarizasyonun yeni endometrium kanseri gelişiminde etkili olduğunu anlıyoruz.

7.2.5 Olasılığa Göre Prior Belirleme

trauma_nm <- read.table("http://www.ics.uci.edu/~wjohnson/BIDA/Ch8/trauma300.txt",header=TRUE)
str(trauma_nm)
'data.frame':   300 obs. of  6 variables:
 $ ID   : int  2979 1167 116 1118 3126 2173 25 917 2828 979 ...
 $ death: int  0 0 0 0 0 0 1 0 0 0 ...
 $ ISS  : int  1 9 29 10 2 9 25 2 11 16 ...
 $ TI   : int  1 0 0 0 1 0 1 0 1 1 ...
 $ RTS  : num  7.84 7.84 2.93 7.84 7.84 ...
 $ AGE  : int  25 63 32 61 43 27 30 25 42 80 ...

y: ölüm (1:ölü, 0:yaşıyor) x1: injurity severity score (0-75) x2: trauma score (çeşitli parametrelerin birleşiminden oluşturulan bir skorlama: 0-75) x3: yaş x4: type of injury (1: penetran, 0: künt)

model: \(logit[P(Y =1)] = \alpha + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_4x_4 + \beta_5(x_3x_4)\)

bstan3 <- stan_glm(death ~ ISS + RTS + AGE*TI , data = trauma_nm,
                   family = binomial(link = "logit"),
                   prior = normal(0,1),
                   prior_intercept = normal(0,1),
                   cores =8, 
                   iter = 10000,
                   seed = 123)

cbind(round(coef(bstan3),2),round(posterior_interval(bstan3, prob = 0.95), 2))
                   2.5% 97.5%
(Intercept) -1.91 -4.78  0.96
ISS          0.07  0.02  0.12
RTS         -0.52 -0.84 -0.22
AGE          0.04  0.02  0.07
TI           0.43 -1.14  1.96
AGE:TI       0.01 -0.04  0.05

Şimdi de Bedrick et al. çalışmasından spesifik priorları alalım. çalışma linki

mean_priors <- normal(location = c(25,7.84,60,0,0),
                      scale = c(10,10,10,10,10))

bstan4 <- stan_glm(death ~ ISS + RTS + AGE*TI , data = trauma_nm,
                   family = binomial(link = "logit"),
                   prior = mean_priors,
                   cores =8, 
                   iter = 10000,
                   seed = 123)

bstan4
stan_glm
 family:       binomial [logit]
 formula:      death ~ ISS + RTS + AGE * TI
 observations: 300
 predictors:   6
------
            Median MAD_SD
(Intercept) -2.6    1.6  
ISS          0.1    0.0  
RTS         -0.6    0.2  
AGE          0.1    0.0  
TI           1.4    1.4  
AGE:TI       0.0    0.0  

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
#launch_shinystan(bstan4)

7.3 Şartlı Logistic Regression

ML ile bulunan lojistik regresyon tahmincileri, örneklem sayısı parametre sayısına oranla büyük olduğunda iyi iş çıkarırlar. Ancak bu durum geçerli değilse modeller şartlı ML (conditional maximum likelihood) ile iyileştirilebilirler.

7.3.1 Conditional Likelihood

\(y_i\) iki olasılıklı i katılımcısına ait deneyin sonucu, \(x_ij\) ise j prediktörünün bu katılımcı için aldığı değer olsun (j = 1,…,p). Her bir \(y_i\) tek bir denemeyi temsil etsin. (yani \(n_i = 1\)). Bu durumda model: \[ P(Y_i = y_i) = \frac{exp[y_i(\alpha + \Sigma^p_{j=1}\beta_jx_{ij})]}{1 + exp(\alpha + \Sigma^p_{j=1}\beta_jx_{ij})} \]

Sonraki adımda ise interceptler ve katsayılar birbirinden ayrıştılır. Özellikle ilgi duymadığımız parametreler(nuisance) yeterli istatistiklerine göre elimine edilebilirler. Interceptler için yeterli istatistik \(\Sigma_iy_i\) yani toplam başarı sayısıdır. Katsayılar için ise bu istatistik \(\Sigma_iy_ix_{ij}\)dir.

\[ logit[P(Y_{it}=1)] = \alpha_i + \beta x_{it} \]

Modeldeki tüm katılımcıların kendi intercepti bulunur. X’e addedilemeyen varyason interceptlerde tutulabilir. Bu interceptler elimine edildikten sonra da \(\beta\) hakkında çıkarsamalar yapılabilir.

7.3.2 Lojistik Regresyon için Küçük Örneklemlerde Çıkarsamalar

Bağımsız değişkenler kategorik olduğunda şartlı likelihood ile “exact” örneklem dağılımları oluşturulabilir. Bu yöntemlerle P-değerleri yaklaşık değil tam olarak bulunabilir. Ancak çok basit modeller dışında hesaplanması zaman ister ve özel yazılım gerekebilir.

7.3.3 2*2 Kontenjans Tabloları için Küçük Örneklemlerde Çıkarsamalar

Temel bir modelle örnek verelim. İki değerlikli tek prediktörü x olan aşağıdaki modeli ele alalım. \(logit[P(Y_i = 1)] = \alpha + \beta x_i\) ve \(i = 1,\cdots,N\) Bu modelde \(x_i=1\) 1. satırı ve \(x_i=0\) 2. satırı temsil etsin. \(\alpha\) için yeterli istatistik \(\Sigma_iy_i\) ilk sütun toplamına eşittir. \(\beta\) için yeterli istatistik ise \(T = \Sigma_iy_ix_{ij}\) ilk satırdaki başarı sayılarının toplamına eşittir. Benzer şekilde iki satır için yeterli istatistik de bu satırlardaki başarı sayılarıdır. Bu iki değişkeni t ve s ile gösterelim. \(n_1\) ve \(n_2\) de satır toplamları olsun.

Sütun toplamları sabit olduğundan, t değişkeninin dağılımı noncentral hipergeometrik olasılık ile gösterilebilir ve sadece \(\beta\) ile ilişkilidir. Bu ilişkiden çıkan \(\beta = 0\) exact testi, daha önce gösterilen 2*2 tablolar için kullanılan Fisher’s exact testi ile aynıdır.

7.3.4 Lineer-Logit Modeller için Küçük Örneklemlerde Çıkarsamalar

\(logit(\pi_i)= \alpha + \beta x_i\) şeklinde gösterilen bu modeller; I*2 hücreli sıralı satırlara sahip tablolarda kullanılabilir. \(\beta = 0\) olduğunda her bir hücredeki gözlem sayısı çok değişkenli hipergeometrik dağılımla açıklanabilir.\(H_0:\beta = 0\) testi Cochran-Armitage istatistiğine denk olduğundan bu test de exact testlerdendir.

7.3.5 22K tablolarda Şartlı Bağımsızlık için Küçük Örneklem Testleri

X ve Ynin şartlı bağımsızlık testi olan Cochran-Mantel-Haenszel testinde \(\beta\)nın yeterli istatistiği \(logit(\pi_{ik})= \alpha + \beta x_i + \beta^Z_k \qquad i = 1,2, \quad k = 1,\cdots,K\)

Exact testler \(T = \Sigma_k n_{11k}\) eşitliğini kullanırlar. \(H_\alpha: \beta >0\) testi için P-değeri, \(T \geq t_{obs}\) olasılığına eşittir. K = 1 olduğunda, Fisher’s exact testine denk gelir.

7.3.7 Exact Test Kullanmakla ilgili İstenmeyen Sonuçlar

Bu tarz testler, aynı Fisher’s exact testi gibi konservatif yaklaşım gösterirler. Özellikle gözlem sayısı düşük olduğunda ya da veriseti dengesiz olduğunda (yani gözlemler homojen dağılmamışsa) bu önerme daha da kuvvetlenir.

Şartlı dağılımın tanımlandığı gözlem değerleri özgünse (başka gruplarda bu değerler bulunmuyorsa) şartlandırma yapılamaz ve P-değeri 1 olarak çıkar. Bu durum özellikle elimine edilen değişkenlerden en az biri sürekliyse ve gözlemlenen değerlerin arası eşit değilse gerçekleşir.

Son olarak, şartlı yaklaşımlarda sadece canonical link kullanılabilir. Yani, logisic modeller bu şekilde incelenebilirken probit modellerde bu durum söz konusu değildir.

ex7<- data.frame(race=c("black","white"),
                 promote=c("yes","no"),
                 month=c("july","august","september"),
                 count=c(0,4,7,16,0,4,7,13,0,2,8,13))

mantelhaen.test(xtabs(count ~ race + promote + month, data = ex7), exact = TRUE)

    Exact conditional test of independence in 2 x 2 x k tables

data:  xtabs(count ~ race + promote + month, data = ex7)
S = 0, p-value < 2.2e-16
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.0000000 0.0136335
sample estimates:
common odds ratio 
                0

7.4 Smoothing(Düzleştirme): Kernel, Cezalandırılmış Likelihood, GAM

7.4.2 Kernel Smoothing

Herhangi bir veri tipinde kullanılabilir. En temel halinde modelden tamamen bağımsızdır. Parametrik dağılım varsayımlarına takılı kalmadan nonparametrik yollarla olasılık yoğunluk fonksiyonlarını tahmin etmede kullanılabilirler.

İlk olarak, çoklu kontenjans tablosunda ortak hücre olasılığı olan \(\pi\)yi hücrelerdeki gözlem oranları olan \(p\)ye göre düzleştirelim. Kernel tahmincileri \(\hat\pi=Kp\) şeklindedirler.

Sütun toplamları olan K değeri 1 olduğuna göre, \(\hat\pi\) toplamı da 1dir. Bu tarz kerneller aşağıda gösterilen şekildeki olasılık tahminleri için kullanılırlar. \[ \hat\pi_i=(1-\lambda)p_i + \lambda(smoother_i) \]

\(\lambda\) değeri düzleştirmeyi ayarlamak için kullanılır. Düşük değerler çok zigzaglı eğrilere neden olurken yüksek değerler modeldeki varyasyonu çok yumuşatıp verinin önemini yitirmesine neden olabilir. Alacağı değerler veri tipine göre değişir. Örnek olarak ordinal datalarda düzleştirici faktörü komşu hücrelere daha yüksek ağırlıklar verir.

Örnekte omurga cerrahisi sonrası kifoz gelişme riskini görüyoruz.

library(gam)
Loading required package: splines
Loading required package: foreach
Loaded gam 1.22-3
Kyph <- read.table("http://www.stat.ufl.edu/~aa/cat/data/Kyphosis.dat",header=TRUE)
head(Kyph)
   x y
1 12 1
2 15 1
3 42 1
4 52 1
5 59 1
6 73 1
k1 <- gam(y ~ s(x,df = 1), family = binomial, data = Kyph) #lineer
k2 <- gam(y ~ s(x,df = 2), family = binomial, data = Kyph) #quadratik
k3 <- gam(y ~ s(x,df = 3), family = binomial, data = Kyph) #cubic
plot(y ~ x, xlab = "Age", ylab = "Kyphosis", data = Kyph)
curve(predict(k1,data.frame(x = x), type = "resp"), add = TRUE, lty = 1, col = "gray")
curve(predict(k2,data.frame(x = x), type = "resp"), add = TRUE, lty = 2, col = "red")
curve(predict(k3,data.frame(x = x), type = "resp"), add = TRUE, lty = 3, col = "blue")

Görüldüğü gibi lineer terimde sürekli bir artış varmış gibi görünürken, kuadratik ya da kübik terimler eklendiğinde aslında çok küçük ve çok büyük yaşlarda kifoz riskinin azaldığını görüyoruz.

7.4.4 En Yakın Komşu Düzleştirme

Kernel tip düzleştiriciler, tahminlerini en yakın noktaya göre yapabilir. Bu yöntemlerinden en basit olanı en yakın komşu düzleştirmedir. Genellikle sınıflandırma için kullanılırlar. Benzer özelliklere sahip diğer kişilere göre ilgili kişiye tahminde bulunulur.

7.4.5 Cezalandırılmış Likelihood

Bu tarz modellerde log-likelihood fonksiyonu \(\lambda\) fonksiyonu ile cezalandırılır. Yöntemine göre katsayıları sıfır yapar ya da sıfıra yakınlaştırır.

kuadratik terimlerle cezalandırma yapan yöntemler L2-norm methodlardır ve ridge regresyonların analogudur.

\(\lambda(\beta)= \lambda\Sigma_j\beta^2_j\)

Buna karşın, L1-norm methodları \(\lambda(\beta)= \lambda\Sigma_j|\beta_j|\) şeklinde bir ceza kullanırlar. Lineer regresyondaki LASSO regresyonun analogudur.

Son olarak L0-norm metodunda \(\lambda(\beta)\) değerleri, sıfır olmayan \(\beta_j\) değerlerinin oranına göre belirlenir.

Kullanılacak \(\lambda\) değerleri cross-validasyon yöntemleriyle belirlenebilir.

7.4.7 Firth Cezalı Likelihood ile Lojistik Regresyon

Bu yöntemle ML tahmincilerinin yanlılığı düşer. Bayesçi yaklaşımda kullanılan Jeffreys prior ile denktir.

library(logistf)
firth1 <- logistf(HG ~ NV2 + PI2 + EH2, family = binomial, data = Endo)
summary(firth1)
logistf(formula = HG ~ NV2 + PI2 + EH2, data = Endo, family = binomial)

Model fitted by Penalized ML
Coefficients:
                  coef  se(coef) lower 0.95 upper 0.95      Chisq            p
(Intercept)  0.3080221 0.7585398 -0.9754702  2.7888475  0.1689951 6.810064e-01
NV2          2.9292730 1.4649738  0.6097244  7.8546317  6.7984572 9.123668e-03
PI2         -0.3474416 0.3788410 -1.2443165  0.4044667  0.7468285 3.874822e-01
EH2         -1.7243006 0.4990014 -2.8903284 -0.8162243 17.7593175 2.506867e-05
            method
(Intercept)      2
NV2              2
PI2              2
EH2              2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Likelihood ratio test=43.65582 on 3 df, p=1.78586e-09, n=79
Wald test = 21.66967 on 3 df, p = 7.641285e-05

Hatırlayacağımız üzere, klasik ML metodların NV2nin güven aralığı sonsuza uzanıyordu. Görüldüğü üzere bu değeri 0’a yaklaştırmış oldu. Diğer değerlerin çok değişmediği de görülüyor.

7.4.9 General Additif Modeller

GLM’ye ek olarak nonnormal dağılımların modellenmesi mümkündür. GLM’deki lineer prediktörü additif düzleştirici fonksiyonla değiştirir(s).

\[ g(\mu_i)= \Sigma_js_j(x_{ij}) \] Sıklıkla kullanılan düzleştiricilerden birisi kübik splinelardır. Kesikli intervallerde ayrı ayrı kübik polinomlar oluşturur ve “knot” adı verilen noktalar ile bu intervaller birleştirilir. Knotların ve fonksiyon derecelerinin belirlenmesi başka bir konudur.

library(gam)
data(kyphosis)
summary(kyphosis)
    Kyphosis       Age             Number           Start      
 absent :64   Min.   :  1.00   Min.   : 2.000   Min.   : 1.00  
 present:17   1st Qu.: 26.00   1st Qu.: 3.000   1st Qu.: 9.00  
              Median : 87.00   Median : 4.000   Median :13.00  
              Mean   : 83.65   Mean   : 4.049   Mean   :11.49  
              3rd Qu.:130.00   3rd Qu.: 5.000   3rd Qu.:16.00  
              Max.   :206.00   Max.   :10.000   Max.   :18.00  
k4 <- gam(Kyphosis ~ s(Age,4) + Number, family = binomial, data = kyphosis)
k4$coefficients
 (Intercept)    s(Age, 4)       Number 
-4.392677317  0.009490128  0.558604472 
summary(k4)

Call: gam(formula = Kyphosis ~ s(Age, 4) + Number, family = binomial, 
    data = kyphosis)
Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-1.7844 -0.7190 -0.3593 -0.1199  2.3198 

(Dispersion Parameter for binomial family taken to be 1)

    Null Deviance: 83.2345 on 80 degrees of freedom
Residual Deviance: 63.1073 on 75 degrees of freedom
AIC: 75.1072 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
          Df Sum Sq Mean Sq F value  Pr(>F)   
s(Age, 4)  1  0.382  0.3818  0.4459 0.50632   
Number     1  7.441  7.4412  8.6914 0.00426 **
Residuals 75 64.212  0.8562                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova for Nonparametric Effects
            Npar Df Npar Chisq P(Chi)
(Intercept)                          
s(Age, 4)         3     5.8993 0.1166
Number                               
k5 <- gam(Kyphosis ~ bs(Age, knots =  5), family = binomial, data = kyphosis)
summary(k5)

Call: gam(formula = Kyphosis ~ bs(Age, knots = 5), family = binomial, 
    data = kyphosis)
Deviance Residuals:
       Min         1Q     Median         3Q        Max 
-9.994e-01 -8.498e-01 -4.649e-01 -2.107e-08  2.200e+00 

(Dispersion Parameter for binomial family taken to be 1)

    Null Deviance: 83.2345 on 80 degrees of freedom
Residual Deviance: 71.9723 on 76 degrees of freedom
AIC: 81.9723 

Number of Local Scoring Iterations: 20 

Anova for Parametric Effects
                   Df Sum Sq Mean Sq F value Pr(>F)
bs(Age, knots = 5)  4  4.660 1.16499  1.3094 0.2741
Residuals          76 67.617 0.88969               
k5$coefficients
        (Intercept) bs(Age, knots = 5)1 bs(Age, knots = 5)2 bs(Age, knots = 5)3 
          -723.5834            720.8748            722.8853            726.0233 
bs(Age, knots = 5)4 
           717.3073

Cezalandırma ya da Düzleştirme yöntemlerinde dikkat edilmesi gereken düzleştirme derecesinin ya da kaçıncı dereceden terim ekleneceğinin genellikle ad hoc gerçekleştirilmesidir. Buna karşın Bayesçi yaklaşımın teorik temelleri daha sağlamdır.

çalışmalinki1

çalışmalinki2

çalışmalinki3