Keyword: 制限つき最尤法, REML, 線形混合効果回帰
概要
本サンプルは制限つき最尤法（REML）を用いた線形混合効果回帰のフィットを行うFortranによるサンプルプログラムです。 本サンプルは以下に示されるデータについて線形混合効果回帰のフィットを行います。
※本サンプルはnAG Fortranライブラリに含まれるルーチン g02jaf() のExampleコードです。本サンプル及びルーチンの詳細情報は g02jaf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はg02jaf のマニュアルページを参照）| このデータをダウンロード |
G02JAF Example Program Data 24 5 3 1 1 1 4 3 2 3 1 3 4 5 3 2 0 1 1 1 56 1 1 1 1 50 1 2 1 1 39 1 3 1 1 30 2 1 1 1 36 2 2 1 1 33 2 3 1 1 32 3 1 1 1 31 3 2 1 1 15 3 3 1 1 30 4 1 1 1 35 4 2 1 1 17 4 3 1 1 41 1 1 2 1 36 1 2 2 2 35 1 3 2 3 25 2 1 2 1 28 2 2 2 2 30 2 3 2 3 24 3 1 2 1 27 3 2 2 2 19 3 3 2 3 25 4 1 2 1 30 4 2 2 2 18 4 3 2 3 1.0 1.0 -1
- １行目はタイトル行で読み飛ばされます。
- ２行目は観測値の数（n=24）、データ行列のカラム数（ncol=5）、固定として処理される独立変数の数（nfv=3）、ランダムとして処理される独立変数の数（nrv=1）、分散成分の数（nvpr=1）を指定しています。
- ３行目は各変数のレベル数(levels)を指定しています。
- ４行目は従属変数をもつデータ行列DATのカラム（yvid=1）、固定独立変数をもつDATのカラム（fvid=3,4,5）、ランダム独立変数をもつDATのカラム（rvid=3）、subject（個体）変数をもつDATのカラム（svid=2）、case weightをもつDATのカラム（cwid=0:無し）、固定切片が含まれるかどうか（fint=1:含まれる）、ランダム切片が含まれるかどうか（rint=1:含まれる）を指定しています。
- ５行目はランダム変数の分散を示すフラグ(vpr)を指定しています。
- ６〜２９行目はデータ行列(dat)を指定しています。
- ３０行目はガンマの初期値(gamma)を指定しています
- ３１行目は反復の最大値(maxit)を指定しています。ゼロより小さい場合は100が使用されます。
出力結果
（本ルーチンの詳細はg02jaf のマニュアルページを参照）| この出力例をダウンロード |
G02JAF Example Program Results
Fixed effects (Estimate and Standard Deviation)
Intercept 37.0000 4.6674
Variable 1 Level 2 1.0000 3.5173
Variable 1 Level 3 -11.0000 3.5173
Variable 2 Level 2 -8.2500 2.1635
Variable 3 Level 2 0.5000 3.0596
Variable 3 Level 3 7.7500 3.0596
Random Effects (Estimate and Standard Deviation)
Intercept for Subject Level 1 10.7631 4.4865
Subject Level 1 Variable 1 Level 1 3.7276 3.0331
Subject Level 1 Variable 1 Level 2 -1.4476 3.0331
Subject Level 1 Variable 1 Level 3 0.3733 3.0331
Intercept for Subject Level 2 -0.5269 4.4865
Subject Level 2 Variable 1 Level 1 -3.7171 3.0331
Subject Level 2 Variable 1 Level 2 -1.2253 3.0331
Subject Level 2 Variable 1 Level 3 4.8125 3.0331
Intercept for Subject Level 3 -5.6450 4.4865
Subject Level 3 Variable 1 Level 1 0.5903 3.0331
Subject Level 3 Variable 1 Level 2 0.3987 3.0331
Subject Level 3 Variable 1 Level 3 -2.3806 3.0331
Intercept for Subject Level 4 -4.5912 4.4865
Subject Level 4 Variable 1 Level 1 -0.6009 3.0331
Subject Level 4 Variable 1 Level 2 2.2742 3.0331
Subject Level 4 Variable 1 Level 3 -2.8052 3.0331
Variance Components
1 62.3958
2 15.3819
SIGMA^2 = 9.3611
-2LOG LIKE = 119.7618
DF = 16
- ３〜１０行目に固定効果が出力されています。
- ５行目に固定切片と標準誤差が出力されています。
- ６〜１０行目に固定独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- １２〜２９行目にランダム効果が出力されています。
- １４行目にsubject変数のレベル１のランダム切片と標準誤差が出力されています。
- １５〜１７行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- １８行目にsubject変数のレベル２のランダム切片と標準誤差が出力されています。
- １９〜２１行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- ２２行目にsubject変数のレベル３のランダム切片と標準誤差が出力されています。
- ２３〜２５行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- ２６行目にsubject変数のレベル４のランダム切片と標準誤差が出力されています。
- ２７〜２９行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- ３１〜３３行目に分散成分が出力されています。
- ３５行目にシグマの２乗が出力されています。
- ３６行目に対数尤度（log-likelihood）が出力されています。
- ３７行目に自由度が出力されています。
ソースコード
（本ルーチンの詳細はg02jaf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM g02jafe
! G02JAF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : g02jaf, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
REAL (KIND=nag_wp) :: reml, tol
INTEGER :: cwid, df, fint, i, ifail, j, k, l, &
lb, lddat, maxit, n, ncol, nff, nfv, &
nrf, nrv, nvpr, rint, svid, warn, yvid
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: b(:), dat(:,:), gamma(:), se(:)
INTEGER, ALLOCATABLE :: fvid(:), levels(:), rvid(:), vpr(:)
! .. Intrinsic Functions ..
INTRINSIC max
! .. Executable Statements ..
WRITE (nout,*) 'G02JAF Example Program Results'
WRITE (nout,*)
! Skip heading in data file
READ (nin,*)
! Read in the problem size
READ (nin,*) n, ncol, nfv, nrv, nvpr
ALLOCATE (levels(ncol),fvid(nfv),rvid(nrv))
! Read in number of levels for each variable
READ (nin,*) levels(1:ncol)
! Read in model information
READ (nin,*) yvid, fvid(1:nfv), rvid(1:nrv), svid, cwid, fint, rint
! If no subject specified, then ignore RINT
IF (svid==0) THEN
rint = 0
END IF
! Calculate LB
lb = rint
DO i = 1, nrv
lb = lb + levels(rvid(i))
END DO
IF (svid/=0) THEN
lb = lb*levels(svid)
END IF
lb = lb + fint
DO i = 1, nfv
lb = lb + max(levels(fvid(i))-1,1)
END DO
lddat = n
ALLOCATE (vpr(nrv),dat(lddat,ncol),gamma(nvpr+2),b(lb),se(lb))
! Read in the variance component flag
READ (nin,*) vpr(1:nrv)
! Read in the Data matrix
READ (nin,*) (dat(i,1:ncol),i=1,n)
! Read in the initial values for GAMMA
READ (nin,*) gamma(1:(nvpr+rint))
! Read in the maximum number of iterations
READ (nin,*) maxit
! Use default value for tolerance
tol = 0.0E0_nag_wp
! Fit the linear mixed effects regresion model
ifail = 0
CALL g02jaf(n,ncol,lddat,dat,levels,yvid,cwid,nfv,fvid,fint,nrv,rvid, &
nvpr,vpr,rint,svid,gamma,nff,nrf,df,reml,lb,b,se,maxit,tol,warn, &
ifail)
! Display results
IF (warn/=0) THEN
WRITE (nout,*) 'Warning: At least one variance component was ', &
'estimated to be negative and then reset to zero'
WRITE (nout,*)
END IF
WRITE (nout,*) 'Fixed effects (Estimate and Standard Deviation)'
WRITE (nout,*)
k = 1
IF (fint==1) THEN
WRITE (nout,99999) 'Intercept ', b(k), se(k)
k = k + 1
END IF
DO i = 1, nfv
DO j = 1, levels(fvid(i))
IF (levels(fvid(i))==1 .OR. j/=1) THEN
WRITE (nout,99995) 'Variable', i, ' Level', j, b(k), se(k)
k = k + 1
END IF
END DO
END DO
WRITE (nout,*)
WRITE (nout,*) 'Random Effects (Estimate and Standard', ' Deviation)'
WRITE (nout,*)
IF (svid==0) THEN
DO i = 1, nrv
DO j = 1, levels(rvid(i))
WRITE (nout,99995) 'Variable', i, ' Level', j, b(k), se(k)
k = k + 1
END DO
END DO
ELSE
DO l = 1, levels(svid)
IF (rint==1) THEN
WRITE (nout,99998) 'Intercept for Subject Level', l, &
' ', b(k), se(k)
k = k + 1
END IF
DO i = 1, nrv
DO j = 1, levels(rvid(i))
WRITE (nout,99997) 'Subject Level', l, ' Variable', i, &
' Level', j, b(k), se(k)
k = k + 1
END DO
END DO
END DO
END IF
WRITE (nout,*)
WRITE (nout,*) ' Variance Components'
WRITE (nout,99996) (i,gamma(i),i=1,nvpr+rint)
WRITE (nout,*)
WRITE (nout,99994) 'SIGMA^2 = ', gamma(nvpr+rint+1)
WRITE (nout,99994) '-2LOG LIKE = ', reml
WRITE (nout,99993) 'DF = ', df
99999 FORMAT (1X,A,2F10.4)
99998 FORMAT (1X,A,I4,A,2F10.4)
99997 FORMAT (1X,3(A,I4),2F10.4)
99996 FORMAT (1X,I4,F10.4)
99995 FORMAT (1X,2(A,I4),2F10.4)
99994 FORMAT (1X,A,F10.4)
99993 FORMAT (1X,A,I16)
END PROGRAM g02jafe