Contoh kasus: pembenihan

library(hglm.data)

## Loading required package: Matrix

## Loading required package: MASS

## Loading required package: sp

data("seeds")

Pembenihan (Crowder, 1978)

Contoh kasus ini diambil dari example 17.6 pada buku All in Likelihood dengan penulis Yudi Pawitan. Data tersebut diambil dari penelitian Crowder (1978), dimana r adalah banyak benih yang tumbuh dari n benih yang disemai. Data merupakan percobaan faktorial 2x2 dengan dua faktor seed dan extract. Faktor seed ada 2 jenis yaitu O75 dan O73. Faktor extract juga ada 2 jenis yaitu bean dan cucumber. Dalam percobaan tersebut, pembenihan dilakukan pada 21 plate. Faktor seed dan extract sebagai fixed effect dan plate sebagai random effect. Berikut ini data yang diambil dari library(hglm.data) dan disesuaikan untuk pengolahan data.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ dplyr::select() masks MASS::select()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

#Mengubah seed O75 menjadi 0 dan O73 menjadi 1, extract bean menjadi 0 dan cucumber menjadi 1. Artinya O75 dan bean menjadi kategori referensi
data <-seeds %>%
  mutate(
    seed = ifelse(seed == "O75", 0, 1),        # 0 jika O75, 1 jika O73
    extract = ifelse(extract == "Bean", 0, 1)  # 0 jika Bean, 1 jika Cucumber
  )
#mendefinisikan data baru dengan menambahkan kolom rr=banyak benih yang tidak tumbuh (gagal)
data_seeds<-data.frame(rr=(seeds$n-seeds$r),data)
data_seeds

##    rr plate seed extract  r  n
## 1  29     1    0       0 10 39
## 2  39     2    0       0 23 62
## 3  58     3    0       0 23 81
## 4  25     4    0       0 26 51
## 5  22     5    0       0 17 39
## 6   8     6    1       0  8 16
## 7  20     7    1       0 10 30
## 8  20     8    1       0  8 28
## 9  22     9    1       0 23 45
## 10  4    10    1       0  0  4
## 11  1    11    0       1  5  6
## 12 21    12    0       1 53 74
## 13 17    13    0       1 55 72
## 14 19    14    0       1 32 51
## 15 33    15    0       1 46 79
## 16  3    16    0       1 10 13
## 17  9    17    1       1  3 12
## 18 19    18    1       1 22 41
## 19 15    19    1       1 15 30
## 20 19    20    1       1 32 51
## 21  4    21    1       1  3  7

Model GLM

Berikut pengolahan data untuk mendapatkan model GLM

fit_glm<-glm(cbind(r,rr)~extract*seed,family=binomial(link="logit"),data=data_seeds)
summary(fit_glm)

## 
## Call:
## glm(formula = cbind(r, rr) ~ extract * seed, family = binomial(link = "logit"), 
##     data = data_seeds)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.5582     0.1260  -4.429 9.46e-06 ***
## extract        1.3182     0.1775   7.428 1.10e-13 ***
## seed           0.1459     0.2232   0.654   0.5132    
## extract:seed  -0.7781     0.3064  -2.539   0.0111 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 98.719  on 20  degrees of freedom
## Residual deviance: 33.278  on 17  degrees of freedom
## AIC: 117.87
## 
## Number of Fisher Scoring iterations: 4

Interpretasi dari model GLM di atas:

Dari nilai p terlihat bahwa extract dan interaksi antara seed dan extrac berpengaruh signifikan terhadap peluang pertumbuhan benih.
Nilai intercept=-0,5582 artinya ketika seed=O75 dan extract=bean maka nilai log-odds benih tumbuh sebesar -0,5582. Sehingga nilai peluang benih akan tumbuh pada kondisi tersebut adalah 36,39%
Nilai penduga extract=1,3182 artinya jika benih berubah dari bean ke cucumber maka nilai log-odds meningkat sebesar 1,3182
Nilai penduga seed=0,1459 artinya jika ekstrak berubah dari O75 ke O73 maka nilai log-odds meningkat sebesar 0,1459
Nilai penduga extract:seed=-0,7781 artinya efek interaksi antara extract=cucumber dan seed=O73 menurunkan nilai log-odds sebesar 0,7781
Contoh prediksi: jika extract=bean (0) dan seed=O73(1) maka nilai log-odds=-0,5582+1,3182*0+0,1459*1-0,7781*0*1=-0,4123 sehingga diperoleh peluang benih tumbuh sebesar 39,84%

Model GLMM : Laplace

Berikut pengolahan data untuk perhitungan GLMM dengan pendekatan Laplace

library(lme4)
glmm_laplace<-glmer(cbind(r,rr)~seed*extract+(1|plate),family=binomial(link="logit"),data=data_seeds)
summary(glmm_laplace)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(r, rr) ~ seed * extract + (1 | plate)
##    Data: data_seeds
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     117.5     122.8     -53.8     107.5        16 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.60042 -0.78762  0.04326  0.72641  1.24275 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  plate  (Intercept) 0.05503  0.2346  
## Number of obs: 21, groups:  plate, 21
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.54848    0.16608  -3.302 0.000958 ***
## seed          0.09743    0.27736   0.351 0.725390    
## extract       1.33681    0.23618   5.660 1.51e-08 ***
## seed:extract -0.81004    0.38417  -2.109 0.034986 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) seed   extrct
## seed        -0.600              
## extract     -0.702  0.412       
## seed:extrct  0.431 -0.705 -0.617

ranef(glmm_laplace)

## $plate
##     (Intercept)
## 1  -0.158499495
## 2   0.009041983
## 3  -0.182688388
## 4   0.241332870
## 5   0.099402378
## 6   0.080617869
## 7  -0.066270464
## 8  -0.117039535
## 9   0.188888936
## 10 -0.081427028
## 11  0.044993839
## 12  0.062779032
## 13  0.166013554
## 14 -0.104336158
## 15 -0.231904303
## 16  0.050760258
## 17 -0.152424836
## 18  0.025502229
## 19 -0.022114516
## 20  0.179552078
## 21 -0.031746814
## 
## with conditional variances for "plate"

Interpretasi dari model GLMM Laplace di atas:

Nilai AIC = 117,5 yang nilainya lebih kecil daripada nilai AIC GLM= 117,87 sehingga model GLMM Laplace dikatakan lebih baik.
Dari nilai p terlihat bahwa extract dan interaksi antara seed dan extrac berpengaruh signifikan terhadap peluang pertumbuhan benih.
Nilai penduga seed=0,09743 artinya jika seed berubah dari O75 ke O73 maka nilai log-odds meningkat sebesar 0,09743
Nilai penduga extract= 1,33681 artinya jika extract berubah dari bean ke cucumber maka nilai log-odds meningkat sebesar 1,33681
Nilai penduga seed:extract=-0,81004 artinya efek interaksi antara seed=O73 dan extract=cucumber menurunkan nilai log-odds sebesar 0,81004
Contoh prediksi: jika extract=bean (0) dan seed=O73(1) pada plate ke-4 maka nilai log-odds=-0,54848+0,09743*1+1,33681*0-0,81004*1*0+0,241332870=-0,20972 sehingga diperoleh peluang benih tumbuh sebesar 44,78%

Model GLMM : Adaptive Gauss-Hermite Quadrature (AGHQ)

Berikut pengolahan data untuk mendapatkan model GLMM : AGHQ

glmm_aghq<-glmer(cbind(r,rr)~seed*extract+(1|plate),family=binomial(link="logit"),data=data_seeds, nAGQ=10)
summary(glmm_aghq)

## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(r, rr) ~ seed * extract + (1 | plate)
##    Data: data_seeds
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##      40.9      46.1     -15.5      30.9        16 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.59639 -0.78239  0.04285  0.72555  1.23519 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  plate  (Intercept) 0.05582  0.2363  
## Number of obs: 21, groups:  plate, 21
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.54842    0.16660  -3.292 0.000995 ***
## seed          0.09693    0.27805   0.349 0.727378    
## extract       1.33703    0.23694   5.643 1.67e-08 ***
## seed:extract -0.81042    0.38518  -2.104 0.035379 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) seed   extrct
## seed        -0.600              
## extract     -0.702  0.412       
## seed:extrct  0.432 -0.705 -0.617

ranef(glmm_aghq)

## $plate
##     (Intercept)
## 1  -0.160058375
## 2   0.009083552
## 3  -0.184039847
## 4   0.243349613
## 5   0.100319833
## 6   0.081645686
## 7  -0.066834744
## 8  -0.118162256
## 9   0.190728591
## 10 -0.082517703
## 11  0.045578316
## 12  0.063128110
## 13  0.167204502
## 14 -0.105365166
## 15 -0.233698076
## 16  0.051354334
## 17 -0.154217671
## 18  0.025945515
## 19 -0.022167918
## 20  0.181307057
## 21 -0.032110801
## 
## with conditional variances for "plate"

Interpretasi dari model GLMM AGHQ di atas:

Nilai AIC = 40,9 yang nilainya lebih kecil daripada nilai AIC laplace= 117,5 sehingga model GLMM AGHQ dikatakan lebih baik.
Dari nilai p terlihat bahwa extract dan interaksi antara seed dan extrac berpengaruh signifikan terhadap peluang pertumbuhan benih.
Nilai penduga seed=0,09693 artinya jika seed berubah dari O75 ke O73 maka nilai log-odds meningkat sebesar 0,09693
Nilai penduga extract= 1,33703 artinya jika extract berubah dari bean ke cucumber maka nilai log-odds meningkat sebesar 1,33703
Nilai penduga seed:extract=-0,81042 artinya efek interaksi antara seed=O73 dan extract=cucumber menurunkan nilai log-odds sebesar 0,81042
Contoh prediksi: jika extract=cucumber (1) dan seed=O75(0) pada plate ke-18 maka nilai log-odds=-0,54842+0,09693*0+1,33703*1-0,81042*0*1+0,02595=0,8145 sehingga diperoleh peluang benih tumbuh sebesar 69,308%

Model GLMM :PQL

Berikut pengolahan data untuk model GLMM dengan pendekatan PQL

glmm_pql<-glmmPQL(cbind(r,rr)~seed*extract, random=~1|plate, family=binomial(link="logit"),data=data_seeds)

## iteration 1

## iteration 2

## iteration 3

summary(glmm_pql)

## Linear mixed-effects model fit by maximum likelihood
##   Data: data_seeds 
##   AIC BIC logLik
##    NA  NA     NA
## 
## Random effects:
##  Formula: ~1 | plate
##         (Intercept) Residual
## StdDev:    0.215445  1.04035
## 
## Variance function:
##  Structure: fixed weights
##  Formula: ~invwt 
## Fixed effects:  cbind(r, rr) ~ seed * extract 
##                   Value Std.Error DF   t-value p-value
## (Intercept)  -0.5454297 0.1826058 17 -2.986925  0.0083
## seed          0.1048122 0.3058526 17  0.342689  0.7360
## extract       1.3233352 0.2586204 17  5.116901  0.0001
## seed:extract -0.7985771 0.4243540 17 -1.881866  0.0771
##  Correlation: 
##              (Intr) seed   extrct
## seed         -0.597              
## extract      -0.706  0.422       
## seed:extract  0.430 -0.721 -0.609
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.59952544 -0.85340590  0.03919434  0.71649291  1.30878355 
## 
## Number of Observations: 21
## Number of Groups: 21

ranef(glmm_pql)

##     (Intercept)
## 1  -0.133830829
## 2   0.006642181
## 3  -0.161176902
## 4   0.205529628
## 5   0.082835921
## 6   0.063874151
## 7  -0.057458067
## 8  -0.099110306
## 9   0.157250418
## 10 -0.064556196
## 11  0.036118561
## 12  0.058642575
## 13  0.147457793
## 14 -0.085440244
## 15 -0.198619343
## 16  0.041840659
## 17 -0.123583606
## 18  0.019021152
## 19 -0.020471760
## 20  0.150852971
## 21 -0.025818757

Interpretasi dari model GLMM PQL di atas:

Pada model ini tidak dihasikan nilai AIC karena pada pendekatan PQL tidak memaksimumkan likelihood marginal yang sebenarnya.
Dari nilai p terlihat bahwa extract berpengaruh signifikan terhadap peluang pertumbuhan benih.
Nilai penduga seed=0,1048122 artinya jika seed berubah dari O75 ke O73 maka nilai log-odds meningkat sebesar 0,1048122
Nilai penduga extract= 1,3233352 artinya jika extract berubah dari bean ke cucumber maka nilai log-odds meningkat sebesar 1,3233352
Nilai penduga seed:extract=-0,7985771 artinya efek interaksi antara seed=O73 dan extract=cucumber menurunkan nilai log-odds sebesar 0,7985771
Contoh prediksi: jika extract=cucumber (1) dan seed=O75(0) pada plate ke-18 maka nilai log-odds=-0,5454297+0,1048122*0+1,3233352*1-0,7985771*0*1+0,019021152=0,796927 sehingga diperoleh peluang benih tumbuh sebesar 69,93%