99999久久久久久亚洲,欧美人与禽猛交狂配,高清日韩av在线影院,一个人在线高清免费观看,啦啦啦在线视频免费观看www

熱線電話：13121318867

登錄
首頁精彩閱讀R語言之縱向數(shù)據(jù)分析：多級線性增長模型2
R語言之縱向數(shù)據(jù)分析：多級線性增長模型2
2016-10-17
收藏

R語言之縱向數(shù)據(jù)分析：多級線性增長模型2

這篇文章已經(jīng)是縱向數(shù)據(jù)分析系列的第三篇了。之前，我們介紹了什么是縱向數(shù)據(jù)，我們?nèi)绾伟验L型的數(shù)據(jù)集轉(zhuǎn)換成縱向的數(shù)據(jù)集，并通過建立相應(yīng)的多級模型進(jìn)行分析。顯然，僅僅介紹基本的多級模型并不足以對虛擬隨機(jī)對照實(shí)驗(yàn)數(shù)據(jù)集進(jìn)行分析。這篇文章則延續(xù)之前的分析，并介紹一種不錯(cuò)的方法來處理多級模型里的治療效果問題。

創(chuàng)建一個(gè)縱向數(shù)據(jù)集并把它轉(zhuǎn)成長型格式

和往常一樣，我們還是先創(chuàng)建一個(gè)縱向數(shù)據(jù)集，然后把它轉(zhuǎn)成長型的格式。

  1. library(MASS)
  3. dat.tx.a<-mvrnorm(n=250,mu=c(30,20,28),
  4. Sigma=matrix(c(25.0,17.5,12.3,
  5. 17.5,25.0,17.5,
  6. 12.3,17.5,25.0),nrow=3,byrow=TRUE))
  8. dat.tx.b<-mvrnorm(n=250,mu=c(30,20,22),
  9. Sigma=matrix(c(25.0,17.5,12.3,
  10. 17.5,25.0,17.5,
  11. 12.3,17.5,25.0),nrow=3,byrow=TRUE))
  13. dat<-data.frame(rbind(dat.tx.a,dat.tx.b))
  14. names(dat)<-c('measure.1','measure.2','measure.3')
  16. dat<-data.frame(subject.id=factor(1:500),tx=rep(c('A','B'),each=250),dat)
  18. rm(dat.tx.a,dat.tx.b)
  20. dat<-reshape(dat,varying=c('measure.1','measure.2','measure.3'),
  21. idvar='subject.id',direction='long')

多級增長模型分析治療效果

這是一個(gè)RCT數(shù)據(jù)集，也就說明，這里潛在存兩個(gè)以療效分組的差異。上一次，我們沒有考慮異質(zhì)性，而僅僅對兩個(gè)組的普通療效進(jìn)行了細(xì)分。從直覺上看，我們可以添加治療變量來修復(fù)這個(gè)模型，從而抓到組間的差異。

  1. library(lmerTest)
  2. m1<-lmer(measure~time+tx+(1|subject.id),data=dat)

這里，我們?nèi)稳皇褂胠merTest包是因?yàn)樗试S修正模型的檢驗(yàn)使用自由度。模型的公式和是和上次的一樣，除非我們把治療變量tx作為獨(dú)立變量。

  1. summary(m1)
  2. Linearmixed model fitbyREML t-testsuseSatterthwaiteapproximations to degrees of freedom[
  3. lmerMod]
  4. Formula:measure~time+tx+(1|subject.id)
  5. Data:dat
  7. REML criterion at convergence:9700.8
  9. Scaledresiduals:
  10. Min1QMedian3QMax
  11. -3.03326-0.656180.080810.658932.64555
  13. Randomeffects:
  14. GroupsNameVarianceStd.Dev.
  15. subject.id(Intercept)6.8342.614
  16. Residual31.9405.652
  17. Numberof obs:1500,groups:subject.id,500
  19. Fixedeffects:
  20. EstimateStd.Errordf t valuePr(>|t|)
  21. (Intercept)30.61470.44461494.000068.856<2e-16***
  22. time-2.34340.1787999.0000-13.112<2e-16***
  23. txB-2.29110.3740498.0000-6.1271.82e-09***
  24. ---
  25. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
  27. CorrelationofFixedEffects:
  28. (Intr)time
  29. time-0.804
  30. txB-0.4210.000

從上述的summary里的信息可以很明顯的看到，time和tx這兩個(gè)變量都相當(dāng)顯著，而接受B治療的人群在抑郁癥得分上比接受A治療的人群低2.34分。那么，我們是否可以認(rèn)為，B治療方案效果更好？絕對不是。模型m1實(shí)際上并不合適對縱向數(shù)據(jù)集里的組間差異進(jìn)行比較。m1里的治療效果，其實(shí)也就是一段時(shí)間內(nèi)的平均治療效果。換句話說，治療方案A和治療方案B的療效所導(dǎo)致的分差在-2.34分，其實(shí)只不過是1,2,3時(shí)間段內(nèi)抑郁癥均分的差。這是RCT，所以，我們所期待的志愿者參與兩組治療方案里時(shí)間段1的差要小于1。我們對兩組療效的均值差異感興趣，而現(xiàn)在的軌跡在時(shí)間段2、3里存在差異。
這個(gè)概率可能對初學(xué)者來說有點(diǎn)復(fù)雜，但是其執(zhí)行過程很簡單，我們只需要添加交互變量time和tx。

  1. m2<-lmer(measure~time*tx+(1|subject.id),data=dat)

在R，我們可以使用:符號來表明變量間的交互關(guān)系。因此，我們可以這樣設(shè)置：time+tx+time:tx。但是，更簡短的寫法是：timetx。操作符要R涵蓋主要變量和治療療效。此時(shí)，使用操作符就能簡單實(shí)現(xiàn)各個(gè)效應(yīng)的展示。請記住，操作符也可用于高維變量交互。比如，ABC就要R包含三個(gè)主要的效應(yīng)，三個(gè)兩兩交互的效應(yīng)，以及一個(gè)三元交互效應(yīng)。

  1. summary(m2)
  2. Linearmixed model fitbyREML t-testsuseSatterthwaiteapproximations to degrees of freedom[
  3. lmerMod]
  4. Formula:measure~time*tx+(1|subject.id)
  5. Data:dat
  7. REML criterion at convergence:9640.8
  9. Scaledresiduals:
  10. Min1QMedian3QMax
  11. -3.0454-0.67760.10410.66222.4296
  13. Randomeffects:
  14. GroupsNameVarianceStd.Dev.
  15. subject.id(Intercept)7.4482.729
  16. Residual30.1005.486
  17. Numberof obs:1500,groups:subject.id,500
  19. Fixedeffects:
  20. EstimateStd.Errordf t valuePr(>|t|)
  21. (Intercept)27.88080.55741421.500050.017<2e-16***
  22. time-0.97650.2454998.0000-3.9807.40e-05***
  23. txB3.17670.78831421.50004.0305.88e-05***
  24. time:txB-2.73390.3470998.0000-7.8798.44e-15***
  25. ---
  26. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
  28. CorrelationofFixedEffects:
  29. (Intr)time txB
  30. time-0.880
  31. txB-0.7070.622
  32. time:txB0.622-0.707-0.880

現(xiàn)在，我們從m2模型中獲得一個(gè)三元交互關(guān)系。在這個(gè)模型當(dāng)中，time效應(yīng)指的就是療效A（這是特定條件下的療效和療效時(shí)間的交互），而txB效應(yīng)就是時(shí)間點(diǎn)為0時(shí)的效應(yīng)（數(shù)據(jù)集里一個(gè)虛的時(shí)間點(diǎn)）。這個(gè)術(shù)語解釋起來有點(diǎn)復(fù)雜：系數(shù)值反映了每個(gè)時(shí)間段內(nèi)兩種治療方案的療效上的差異，其條件是取平均時(shí)間和治療效果?？陬^上說，人們可以理解的是，這表明兩組治療方案是如何隨著時(shí)間的推移而不同。對于每個(gè)單位時(shí)間內(nèi)的增長（即，從時(shí)間段1到時(shí)間段2），在完成平均治療效應(yīng)差的分析以后，接受B治療方案的參與者理論上在抑郁癥得分方面比接受A治療方案的要低2.88分。所以，在時(shí)間段2內(nèi)，接受B治療方案的參與者得分差是 -2.88 x 2 + 3.43 = -2.33，而在時(shí)間段3內(nèi)，接受B治療方案的參與者得分差是-2.88 x 3+ 3.43 = -5.21。我們應(yīng)當(dāng)時(shí)刻注意，在一個(gè)給定的時(shí)間段里對比療效差的時(shí)候，主要療效以及療效變量間的交互。但是，為什么要這樣？你不妨設(shè)置分析一下，當(dāng)你忽略了主要治療效應(yīng)的時(shí)候，你在時(shí)間段1內(nèi)所得到的療效差是怎樣的情形。

選擇最后的模型

我可以告訴你，模型m2比m1和m0（上一篇文章所采用的模型）都要好，原因在于，這個(gè)數(shù)據(jù)集是我們自己創(chuàng)建的。但在實(shí)際生活中，我們會在選擇最好的一個(gè)模型的時(shí)候會遇到很多困難。其中一種選擇的方式就是，按照我們之前所講的，使用summary()函數(shù)，以此算出額外變量的統(tǒng)計(jì)顯著性。其它一個(gè)常用的方法就是使用方差分析（ANOVA)表（僅僅針對嵌套模型），而這種方法適合用于這個(gè)模型。

  1. m0<-lmer(measure~time+(1|subject.id),data=dat)
  2. anova(m0,m1,m2)
  4. Data:dat
  5. Models:
  6. object:measure~time+(1|subject.id)
  7. ..1:measure~time+tx+(1|subject.id)
  8. ..2:measure~time*tx+(1|subject.id)
  9. DfAIC BIC logLik devianceChisqChiDfPr(>Chisq)
  10. object49741.99763.1-4866.99733.9
  11. ..159707.69734.1-4848.89697.636.33411.663e-09***
  12. ..269649.29681.1-4818.69637.260.34517.962e-15***
  13. ---
  14. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

這里，我們使用anova()對m0、m1和m2模型進(jìn)行方差分析，原因在于我們把m0嵌套于m1，然后把m1嵌套于m2。如果還有另外的子集，我們就把模型嵌套到那個(gè)子集上。換句話說，含有獨(dú)立變量A、B的模型則嵌套于含有變量A、B、C的模型中。
如果你做過方差分析，那這個(gè)方差分析表對于你來說應(yīng)該是很熟悉了。Df列表明了模型里的自由度，它簡單的反映了這個(gè)案例的參數(shù)估計(jì)。例如，自由度為4的模型m0表明了這個(gè)模型預(yù)測了4個(gè)參數(shù)（1.修正截距效應(yīng)系數(shù)，2.時(shí)間效應(yīng)，3.隨機(jī)截距方差，4.殘差）。
AIC和BIC是兩個(gè)常用的擬合指數(shù)。AIC和BIC都考慮到兩個(gè)因素：模型的數(shù)據(jù)擬合效果有多好，和這個(gè)模型復(fù)不復(fù)雜。AIC和BIC的不同之處在于模型復(fù)雜度的罰分上（BIC對模型復(fù)雜度的罰分影響更大）。我發(fā)現(xiàn)，AIC理解起來更容易，因?yàn)樗鼭u進(jìn)的使用留一交叉檢驗(yàn)法（LOOCV)來預(yù)測運(yùn)行效果。我們應(yīng)當(dāng)理解交叉檢驗(yàn)和預(yù)測性能的概念，隨后，如果你對此不太理解，你可以先把它跳過。

LogLik列模型參數(shù)的似然對數(shù)。從根本上說，它是模型里數(shù)據(jù)擬合效果的指標(biāo)。其偏差是?2×logLik?2×logLik。為什么我們需要它？因?yàn)閮蓚€(gè)嵌套模型的偏差在原假設(shè)的條件下遵循卡方檢驗(yàn)（假設(shè)兩個(gè)嵌套模型的偏差一樣）。此時(shí)，我們可以把它用作模型擬合差的測試。偏差的差則在Chisq列顯示?？ǚ綑z驗(yàn)的統(tǒng)計(jì)相關(guān)自由度則在Chi Df（簡單來說，就是額外參數(shù)的個(gè)數(shù)）。最后一行的結(jié)果顯而易見：它算出了測試兩個(gè)嵌套模型的差的p值。雖然是這樣，我們還是要警惕ANOVA的結(jié)果，因?yàn)樵诒姸嗟暮蜻x模型中，它會導(dǎo)致錯(cuò)誤類型1上升極快。

由于模型m2比其它兩個(gè)候選模型AIC和BIC都要好，而對于方差分析，我們可以認(rèn)為這個(gè)模型是目前為止我們擬合的最好的了。當(dāng)然，這個(gè)模型還有許多要改進(jìn)的地方。后續(xù)，我們可能還可以在模型的性能上進(jìn)行提升，并作圖進(jìn)行預(yù)測。

數(shù)據(jù)分析咨詢請掃描二維碼

若不方便掃碼，搜微信號：CDAshujufenxi

數(shù)據(jù)分析師資訊
更多

OK
客服在線
立即咨詢
客服在線
立即咨詢
') } function initGt() { var handler = function (captchaObj) { captchaObj.appendTo('#captcha'); captchaObj.onReady(function () { $("#wait").hide(); }).onSuccess(function(){ $('.getcheckcode').removeClass('dis'); $('.getcheckcode').trigger('click'); }); window.captchaObj = captchaObj; }; $('#captcha').show(); $.ajax({ url: "/login/gtstart?t=" + (new Date()).getTime(), // 加隨機(jī)數(shù)防止緩存 type: "get", dataType: "json", success: function (data) { $('#text').hide(); $('#wait').show(); // 調(diào)用 initGeetest 進(jìn)行初始化 // 參數(shù)1：配置參數(shù) // 參數(shù)2：回調(diào)，回調(diào)的第一個(gè)參數(shù)驗(yàn)證碼對象，之后可以使用它調(diào)用相應(yīng)的接口 initGeetest({ // 以下 4 個(gè)配置參數(shù)為必須，不能缺少 gt: data.gt, challenge: data.challenge, offline: !data.success, // 表示用戶后臺檢測極驗(yàn)服務(wù)器是否宕機(jī) new_captcha: data.new_captcha, // 用于宕機(jī)時(shí)表示是新驗(yàn)證碼的宕機(jī) product: "float", // 產(chǎn)品形式，包括：float，popup width: "280px", https: true // 更多配置參數(shù)說明請參見：http://docs.geetest.com/install/client/web-front/ }, handler); } }); } function codeCutdown() { if(_wait == 0){ //倒計(jì)時(shí)完成 $(".getcheckcode").removeClass('dis').html("重新獲取"); }else{ $(".getcheckcode").addClass('dis').html("重新獲取("+_wait+"s)"); _wait--; setTimeout(function () { codeCutdown(); },1000); } } function inputValidate(ele,telInput) { var oInput = ele; var inputVal = oInput.val(); var oType = ele.attr('data-type'); var oEtag = $('#etag').val(); var oErr = oInput.closest('.form_box').next('.err_txt'); var empTxt = '請輸入'+oInput.attr('placeholder')+'！'; var errTxt = '請輸入正確的'+oInput.attr('placeholder')+'！'; var pattern; if(inputVal==""){ if(!telInput){ errFun(oErr,empTxt); } return false; }else { switch (oType){ case 'login_mobile': pattern = /^1[3456789]\d{9}$/; if(inputVal.length==11) { $.ajax({ url: '/login/checkmobile', type: "post", dataType: "json", data: { mobile: inputVal, etag: oEtag, page_ur: window.location.href, page_referer: document.referrer }, success: function (data) { } }); } break; case 'login_yzm': pattern = /^\d{6}$/; break; } if(oType=='login_mobile'){ } if(!!validateFun(pattern,inputVal)){ errFun(oErr,'') if(telInput){ $('.getcheckcode').removeClass('dis'); } }else { if(!telInput) { errFun(oErr, errTxt); }else { $('.getcheckcode').addClass('dis'); } return false; } } return true; } function errFun(obj,msg) { obj.html(msg); if(msg==''){ $('.login_submit').removeClass('dis'); }else { $('.login_submit').addClass('dis'); } } function validateFun(pat,val) { return pat.test(val); }