Keyword: 制限つき最尤法, REML, 線形混合効果回帰
概要
本サンプルは制限つき最尤法（REML）を用いた線形混合効果回帰のフィットを行うFortranによるサンプルプログラムです。 本サンプルは以下に示されるデータについて線形混合効果回帰のフィットを行います。
※本サンプルはnAG Fortranライブラリに含まれるルーチン g02jaf() のExampleコードです。本サンプル及びルーチンの詳細情報は g02jaf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はg02jaf のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
| このデータをダウンロード |
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 のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
| この出力例をダウンロード |
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ライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
| このソースコードをダウンロード |
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