请选择 进入手机版 | 继续访问电脑版
楼主: 大叔肥肥
736 1

[面板数据求助] 帮忙看一下以下STATA代码是什么玩意 [推广有奖]

  • 1关注
  • 0粉丝

本科生

40%

还不是VIP/贵宾

-

威望
0
论坛币
0 个
通用积分
1.0000
学术水平
0 点
热心指数
1 点
信用等级
0 点
经验
290 点
帖子
5
精华
0
在线时间
184 小时
注册时间
2012-9-15
最后登录
2021-10-11

大叔肥肥 在职认证  发表于 2018-4-11 11:05:06 |显示全部楼层

+2 论坛币
k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群
赵安豆老师微信:zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

感谢您参与论坛问题回答

经管之家送您两个论坛币!

+2 论坛币
// Perform fixed effects analyses and all REML random effects analyses
// There are 5 variants of random effects we are fitting here with REML. They differ
// in the way the T matrix is structured (covariance matrix of the random effect)
// case A: common sigma and common correlation (2 parameters)
// case B: different sigmas and a common correlation (5 parameters)
// case C: common sigma and a banded correlation structure (2 parameters)
// case D: different sigma and banded correlation structure (5 parameters)
// case E: unstructruted (10 parameters)

use data, clear
set seed 12345

local p 4

qui count
local n = r(N)

global restricted = 1 // do REML rather than ML


// These globals and matrices are needed to pass information to the likelihood
// optimization programs.
//这些全局变量和矩阵是用于将信息传递给似然性的
//优化程序。
global n `n'
global ymat  "ymat"
global Smat  "Smat"
global p `p'

forval i=1/`n' {
        mat ymat`i' = J(1, `p', 0)       
        forval j =1/`p' {
                qui summ b`j' in `i' , meanonly                 
                mat ymat`i'[1,`j'] = r(mean)
        }
        if (matmissing(ymat`i')) {
                noi di in yellow "Study `i' has missing values:"
                noi mat li ymat`i'
        }

        mat Smat`i' = J(`p', `p' , 0)
        forval k =1/`p' {
                forval m =`k'/`p' {
                        qui summ  V`k'`m' in `i' , meanonly
                        mat Smat`i'[`k', `m'] = r(mean)
                        mat Smat`i'[`m', `k'] = r(mean)
                }
        }       
}


// Note that study 17 has missing values. Its likelihood contribution is different
// than that of the other studies.  
// Here we replace the missing values in time point 1 and 3 (6 and 18 months)
// with zeros.  
// This will work fine for the fixed effects GLS -- see below.
// The REML programs however will have to do some footwork to accommodate the 17th study

forval i=1/`p' {
        if (el(ymat17,1,`i')>=.)  mat ymat17[1,`i']=0
        forval j=1/`p' {
                if (el(Smat17,`i',`j')>=.)  mat Smat17[`i',`j']=0
        }
}


// All matrices must be non-singular. In this example, Smat14 (study 14) is singular
// We must assert that these matrices are positive definite; the fastest way
// would be to try-catch a Cholesky decomposition (See G Strang Linear Algebra)
// but here i will simply check the eigenvalues, as described in the paper.
// correctmemat is a short program that does this correction quickly.

// First do the 16 studies with complete data
forval i=1/`=`n'-1' {
        correctmemat, matname(Smat`i') epsilon(0.08)
        if (r(corrected)) mat Smat`i' = r(M)
}

// now check the 17th study for the outcomes that are non-missing
// Calculate a permutation matrix P that will permute Smat17 so that
// the nonmissing timepoints are first and the missing follow
// this is done with a small program I wrote
gimmepermmat , matname(Smat17)
mat P = r(P)  // this is the permutation matrix

mat A = P*Smat17*P'  // save rearranged matrix as A
mat A= A[1..2, 1..2] // only the nonmissing values
di det(A)             // determinant is positive, go on
                     // if it were not, we would correct the 2x2 A matrix as above,
                     // pad the corrected 2x2 with 0's to get a 4x4
                     // and restore the order of the rows and columns
                     // using P'*corrected_padded_matrix*P
mat drop A             


// This is the fixed effects model, fit with GLS. The same can be fit using ML
// but I want to show this coding also. So bear with me.
// Because i have no covariates, there is no design matrix here (it's I(4), omitted)

if (0==0) {  
        mat Wsum = J(`p', `p', 0)
        mat Wysum = J(1, `p', 0)

        // The 16 studies with all 4 time points are OK
        // turns out that for the fixed effects model,  
        // i can use the Smat17 which has 0's in the rows and columns
        // of the missing time points; and the means vector
        // ymat17, which also has 0's for the missing timepoints.
        // This works out to be the same as the method by
        // Gleser and Olkin referenced in the paper.
        // this strategy will not work OK for the REML calculations!

        forval k=1/`n' {
                mat W`k' = syminv(Smat`k')
                mat Wysum = Wysum + ymat`k' * W`k'
                mat Wsum = Wsum + W`k'
        }

        mat covbetaFixed = syminv(Wsum)
        mat betaFixed = Wysum * covbetaFixed
       
        // by definition TauFixed =0; I will use this when gathering results
        mat TauFixed = J(`p', `p', 0)
}

// This is the fixed effects model with a different coding (REML)
// The results are identical with those from the GLS approach above

if (0==0) {
        // Fixed effects model with REML
        program drop _all
        global be_verbose = 1  // force ll program to identify itself - avoid blunders

        ml model d0 ll_fixed_miss (mu: b1 b2 b3 b4 , nocons), obs(`n')  collinear

        mat b0 = [betaFixed]
        ml init b0 , copy
        //ml check

        ml max , difficult iterate(30) ltolerance(1e-6)// invalid syntax
        est store Fixed
}


// These are the five examples of random effects that are discussed in the paper
if (1==1) {
        // case A: common variance and common correlation (5 parameters)
        program drop _all
        global be_verbose = 1  // force ll program to identify itself - avoid blunders

        ml model d0 ll_caseA_miss (mu: b1 b2 b3 b4 , nocons) (S:) (rho:), obs(`n')  collinear
       
        ml search S: -10 10 rho: 0 1

        mat b0 = [betaFixed, 0, 0]
        ml init b0 , copy
        //ml check

        ml max , difficult iterate(30) ltolerance(1e-6)
        est store A
        mat betaA = BETA
        mat TauA = T
        mat covbetaA = COVBETA
}

if (2==2) {
        // case B: different variances and common correlation (5 parameters)
        program drop _all
        global be_verbose = 1  // force ll program to identify itself - avoid blunders

        ml model d0 ll_caseB_miss (mu: b1 b2 b3 b4 , nocons) (S1:) (S2:) (S3:) (S4:) (rho:), ///
                obs(`n')  collinear

        ml search S1: -10 10 S2: -10 10 S3: -10 10 S4: -10 10 rho: 0 1
       
        mat b0 = [betaFixed, 1, 1,1,1,0]
        ml init b0 , copy skip
        //ml check

        ml max , difficult iterate(30) ltolerance(1e-6)
        est store B
        mat betaB = BETA
        mat TauB = T
        mat covbetaB = COVBETA
}

if (3==3) {
        // case C: common variance - autoregressive correlation (2 parameters)
        program drop _all
        global be_verbose = 1  // force ll program to identify itself - avoid blunders

        ml model d0 ll_caseC_miss (mu: b1 b2 b3 b4 , nocons) (S:) (rho:) , obs(`n')  collinear
       
        ml search S: -10 10  rho: 0 1

        mat b0 = [betaFixed, 0, 0]
        ml init b0 , copy skip
        //ml check

        ml max , difficult iterate(30) ltolerance(1e-6)
        est store C
        mat betaC = BETA
        mat TauC = T
        mat covbetaC = COVBETA
}

if (4==4) {
        // case D: different variances - autoregressive correlation (5 parameters)
        program drop _all
        global be_verbose = 1  // force ll program to identify itself - avoid blunders

        ml model d0 ll_caseD_miss (mu: b1 b2 b3 b4 , nocons) (S1:) (S2:) (S3:) (S4:) (rho:), ///
                obs(`n')  collinear

        ml search S1: -10 10 S2: -10 10 S3: -10 10 S4: -10 10 rho: 0 1
       
        mat b0 = [betaFixed, 1, 1,1,1,0]
        ml init b0 , copy skip
        //ml check

        ml max , difficult iterate(30) ltolerance(1e-6)
        est store D
        mat betaD = BETA
        mat TauD = T
        mat covbetaD = COVBETA
}

if (5==5) {
        // case E: completely unstructured (10 parameters)
        program drop _all
        global be_verbose = 1  // force ll program to identify itself - avoid blunders

        ml model d0 ll_caseE_miss (mu: b1 b2 b3 b4 , nocons) ///
                (S11:) (S12:) (S13:) (S14:) ///
                       (S22:) (S23:) (S24:) ///
                                     (S33:) (S34:) ///
                                     (S44:) , obs(`n')  collinear

        ml search S11: -10 10 S12: -10 10 S13: -10 10 S14: -10 10 ///
                                S22: -10 10 S23: -10 10 S24: -10 10 ///
                                              S33: -10 10 S34: -10 10 ///
                                                      S44: -10 10

       
        mat b0 = [betaFixed, 1,1,1,1,1,1,1,1,1,1]
        ml init b0 , copy skip
        //ml check

        ml max , difficult iterate(30) ltolerance(1e-6)
        est store E
        mat betaE = BETA
        mat TauE = T
        mat covbetaE = COVBETA
}


// Use the betaX, covbetaX and TauX matrices from the multivariate models
// and reports the results.

tempfile Results2
tempname multi
postfile  `multi' sorter logor se Q df str20 ( timepoint fix_ran uni_multi ) using `Results2', replace

local k = 0
foreach X in Fixed A B C D E {
       
        // calculate Q per model
        mat qmat = J(1, 1, 0)
        forval j=1/`n' {
                mat qmat = qmat + (ymat`j'-beta`X')*syminv(Smat`j' + Tau`X') * ///
                        (ymat`j' - beta`X')'
        }

        forval i =1/`p' {
                local k = `k' +1

                // the contrast matrix
                mat a = J(1, `p', 0)
                mat a[1, `i'] = 1

                // the effect size and its variance
                mat ES  = a*beta`X''
                mat VAR  = a*(covbeta`X')*a'

                post `multi' (`k') (`=el(ES, 1, 1)') ///
                (`=sqrt(el(VAR,1,1))') (`=el(qmat,1,1)') ///
                        (`=`p'*(`n'-1)') ("`=`i'*6' months") ("`X'") ("multivariate")
               
        }
}

postclose `multi'


use Results, clear // we saved the univariate results here
append using `Results2'
save Results, replace


二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

关键词:Stata tata fixed effect Multivariate Optimization

stata SPSS
Zyzjerryzhang 在职认证  学生认证  发表于 2018-9-24 21:34:40 |显示全部楼层
大姐,你这个我没看到你问题呀?

使用道具

您需要登录后才可以回帖 登录 | 我要注册

京ICP备16021002-2号 京B2-20170662号 京公网安备 11010802022788号 论坛法律顾问:王进律师 知识产权保护声明   免责及隐私声明

GMT+8, 2021-10-20 05:01