医咖会stata 笔记（自己能看懂版
感謝brain 老師和 所有跟制作教程有關 評論區熱心解答問題 的人
?
results：輸入的代碼+結果? 如果運行出錯，也會顯示
Command：輸入代碼 regrets y x 課程之后會用do file 可以替代 command
Review：回顧之前的操作? 第一欄command 之前用過的代碼
第二欄RC error code 如果之前的代碼錯了，RC 會是紅色
Variables：變量
Properties：變量的類型是數值型、邏輯型
Note：還有注意的
Log：將所有命令和結果記錄
Viewer：包含了所有的help file，
Graph：對繪制的圖像進行更改
Do-file：可以將代碼保存至此，具有可重復性
第三講 初步數據觀測
拿到數據集之后：
有多少個變量-列，觀測值-行
符合條件的觀測值有多少個？
一人一條數據？一人多條數據？
Describe：描述 count：計數 isid：是不是id unique：特殊的
選擇內置數據集 1978 automobile data
[command]? sysuse auto, clear
1.describe 種類（數值型），格式（保留幾位小數），標簽
syntax（語法）：describe [varlist] 變量名單=一個/多個變量
attention：[] 意味著 varlist 不是必須的，可以只輸入 describe
如果只想看數據集的整體信息，不想看詳細的 [command] describe, short
看每個變量的信息 [command] describe price
2.count 直接告訴該數據集有多少行rows?? observations
syntax：count
可以與logic 連用
syntax：count if
3.isid (is ID or not?) ID 是獨特的，可以區分每一行觀測值的；是否是這個數據的可識別標簽
syntax：isid varlist 沒有[ ] 因此 varlist是必須輸入的內容
如果數據集，每人只有一個觀測值，通過ID（每人的）就能區分每一行是哪個人的信息，若一人有兩條或多條觀測值，只看ID是區分不了一人的每一行數據，只能知道這些一行屬于一個人。此時ID不是獨特的，可以區分每一行數據的值
意思是有些觀測值的mpg是相同的
mpg 和 weight相結合 是unique 的 identifier
也沒有報錯，證明大家的 price都是獨特的
4.unique
盜版的 java installation not found 外部命令手動安裝
網址“Index of /RePEc/bocode/u (bc.edu)”
ado pkg sthlp hlp 四種格式的后綴 如果有hlp 盡量把4個全部下載；沒有hlp 前三個即可，【將鏈接另存為】
然后將【保存的鏈接】挪到【D:\stata17\ado\base\u】對應的首字母目錄下，重啟stata 即可
syntax：unique varlist
74 record 中 只有 21 value，所以 not uniquely identify，isid mpg就會報錯
第四講 統計描述指標1
上節課：變量數量、觀測值數量、觀測值篩選、具有獨特值的數據量 ?describe、count、isid、unique
sysuse auto, clear
1.codebook 數據字典
一般拿到一個新的data set后，看感興趣的變量的數據字典
overview of variable type, stats, number of missing/ unique values 類型，統計量，缺失值，特異值
syntax: codebook [varlist] [if] [in] [,options] ??[ ]內的不是必須輸入項目
codebook該數據集中所有變量的信息
varlist：變量名單=一個或多個變量 if：邏輯判斷 in：第幾到第幾個觀測值 options：跟在逗號后面，一些可以自定義的選項
codebook price
codebook price if price >5000
codebook price in 10/20
in 10 ??看第10個觀測值的codebook
in 10/30 …第10到30個…
in 10/l? …第10到last（最后一個）觀測值…
in f/10? …第1（first）到10…
help codebook
codebook, mv
不過一般codebook很少用到options，因為一般只是初步看看觀測值
2.summarize
print summary statistics（mean，stdev，min，max）for variables
syntax：summarize [varlist] [if] [in] [weight] [, options]
推薦使用 sum 或 sum，完全等同于summarize
summarize ( sum / sum) price
detail 看一些額外的統計量，例子見下
sum price, detail （【detail】跟在【，】后面因為它是option）
| 偏度 豐度 |
| 最便宜的，第二便宜的… |
第五講 統計描述指標2
histogram 直方圖 graph box（箱形圖；箱形圖box-and-whisker plot） violin plot 小提琴圖
sysuse auto, clear
1.histogram 直方圖
varname：不是varlist，因此只能有一個變量
histogram可以簡寫成hist（觀察下劃線
syntax: histogram varname [if] [in] [weight] [, [continuous_opts | discrete_opts] options]
hist price
分成了8組，縱軸默認是density數據的密度，
hist price, freq?????????????????????? hist price, percent?????????????? hist price, frac
?
frequency 頻數 percent 百分比 fraction 小數（也就是percent/100 總數是1，而不是100%）
可以改變bin的寬度（Bins 是您想要將所有數據劃分成的間隔數，以便它可以在直方圖上顯示為條形。）
continuous 連續性變量，指定多少個bin，指定width寬度
discrete 離散型變量
hist price, freq bin(5) 設定5個bin，stata自動計算每個bin的寬度
還可以加一些密度曲線，如果樣本量足夠大，組距就會變得很短，階梯折線變為曲線
hist price, freq bin(5) normal 擬合正態曲線
hist price, by(foreign) 按照foreign分組繪制
2.graph box 或者 graph hbox 箱線圖（horizontal box橫著的箱子）
syntax：graph box yvars [if] [in] [weight] [, options] ??graph hbox yvars [if] [in] [weight] [, options]
【復習】箱線圖中，median 中位數，upper hinge是75百分位數，lower hinge是25百分位數
upper adjacent value為upper hinge+1.5*IQR（IQR=75百分位數-25百分位數）
lower adjacent value 為lower hinge-1.5*IQR
若數據落在upper/lower adjacent value之外，稱之為outside values離群值
graph box price????????????????????????? ?graph hbox price
graph box price, over(foreign) 通過foreign分層
3.vioplot 小提琴圖（不是自帶的 需要安裝
vioplot price
| 75百分位有這么多觀測值 |
| 25百分位有這么多觀測值 |
vioplot price, over (foreign)
第六講 散點圖 scatter plots
調整stata的主題
graph query, schemes
推薦的主題 s1mono 干凈整潔適合雜志發表
set scheme s1mono （stata重啟前設置為s1mono主題
set scheme s1mono, perm （永久設置為s1mono主題 permanent 永久的
1.scatter plot散點圖
[twoway] scatter varlist [if] [in] [weight] [, options]
twoway -- xy兩個變量繪制圖 ?
varlist -- 最基礎的形式 twoway scatter y x
進階形式 twoway scatter y1 y2 y3 … x（x與不同y變量之間的關系
twoway scatter mpg weight
twoway scatter weight length price
改變散點圖的形狀、顏色、大小等等
twoway scatter mpg weight, msymbol( ) mcolor( ) msize( )
msymbol 改變形狀（help symbolstyle
mcolor 改變顏色（help colorstyle
msize 改變大小（help markersizestyle
twoway scatter mpg weight, msymbol(D) mcolor(blue) msize(medium)
option：by 按照某分類變量分別繪制散點圖
twoway scatter mpg weight, by(foreign)
想把domestic 和 foreign 整合在一張圖上
twoway(scatter mpg weight if foreign==0)(scatter mpg weight if foreign==1),legend(label(1"Domestic")label(2"Foreign")) 報錯 原因1、2后面沒寫空格
twoway (scatter mpg weight if foreign==0)(scatter mpg weight if foreign==1), legend(label(1 "Domestic") label(2 "Foreign"))
if foreign==0 -- 國內生產的車domestic
每個scatter 用括號括了起來，是一個單獨的plot，都跟在twoway后面，合并到一起
legend 設置兩個label，第一個scatter plot設置為domestic，第二個是foreign
第七講 雙變量圖像2
sysuse uslifeexp, clear
help twoway 查看 stata可以繪制的雙變量圖像
twoway 命令入門
twoway plot 變量 [if] [in] [, twoway_options]
plot：選擇圖像的種類，比如 scatter, line, connected, area, bar等
變量：這里可以寫一個或多個y變量，一個x變量。最后一個是x變量，之前的為y變量。
if：定義所取某一個自變量的范圍
in：定義所取觀測值的范圍
twoway_options：可以定義圖像的“美觀”部分，例如坐標軸范圍、標題、注釋、標簽等等
【上節課】
twoway scatter le year
twoway scatter le_male le_female year
le_male：y1? le_female：y2? year：x
折線圖 twoway line y1 y2 … x
twoway line le year ????????????????????twoway line le_male le_female year
帶數據標記的折線圖 twoway connected y1 y2 … x
twoway connected le year ??????????????twoway connected le_male le_female year
?
垂直線圖 twoway dropline y1 y2 … x
twoway dropline le year ????????????????twoway dropline le_male le_female year
?
脈沖圖 twoway spike y1 y2 … x
twoway spike le year ???????????????????twoway spike le_male le_female year 圖二
圖二 twoway spike le_male le_female year
為什么只顯示了female（灰色）？ twoway plot y1 y2 x，y2會覆蓋在y1之上
twoway spike le_female le_male year 處理方式，將順序調換
twoway spike le_female le_male year, scheme(s2color)
面積圖 twoway area y1 y2 … x
twoway area le year????????????????????? ?twoway area le_female le_male year
?
LOWESS圖（LOESS圖）twoway lowess y x
twoway lowess le year ??????????????????lowess le year
?
第八講 統計描述指標3
1.ci mean：連續變量mean（平均值）的置信區間
compute standard errors and confidence intervals
ci means [varlist] [if] [in] [weight] [, options]
cii means #obs#mean#sd [, level(#)]
默認95%置信區間（可通過 level(99)等進行更改）
ci mean mpg price, level(95)
cii mean 166 19509 4379, level(95)
直接告訴stata 觀測值的數量=166 # 平均值=19509 # 標準差=4379standard deviation
2.proportion（可簡寫為prop）：分類變量mean的置信區間
分類變量的置信區間
方法1 Choice 1：ci proportions [varlist] [if] [in] [weight] [, prop_options options]
ci prop foreign
缺點：ci prop只能用于binary二分類變量
codebook rep78
value 1 2 3 4 5 ·(缺失值)，有這么多種分類，所以不是二分類
ci prop rep78
方法2 Choice：proportion分類變量的置信區間
proportion varlist [if] [in] [weight] [, options]
默認95%置信區間
prop foreign
二分類變量，所以proportion 加起來=1
prop rep78
number of observations = 69，去除了缺失值后的數量
proportion1 2 3 4 5加起來=1
prop foreign rep78
prop foreign ??????????????????prop foreign rep78
觀測值數目=74 ???????????????觀測值數目=69
Domestic proportion =.6956522 ?Domestic proportion =.7027027
為什么會發生變化？
74-69=5，有5個缺失值
執行 prop foreign rep78【prop多個變量時，stata默認去除缺失值之后再計算proportion
【解決辦法】prop foreign rep78, miss
但是！！報錯了 沒有出現講課截圖中的樣子（下圖）
3.pwcorr：變量的配對相關性
pwcorr [varlist] [if] [in] [weight] [, pwcorr_options]
pwcorr price headroom mpg displacement
對角線上的數據都是1，因為price-price
price-headroom：0.1145 正相關 price-mpg：-0.4686負相關
pwcorr price headroom mpg displacement, sig【sig 顯著性 即p值，在每個相關性的下方顯示出來】
不想看所有的相關性的p值，只想看p值<0.05的。
pwcorr price headroom mpg displacement, star(0.05)
加*的為p值<0.05者
4.graph matrix：繪制相關性矩陣（圖像）
graph matrix varlist [if] [in] [weight] [, options]
graph matrix price headroom mpg displacement
|
可以觀察到 對角線兩側是完全對稱的，因此我們可以只看一半
graph matrix price headroom mpg displacement, half
第九講 可重復性：引入 reproducibility
Log File 所有顯示在results上的結果都保存為文檔形式
類似于R markdown文件
Do file 合作者、提醒自己、存檔
Do flie寫一個read me file 需要安裝fre
第十講 可重復性：Log文件
Log command
log using filename [, append replace [ text | smcl ] name( logname) ]
filename：Log file 的名字
任選其一 ?append：如果文件存在，附加append在文件上
replace：如果文件存在，替換replace該文件
如果文件不存在，append和replace都會創造新的文件
log using yikahui, append
smcl：stata默認的log file格式，可以保存各種顏色
text：單色，便于在文本編輯器中打開
作者偏向于text，比較干凈，體積較小
name(logname) 通常在一個stata運行的過程中，只有一個log file
已經有一個log file在運行了
有時，想創建不同log file給不同的合作者，又想在一個代碼中完成并記錄
因此：如果給log file起一個名字name( logname)，那么可以同時打開好幾個log file
不同的log file給不同的合作者
log using Brian, name( Log_For_Brian)
log using Allen, name( Log_For_Allen)
log using Jimmy, name( Log_For_Jimmy)
于是打開了3個不同的log file
關閉log file
沒有起名字（最常見）
log close
如果起了名字，需要指定名字
log close logname
log close Log_For_Brian
一口氣關閉所有log file，包括起了名字的
log close _all
第十一講 可重復性：Do file
打開文件 Lecture11.do
綠色文字為注釋，在stata中并不會運行，可以協商【該command是做什么用的、提醒自己注意的事項】
*只能作為一行的標記，/*? */ 可以跨越許多行，建議在每個do file的開頭做一個【Purpose、 Author、 Date created、 Last modified】，告訴自己該 Do file的用處
選中 clear all（清除掉stata內存中的所有數據），點擊do
打開 data browse，即發現是空的，代表stata內存里沒有數據
建議創建一個文件夾，里面保存所有數據log file和do file，便于數據的可重復性
pwd 查看現在的工作路徑
如果想改變工作路徑 cd “D:\stata17\外部命令“
定位好工作路徑之后，stata便會在該路徑下尋找操作所需的data，
log using "Lecture11", replace
use auto.dta, clear 調用一個 .dta格式的數據集
（該圖片可以看出 工作路徑 的意義）
如果數據集是其他格式的，如excel、csc等，之后的課程會講解如何進行數據導入
, clear 是一個比較保險的option，告訴stata：如果內存里已有數據，清除之后再導入
除了使用command，還可以使用窗口菜單，見下圖
通過窗口菜單使用的command，可以放入 do file中：全選 復制 粘貼 do
do file中，一行command可以只選中一部分，就運行do
可以不選擇任何一個command，直接運行該do file中全部command
第十二講 t檢驗
t檢驗包括 單樣本t檢驗 One-sample t test
獨立樣本t檢驗 Two-sample t test
配對樣本t檢驗 Paired t test
因為是盜版的，不能用 webuse
想看price 的 mean 是否等于6000，如下圖設置（by/if/in 設置各種條件判斷）
單樣本t檢驗 One-sample t test
ttest varname == # [if] [in] [, level (#)]
ttest price == 6000
level 默認95%水平
ttest price == 6000, level(99)
ttest price == 6000 if foreign == 0, level(90)
獨立樣本t檢驗 Two-sample t test
EG1：
EG2：
webuse fuel.dta, clear（因為是盜版 沒法使用聯網數據庫
在收藏夾里有stata官網，檢索相應的名字，即可找到.dta
>0.05不拒絕原假設
配對樣本t檢驗 Paired t test
雙樣本t檢驗
一個變量分組比較 ttest varname [if] [in], by(groupvar) [level(#)]
例子：ttest price, by(foreign)
兩個變量比較
非配對樣本：ttest varname1 == varname2 [if] [in], unpaired [level(#)]
例子：ttest mpg1 == mpg2, unpaired
配對樣本：ttest varname1 == vaarname2 [if] [in], [level(#)]
例子：ttest mpg1 == mpg2
第十三講 卡方檢驗/Fisher精確檢驗
sysuse nlsw88, clear
卡方檢驗 Chi-squared test
Pr<0.005 不同婚姻的race分布不同；不同race的婚姻狀況不同
frequency 487…
Chi-square contribution 16.7 9.3 …
expected frequency期望頻數（有時又叫理論頻數 586 1051…
在 key 可以看排列的內容種類
Fisher精確檢驗 Fisher’s exact test
通常，當列聯表中理論頻數（期望頻數）<5時，我們可以增加樣本量、刪去理論頻數太少的行或列、或者合并某些行或列；也可以使用Fisher精確檢驗
任何樣本量都可以使用Fisher精確檢驗
使用代碼進行卡方檢驗、Fisher精確檢驗
卡方檢驗：tabulate var1 var2, chi2
每個單元格對卡方檢驗的貢獻：tabulate var1 var2, cchi2 chi2 （cell chi-square==cchi cchi 和chi 順序可以互換）
理論頻數（期望頻數）：tabulate var1 var2, chi2 expected
Fisher精確檢驗：tabulate var1 var2, exact
注意：
1 總例數≥40，所有理論頻數≥5，看Pearson Chi-Square結果
2 總例數≥40，出現1個理論頻數≥1且<5，X2檢驗需進行連續性校正，以Continuity Correction結果為準
3 總例數≥40，至少2個理論頻數≥1且<5，看Fisher’s Exact Test結果
4 總例數<40或者出現理論頻數<1，看Fisher’s Exact Test結果
為什么可以不用連續性校正？
- stata不自帶卡方檢驗的連續性校正
- stata有用戶自寫的package，但不推薦；卡方檢驗的連續性校正并不是必須的，也不是最推薦的方法
- 在樣本量足夠大的時候，使用卡方檢驗時，是否使用卡方檢驗的連續性校正區別很小；使用Fisher精確檢驗也是沒有問題的
- 在樣本量小的時候（通常是某個格子期望頻數<5），可直接使用Fisher精確檢驗，亦不需要使用“卡方檢驗+連續性校正”
第十四講 RR值計算
同樣 去stata database 中檢索 csxmpl.dta 下載（如遇安全問題 選擇 “保留”）
sysuse csxmpl, clear
list
數據的結構
具有三個變量：case：1==case 0==non-case
exp：1==exposed 0==non-exposed 暴露/非暴露
pop 人數population
選擇pop人數進行加權，否則 stata四行每一行都會當作觀測值，但我們知道每一行是許多觀測值的合并
fweight 頻數加權
risk difference（相減）與0參照 置信區間 -0.7~-0.1 不包括0 所以p<0.05
risk ratio（相除）與1參照 置信區間0.2~0.9 不包括1 所以p<0.05
隊列研究計算器
使用代碼進行RR值的計算
cs var_case var_exposed [if] [in] [weight] [, cs_options]
cs case exp [fweight = pop]
csi #a #b #c #d [, csi_options] ?zhuyi注意 #不輸入
第十五講 OR值計算
OR值和RR值的區別？
RR值：Cohort study或者RCT中，研究者前瞻性地觀察“暴露組”和“非暴露組”的發病情況，之后通過RR來評價“暴露組”研究對象的發病風險是“非暴露組”研究對象的多少倍？
OR值：在回顧性研究（如case-control）中，研究對象是已經患病的“病例組”和未患病的“對照組”，研究者回顧性地調查病例組和對照組的暴露情況，因此無法計算發病率等指標。
最終目的：知道相對風險
OR是否可以估計RR值？
當終點事件發生率較低時，OR可以近似為RR（<15% 不是金標準）
當終點事件發生率較高時，OR會“夸大”RR“值
OR值相對于RR值“更遠離1“
當RR值大于1時，OR大于RR（1<RR<OR）
當RR值小于1時，OR小于RR（OR<RR<1）
終點事件發生率越高時，OR越會overestimate
對于隊列研究/RCT，可以報告OR么？
可以，但是：RR對于效應值的估計更加準確，對于臨床意義的解釋更加明確（病人理解，吃藥 比 不吃藥 使疾病風險降低50% √ 吃藥比值比是50% ×）
regression model回歸模型 中：對于結局是二分類變量的研究，logistic回歸只能提供OR值，不能提供RR值（之后會講：當結局發生率高時，應該使用log-binomial回歸或者使用帶有穩健方差估計的泊松回歸，直接提供RR值
接上節課RR內容
如果看risk ratio 該暴露降低了49%的風險，odds ratio 降低了87%，因此odds ratio是夸大了暴露因素的影響
強調：能提供RR值時（隊列研究、RCT）請不要提供OR值
回顧結束，本節課內容————————————————————
調入并觀察數據（匯總數據）
webuse ccxmpl, clear （依舊是手動 sysuse ccxmpl, clear）
對話框操作（使用數據集算OR）
?
病例對照優勢比計算器
使用代碼進行OR值的計算
cc var_case var_exposed [if] [in] [weight] [, cc_options]
cc case exp [fweight = pop]
cci #a #b #c #d [, cci_options]
第十六講 方差分析ANOVA
單因素方差分析 F值 與 單因素回歸分析 t值 是相同的 t2即F，p值也一樣
多因素方差分析 不如直接使用 多因素回歸分析 后者能知道每一個變量對結局的貢獻
回歸分析中 能更加方便的看兩個變量是否存在交互作用 ANCOVA（co-variance
下節課：線性回歸 logistic回歸 log-binomial回歸
方差分析的假設
假設1 y變量為連續變量
2 有一個包含2個及以上分類、且組別間相互獨立的x變量
3 每組間和組內的觀測值相互獨立（123在實驗設計階段盡力）
4 每組內沒有明顯異常值
5 每組內y變量符合正態分布
6 進行方差齊性檢驗，觀察每組的方差是否相等
sysuse systolic, clear
假設4：每組內有無明顯異常值？
方法：繪制 Boxplot（violin plot也可以
假設5：每組內y變量符號正態分布？
點擊 data browser看一下數據結構
systolic作為結局變量 drug分1234
codebook drug
sum systolic, detail
假設4 驗證
increment in systolic可能為負 第三個的下限能包括進該值，故暫不認為是異常值
假設5
- 偏度和峰度正態性檢驗
p值 0.1331 > 0.05 不能拒絕 符合正態分布 這一原假設
②Swilk-Shapiro-Wilk正態性檢驗
p值 0.37331 > 0.05 不能拒絕 符合正態分布 這一原假設
③sfrancia- Shapiro-Francia正態性檢驗
同上
三種檢測方式任選其一即可 很多人使用②，提前說好用哪一個，
完成假設4 5 的檢驗之后，可以開始oneway ANOVA
df 3 分子的自由度 54 分母的自由度 F 3 54分別是多少？肯定遠遠小于9.09 因為p 0.0001
所以四個組里至少有兩個組的平均值是有顯著差異的
假設6 方差齊性檢驗 觀察每組方差是否相等
Bartlett’s test for equal variances：p值0.800 >0.05 通過方差齊性檢驗，因此四個組方差齊
到底是哪兩個組呢？使用兩兩比較 示范用的Bonferroni
p值1.000 已經經過Bonferroni 處理 所以p值跟0.05比較
1-3 2-3 4-1 4-2 的p值<0.05 都是有顯著差別
row mean- col mean 組2-組1=-053 組3-組2=-17.31 等等
?
使用代碼
oneway systolic drug oneway y x 與之后的 regression 相似 regress y x? regress y x1 x2
oneway systolic drug, bonferroni
oneway systolic drug, bonferroni tabulate
twoway方差分析
?
第十七講 簡單線性回歸 simple linear regression
線性回歸的假設
假設1：y是連續變量
2：x可以被定義為連續變量（這一條不嚴格 x可以是啞變量
3：y和x之間存在線性關系
4：具有相互獨立的觀測值
5：不存在顯著的outlier
6：等方差性
7：residual近似正態分布
假設3：因變量和自變量之間存在線性關系
散點圖/Lowess圖（分別對應 第六講：散點圖 第七講：雙變量圖像2）
如果不符合線性關系怎么辦？
Spline：一次方向，分段fit直線
Quadratic：二次方項
Cubic：三次方項
Restricted cubic：三次方項（頭、尾近乎直線）
…
假設4：具有相互獨立的觀測值
可以使用杜賓-瓦特森（Durbin-Watson）統計量
stata對于非time series（時間序列）數據不設有這個統計量的檢測
更多情況下，不需要測量這個假設是否成立（實驗設計階段 即完成
如果是相互關聯的觀測值：GEE模型，Multi-level模型
假設5：不存在顯著的異常值
Boxplot或者violin plot（第五講
假設6：等方差性
使用Residual-versus-fitted plot（縱坐標residual橫坐標fitted
代碼：rvfplot
窗口菜單：Statistics → Postestimation
假設7：回歸殘差近似正態分布
step1：得到殘差
step2：使用正態分布的檢驗方法
直方圖
使用qq plot觀測
偏度峰度、Shapiro-Wilk檢驗、Shapiro-Francia檢驗
??
data browser
假設3：因變量和自變量之間存在線性關系
set scheme s1mono
twoway scatter mpg weight
基本是一個線性關系
twoway lowess mpg weight
lowess mpg weight 將lowess 和 散點圖 顯示在一起
?
?
假設4：具有相互獨立的觀測值（假定符合，因為都是不同的車嘛
假設5：不存在顯著的異常值
Boxplot或者violin plot（第五講
??最上面的橫線：75分位+IQR 并沒有多出去多少，所以姑且不算異常值
假設6：等方差性
?
regress mpg weight 建議回歸分析直接敲代碼 因為簡單：regress y x / regress y x1 x2…
F（1，72）整個模型有無共線：比如有好幾個x，測所有變量的β值都等于0，（？？）
R-squared納入模型的變量能夠解釋y變量的多少 65% 64%
當納入的變量越多，R-squared↑，但大多時候這對模型沒有意義，因為可能有共線性、過擬合等等；Adjusted R-squared 進行過校正更好顯示貢獻
截距項39.44（weight=0時 mpg=39.44）
rvfplot
fitted value 即 predicted y value
看residual是否分布在0兩側
???
假設7：回歸殘差近似正態分布
predict cancha, residual
cancha 隨便取一個新變量的名字（cancha殘差
data browser
做一個cancha殘差的直方圖
??jinsi
近似正態分布
也可以用qq plot觀測
?
如果貼合這條直線就是完全正態分布，
使用偏度峰度，Shapiro-Wilk檢驗、Shapiro-Francia檢驗
?
p值<0.05
第十八講 多元線性回歸
sysuse auto.dta, clear
price=β0+β1weight? regress price weight
β1：汽車的重量每增加一個單位，售價增加2.04個單位（95%CI：1.29，2.80）
p=0.000？×? p<0.001 √
β0怎么解釋？：當汽車的重量=0時 汽車的售價=-6.70，但這不合理：汽車重量不可能=0，售價也不可能為負數
price=β0+β1weight+β2length regress price weight length
β1：在控制了汽車的長度后，汽車的重量每增加一個單位，售價增加4.70個單位（95%CI：2.46，6.94）
β2：…重量…，長度…，售價降低97.96個單位（95%CI：-176.08，-19.85）
price=β0+β1weight+β2length+β3mpg regress price weight length mpg
汽車的價格是否和修理次數有關？
codebook rep78
5個缺失值
price=β0+β1weight+β2length+β3mpg+β4rep78 regress price weight length mpg rep78
β4：在控制了汽車的重量、長度、里程后，汽車的修理次數每增加一個單位，售價增加910.99個單位（95%CI：302.62，1519.35）
這里rep78當作了連續變量，但應該當作多分類變量--
Dummy variable：啞變量
當自變量x為多分類變量時，例如職業、學歷、血型、疾病嚴重程度等等，此時僅用一個回歸系數來解釋多分類變量之間的變化關系，及其對因變量y的影響，就顯得不太理想
在多分類變量前加上“i.“，告訴Stata這不是連續變量
regress price weight length mpg i.rep78
rep78=1 被當作對照組，
Irep78_2：在控制了汽車的重量、長度、里程后，汽車的修理次數兩次和一次相比，售價增加1137.28個單位（95%CI：-2468.70，4747.27）
回歸分析中需要注意的事情！
1、時刻注意納入模型的觀測值數量
?rep78中有5個缺失值，
Case-wise Deletion 只用納入模型的變量都沒有缺失值，這個觀測值才會納入回歸模型分析！
也就是 該observation 在 price/ weight/ length/ mpg/ rep78五個變量中都沒有缺失值，才能納入
2、如何讓β0有意義？
汽車平均重量：3019.459（sum weight）
轉化為 price=β0+β1（weight-3019.459）
β0：當weight=平均重量時，汽車的售價
gen weight_center=weight-3019.459
regress price weight_center
窗口菜單操作
?
第十九講 二分類Logistic回歸 Logistic regression(binary outcome)
回顧第十五講 OR是否可以估計RR值？“當終點事件發生率較高時，OR會夸大RR值”
橫坐標 終點事件發生率 縱坐標 OR值
可以看到，終點事件發生率20+時，RR值為3（真實的relative risk），OR值為8，夸大了
（偏離RR值越大 不論是減少 或增加 都稱之為overestimate，因為相關性 也有正負 值的絕對值越大 相關性越強）
當outcome發生率>15%時，logistic regression得出的OR值會overestimate實際的RR值
解決辦法
做Logistic回歸之前檢查
條件1：Y是二分類分類變量
條件2：Y的發生率<15%
（條件1、2就是區分 Logistic 和 log binominal 的區別
使用二分類Logistic模型前，需判斷是否滿足以下7項假設
假設1：因變量（結局) 是二分類變量
2：有至少1個自變量，自變量可以是連續變量，也可以是分類變量。
3：每條觀測間相互獨立。分類變量（包括因變量和自變量) 的分類必須全面 且 每一個分類間互斥。
4: 最小樣本量要求為自變量數目的15倍，但一些研究者認為樣本量應達到自變量數目的50倍。
5：連續的自變量與因變量的logit轉換值之間存在線性關系。
6：自變量之間無多重共線性。
7：沒有明顯的離群點、杠桿點和強影響點。
sysuse lbw, clear
變量low（<2500g）是結局事件，想了解什么因素和孩子的low birthweight 相關
codebook low
tab low
發生率31.22%>15%
disclaimer免責聲明：本次課的數據集中，結局事件的發生率遠大于15%，應使用Log-binomial模型進行分析。本次課使用Logistic regression進行分析僅為了講解如何使用stata進行操作，以及為下節課講解Log-binomial進行鋪墊
Logistic regression
Stata code：logistic y x1 x2 x3…
窗口菜單：
Model1：low=β0+β1age（母親age與孩子low weight是否有關
logistic low age
β1：母親年齡每增加1歲，孩子低體重風險是之前的0.95倍（95%CI：0.89，1.01）
Model2：low=β0+β1age+β2smoke
logistic low age i.smoke
β1：控制了母親吸煙狀況后，母親年齡每增加1歲，孩子低體重風險是之前的0.95倍（95%CI：0.89，1.01）
Model3：low=β0+β1age+β2smoke+β3race
logistic low age i.smoke i.race
β1：控制了母親吸煙狀況和種族后，母親年齡每增加1歲，孩子低體重風險是之前的0.97倍（95%CI：0.90，1.03）
β2：控制age和race后，吸煙者與不吸煙者相比，孩子低體重風險是3.01倍（95%CI：1.45，6.23）
β3：控制age和smoke后，black與white相比，孩子低體重風險是2.74倍（95%CI：1.04，7.23）
【stata會默認將值最小的作為多分類變量中的reference，此處即race中的white】
窗口操作
第二十講 廣義線性模型 generalized linear model, GLM
回顧：大于15%會產生bias，這個問題只出現在我們想用odds去estimate risk的時候才會發生
在一個研究中，我們只想報告odds，而不去estimate risk，無論結局事件發生率是多少，logistic回歸都是valid。
線性回歸（lecture17-18）：Expected Y=β0+β1*X1+…+βn*Xn（不用對Y進行轉換
廣義線性模型：f（expected Y）=β0+β1*X1+…+βn*Xn（f 意味著function）
simple GLM 與 Multiple GLM
簡單廣義線性模型 f(expected Y)=β0+β1*X1
多元廣義線性模型 f(expected Y)= β0+β1*X1+…+βn*Xn
x變量可以是各種類型的變量（連續、分類）
in GLM，we have to specify 2 things
family ?== Y變量如何分布？ 比如正態分布、泊松分布、二項分布
link ?==把Y怎么轉化才能和X成linear關系？
Linear regression 最簡單的線性回歸
Y variable? ①continuous variable
②如何分布？ 高斯分布（Gaussian Distribution）
它的family是：我們有一個assumption, 在線性回歸中Y變量
應該是正態分布或高斯分布
model ?①Expected Y=β0+β1*X1+…+βn*Xn
②Y怎么轉化才和X成linear關系？不用轉化！link: Identity(恒等 意思：不轉化它便相等，即為恒等（？)
Logistic regression
Y variable? ①Binary variable
②如何分布？ 二項分布（Binomial）
model? ①ln [ P ( Y=1 )/ P ( Y=0 ) ]=β0+β1*X1+…+βn*Xn? （ P(Y=0)可以寫作 1-P( Y=1 ) 因為是二項分布）
②Y怎么轉化才和X成linear關系？link: Logit
Poisson regression
Y variable? ①count variable：整數（1, 2, 3…n）
②如何分布？泊松分布（Poisson）
model? ①ln [risk of event]= β0+β1*X1+…+βn*Xn
②Y怎么轉化才和X成linear關系？link: Log
Log-binomial regression
Y variable? ①binary variable
②如何分布？二項分布（Binomial）
model? ①ln [ P ( Y=1) ]= β0+β1*X1+…+βn*Xn
②Y怎么轉化才和X成linear關系？link: Log
help glm 可以查看該表內容
Linear regression: lecture 17-18? regress y x1 x2 x3
Logistic regression: lecture 19? logistic y x1 x2 x3
Poisson regression 右側二維碼 ?poisson y x1 x2 x3
glm command in stata 就像一把雨傘 什么都可以裝在里面
umbrella command
linear regression: family ( gaussian ) link( identity)
logistic regression: family ( binomial ) link( logit)
poisson regression: family ( poisson ) link( log )
log-binomial regression: family ( binomial) link( log )
Linear regression
regress price weight length mpg i.rep78
glm price weight length mpg i.rep78, family(gaussian) link(identity)
Logistic regression
logistic low age i.smoke i.race
glm low age i.smoke i.race, family( binomial) link( logit)? 【logistic regression 報告odds ratio，e^n 報告的是n】
glm low age i.smoke i.race, family( binomial) link( logit) eform 【但我們想知道odds ratio 即e^n ，本行command報告的是e^n】
對比
regress price weight length mpg i.rep78 (見lecture 18)
glm price weight length mpg i.rep78, family(gaussian) link(identity)
logistic low age i.smoke i.race (見 lecture 19)
glm low age i.smoke i.race, family( binomial) link( logit) eform
Log-binomial Model
glm low age i.smoke i.race, family( binomial) link( log)
alternative to Logistic regression【是logistic regression 的替代
問題：fail to converge (對于一個command 一直跑代碼 停不下來
解決辦法solution：Poisson regression with robust variance estimate 使用帶有穩健方差估計的泊松回歸
glm low age i.smoke i.race, family( poisson) link( log) robust ×
glm low age i.smoke i.race,family(poisson) link(log) robust eform √
【此處老師講的代碼和圖不匹配 應該是eform這一行 對應視頻中的IRR的圖。感謝評論區
【incidence rate ratio 即IRR，此處當成relative risk也可以
窗口菜單
?
第二十一講 Kaplan-Meier曲線和Log-Rank檢驗
回顧：Survival analysis in use 做生存分析時的目標：
①描述一個組內個體的生存時間：壽命表法（life tables methods）、非參數Kaplan-Meier曲線（lecture 21）（后者逐漸成為主流）
②比較兩個或多個組的生存時間：Log-rank test（lecture 21）
③研究生存時間和變量之間的關系：半參數Cox比例風險模型（lecture 22）、參數生存分析模型（too advanced, not covered in these 30 lectures）
K-M曲線的歷史：
介紹了一種全新的、解決隨訪期間Right Censoring的問題的生存分析方法
特點：精確地記錄并利用每個個體發生終點事件的具體時間，在任何一個終點事件發生的時間點計算出一個新的、基于之前所有信息的Cumulative survival
優點：①相比于life-table method，更加充分地利用了信息，給出更準確的統計量
②非參數估計方法：不要求總體的分布形式，因此非常適合（初期的）生存分析時使用
③K-M曲線可以很直觀的表現出兩組或多組的生存率或死亡率，適合在文章中展示
手動下載drugtr.dta
導入數據：sysuse drugtr, clear
恢復成為普通的數據格式：stset, clear
stata已經將這個數據集設置成了生存數據的格式，導入數據集后，將數據集恢復成普通的數據格式，才是我們在臨床研究中見到的數據結構
step1 數據集的初步觀測
list values of variables
list
list in 5/10 呈現5 6 7 8 9 10這六個觀測值
有四個變量，每個的觀測值可以看見
codebook drug
0=placebo 安慰劑 所以1=試驗藥
codebook studytime
step2 聲明數據為生存分析數據代碼操作
告訴stata：終點事件（failure variable），隨訪時間（time variable）
stset timevar, failure ( failvar [==numlist ] )
timevar：隨訪時間變量
failvar：終點事件變量
numlist：終點事件變量中，哪些值代表發生了終點事件（例如1代表終點事件，numlist即為1）
stset studytime, failure ( died==1)
窗口菜單
?
多記錄的ID變量multiple-record ID variance：
step3 生存數據再觀測
注：必須要在指定數據集為生存分析數據集之后（stset之后）才能使用任何其他的st開始的命令。
stsum
50%中位生存時間，很在乎
窗口菜單
stdescribe
窗口操作
sts list
stata會在每一個時間點上計算出一個新的tumor list survival
窗口操作
step4 K-M曲線的繪制
可以畫生存曲線 或者 死亡曲線
代碼：sts graph, by( var )
?
sts graph, by( drug )
sts graph if age<50, by(drug) 繪制age<50的人
sts graph, by(drug) risktable
?
期待：高級繪圖1-2 lecture29、30 講解
fancy 的圖捏！
step4 檢驗組間差別（Log-Rank Test）
sts test drug 比較drug=1 和 drug=0這兩個組的survivor functions是否一致
p<0.0001 因此有顯著差別，得出結論：實驗藥能夠顯著提高生存率
第二十二講 Cox回歸和比例風險假定檢驗
回顧：Kaplan-Meier曲線： 描述一個組內個體的生存時間
對于兩組或多組患者的生存率、死亡率進行直觀的比較
Log-rank檢驗：比較兩個或多個組的生存時間
得到統計學上的相應證據，即p值
但：Kaplan-Meier曲線是一種非參數的估計方法，不能控制潛在的混雜因素對于結局事件的影響，也不能看出某種因素，如患者的年齡、性別、種族等，是否是結局發生的獨立因素或危險因素等，
question：在一個抗癌藥物的Clinical Trial中，48名患者被隨機分配到新藥組（28人）和安慰劑組（20人），研究人員想知道新藥是否影響患者的生存情況
sysuse drugtr, clear
提示：數據初步觀測 和 轉化為生存數據格式 在上節課已講，在此不贅述
step1 Cox回歸
窗口菜單
statistics → survival analysis → regression models → Cox proportional hazards model
我們將數據轉化為survival data時（stset），指定過終點事件（failure variable）、時間變量（time variable），我們在此只需告訴stata 放入哪些x變量
drug的coefficient 新藥組終點事件發生風險是安慰劑組的13.3%（95%CI：5.6%，31.4%）
“與安慰劑組相比，新藥組發生終點事件的風險降低了86.7%”
drug的coefficient 在控制了患者年齡age后，新藥組終點事件發生風險是安慰劑組的10.5%（95%CI：4.3%，25.6%）
age的coefficient 在控制了治療方法drug后，患者年齡每增加1歲，發生終點事件的風險增加12.0%（95%CI：4.1%，20.5%）
代碼操作
stcox va1 var2 var3 … [if] [in] [, options]
再提醒：①必須指定data為survival data后（stset）才能使用任何st開頭的命令
②因一開始已將data轉化為survival data時，指定過終點事件failure variable、時間變量time variable，我們在此只需設置回歸方程中independent variable即可
example：stcox drug age
stcox drug age if age<50
step2 PH假定檢驗
窗口菜單
statics → survival analysis → regression models → test proportional-hazards assumption
PH是一個必須要緊跟在回歸模型之后的，一個對于模型的評價和估計
其零假設：納入Cox回歸模型的變量滿足PH假定
p=0.8064>0.05，不能拒絕零假設，所以滿足PH
注：若前面沒有進行Cox回歸，輸入代碼 estat phtest，會報錯
通過圖像觀測變量是否滿足PH假定
statics → survival analysis → regression models → graphically assess proportional-hazards assumption
??set scheme s1mono
如果兩條線近乎平行，則滿足PH假定
?
調整了age之后，變化不大，依舊近乎平行
代碼操作
eatat phtest
再強調：該命令需緊跟在Cox回歸的命令之后，否則stata不知道檢驗哪個回歸的PH假定
stphplot, by ( var1 ) adjust ( var3 var3… )
var1是自變量名，var2是希望控制的變量。注：這個命令不一定要跟在cox回歸之后
stphplot, by ( drug ) adjust ( age )
如果不滿足PH假定怎么辦？
比如本例之中的drug）的交互項（interaction term），即time*exposure
【生成一個新變量放入模型，觀察藥物在不同時間的作用效果放松PH假定
第二十三講 廣義估計方程 generalized estimating equation GEE
examples of longitudinal data analysis LDA，縱向數據分析
child BP measured at each annual visit from 3 to 9 years old
infant weight measured at 3,6,9,12
in a 3-arm cross-over trial, short chain fatty acid measured at the end of each other
(cross-over trial 理解：實驗中三個組 高糖高脂高蛋白，先把人進組，吃一段時間高脂，洗脫，吃一段時間高糖，洗脫吃一段時間高蛋白)
縱向數據分析到底比橫斷面分析強在哪兒？
a classical example: Association between age and reading ability? 孩子的年齡與閱讀理解能力
若是橫斷面研究，看起來似乎年齡↑閱讀能力↓
不論基線水平，對一個孩子而言，年齡↑閱讀能力↑
更有說服力，雖然和橫斷面一致
advantage of LDA 優勢
establish temporal trends (時序：A在B之前發生)
separate cohort effects from aging effects
children have different baseline values in reading abilities
the trajectories of reading abilities differs by people
from mathematical perspective
若不視作一個個體的mpg1 mpg2，而分開
?
p值變小了，差值diff沒變，故是standard error變小了（1.23→0.78），
Var ( X-Y ) = Var ( 1 x X +(-1) Y )=( 12 ) Var ( X ) + ( -12 ) Var ( Y ) + 2 ( 1 ) ( -1 ) Cov (X, Y)
非配對t檢驗中，紅色=0，配對t檢驗中，紅色≠0，因此Var↑
LDA中，when estimating the “change”, taking into account the correlation between observations will give us small SE ( standard error )
線性回歸的假設中，觀測值需要有相互獨立（lecture17 假設4）
within an individual, the observations are correlated
the next question is: how are they correlated?（觀測值如何關聯）
independent correlation structure 最常用
如果實際情況不符合independent correlation structure，Robust variance estimate if worried about mis-specification of the correlation structure
use QIC(an alternative for AIC for GEE models) to see which model works best
data: NLS women 14-24 in 1948 (national longitudinal surveys)
sysuse union, clear
觀察數據list in 1/10
ID=1的人在72年隨訪，在71年沒有，id=2在71年隨訪
告訴stata：I want to use GEE!
xtset studyid timevar
studyid: unique ID for each participant
timevar: timevar will usually be a variable that counts 1, 2, …, and is to be interpreted as first year of survey, second year, …, or? first month of treatment, second month.
xtset id year
GEE command
xtgee y x1 x2 x3 … [if] [in] [weight], family (family) link (link) corr(correlation structure) robust
family: lecture 20 ; or help xtgee
link: lecture 20 ; or help xtgee
correlation structure
xtgee union age grade not_smsa south, family(binomial)link (logit) corr(ind) robust
第二十四講 有序多分類Logistic回歸 logistic regression(ordinal)
多分類變量分為有序多分類和無序多分類變量
有序多分類變量：疾病分期（1-4期）；疼痛情況（1-10：1最輕，10最痛）
無序多分類變量：居住地（1東/2南/3西/4北）；手機品牌（1蘋果/2三星/3華為/4小米/5其他）
有序多分類的原理
將y變量的n個分類拆成n-1個二分類Logistic回歸
今天例子中的y變量有5個level：Excellent; good; average; fair; poor拆分成：
poor(1) vs excellent + good + average + fair(0)
fair + poor(1) vs excellent + good + average(0)
average + fair + poor(1) vs excellent + good(0)
good + average + fair + poor(1) vs excellent(0)
proportional odds假定
odds ( poor) /odds ( excellent + good + average + fair )? OR1
odds ( fair + poor) /odds ( excellent + good + average )? OR2
odds ( average + fair + poor ) /odds ( excellent + good)? OR3
odds ( good + average + fair + poor) /odds ( excellent )? OR4
本節課使用數據 1977年汽車修理記錄數據
sysuse fullauto, clear
codebook rep77 outcome：汽車維修狀況
?結局變量5個level
codebook foreign exposure：是否為進口車
step1 Chi-squared test(卡方檢驗)
H0：車輛是否為進口車和車輛維修狀況沒有關系
tab foreign rep77, chi2
p=0.008<0.05 拒絕零假設，得出結論：車輛是否為進口車和車輛維修狀況有關系
→下一步：如何量化這種關系？
step2 Ordinal Logistic Regression
ologit y x1 x2 x3 [if] [in] [weight] [, options]
常用的options包括or（若不加or，stata給出結果為coefficient，得不到OR
examples ?ologit rep77 foreign
ologit rep77 foreign, or
ologit rep77 foreign length mpg, or
ologit rep77 foreign
1.46 進口車（foreign=1）和國產車（foreign=0）比：
odds (poor) / odds (excellent + good + average + fair)
= odds (fair + poor) / odds ( excellent + good + average)
= odds (average + fair + poor) / odds ( excellent + good)
= odds (good + average + fair + poor) / odds ( excellent)
=e^(-1.46)=0.23 （更高維修狀況等級為reference，在更低維修狀況的odds為0.23）
stata實際進行的分析操作 ?不夠直觀
此處stata已經幫我們轉化成 比較差的維修狀況=0，為reference
=e^(1.46)=4.29（更低維修狀況等級為reference，在更高維修狀況的odds為4.29）
進口車和國產車相比，在“更低的車輛維修狀況等級”的odds 是在“更高維修狀況等級“的0.23倍
進口車和國產車相比，在“更高的車輛維修狀況等級”的odds 是在“更低維修狀況等級“的4.29倍
ologit rep77 foreign, or
一定要記住stata給出的OR是以較低的為reference
ologit rep77 foreign length mpg, or
在控制了汽車的長度、里程之后，進口車有著更高車輛維修狀況等級的odds是國產車的18.12倍（95%CI：3.85，85.32）
在控制了汽車的產地、里程之后，車輛每增加1 inch，有更高車輛維修狀況的odds增加8.64倍（95%CI：3.90，13.58）
mpg？如何解釋
窗口菜單
statistics → ordinal outcomes → ordered logistic → regression
?
第二十五講 無序多分類Logistic回歸 logistic regression(multinomial)
sysuse fullauto, clear
ologit rep77 foreign, or
進口車（foreign=1）有著更高車輛維修狀況等級 的odds是國產車（foreign=0）的4.29倍（95%CI：1.51，12.13）
安裝gologit2命令
gologit2命令
gologit2 y x1 x2 x3…, pl or(pl 又叫parallel line 假定)
這個command 和 ologit command 給出的結果相同
gologit2 y x1 x2 x3…, npl or
likelihood-ratio test: lrtest(檢驗p值是否小于0.05)
gologit2 rep77 foreign, pl or
原理：5個level分成4個二分類變量
category>poor vs cat≤poor
cat>fair vs cat ≤fair
cat>average vs cat ≤average
cat>good ??vs cat ≤good
有序多分類 當proportional odds假定成立時
進口車（foreign=1）和國產車（foreign=0）比：
odds (excellent + good + average + fair) / odds (poor)
= odds ( excellent + good + average) / odds (fair + poor)
= odds ( excellent + good) / odds (average + fair + poor)
= odds ( excellent) / odds (good + average + fair + poor) =4.29
gologit2 rep77 foreign, npl or
3.94*10^7這個estimate是很不準，上節課 有一個level-excellent的車輛數量是0的，所以estimate很大，standard error也很大，但并不影響interpret results
同樣是 category>poor vs cat≤poor
cat>fair vs cat ≤fair
cat>average vs cat ≤average
cat>good?? vs cat ≤good
有序多分類變量 當proportional odds假定不成立時
進口車（foreign=1）和國產車（foreign=0）比：
odds (excellent + good + average + fair) / odds (poor)=0.93
odds ( excellent + good + average) / odds (fair + poor)=3.45
odds ( excellent + good) / odds (average + fair + poor)=3.28
odds ( excellent) / odds (good + average + fair + poor) =3.94*10^7
檢查proportional odds假定是否成立
H0：non-proportional odds 模型可以更好解釋結局變量各個等級之間關系
gologit2 rep77 foreign, pl or
est store A (將上面的model存到A中
gologit2 rep77 foreign, npl or
est store B
lrtest A B
p>0.05，拒絕H0，non-proportional odds并沒有更好解釋結局變量各個等級之間關系，因此可以說滿足proportional odds假定
以上為有序多分類變量不滿足proportional odds假定的情況，接下來為無序多分類logistic回歸
無序多分類logistic回歸
把結局變量的某個分類作為reference，然后比較結局變量其他分類相對于reference的相對風險（relative risk）
?j指的是結局分類其他變量
有序和無序多分類比較
有序多分類logistic回歸：ORj =Pr (cat > j) / Pr (cat ≤ j)
ologit y x1 x2 x3…, or
無序多分類logistic回歸：RRj =Pr (cat = j) / Pr (reference cat)
mlogit y x1 x2 x3…, rrr baseoutcome ( j ) ( baseoutcome 不是必須的，若不指定，stata會自動給一個reference level進行比較
mlogit rep77 foreign, rrr baseoutcome (1)? poor=1即baseoutcome
無序多分類
進口車（foreign=1）和國產車（foreign=0）比：
Risk (fair) / Risk (poor) =0.20
Risk (average) / Risk (poor) =0.70
Risk (good) / Risk (poor) =1.08
Risk (excellent) / Risk (poor) =1.32*10?7
窗口菜單
statistics > categorical outcomes > multinomial logistic regression
?
注意gologit2非stata自帶command 因此沒有窗口菜單操作
26-28講：數據清理（包括數據導入
29-30講：高級繪圖
想學的沒有學到？
第二十六講 數據整理1 intro to data management 1
26-28講，我們將使用Do-File（第十講 第十一講）
Stata的須知
Stata區分大小寫！Stata is case-sensitive（無論是變量的名字 還是command
常見符號
= 是 賦值（e.g. gen x=1
== 是 恒等（e.g. gen x=1 if y==2
| 是 或者（e.g. gen x=1 if y==2 | y==3
& 是 且（e.g. gen x=1 if y==2 & z==3
數據整理常見步驟（stata command）
查看工作路徑（pwd）；改變工作路徑（cd）
導入excel文件（import excel），CSV文件（import delimited），dta文件（use）
添加變量或變量數值標簽（label）
生成新變量（gen）或統計量變量（egen）
講觀測值按照變量數值大小排序（sort；gsort）
改變變量前后順序（order）
將數據集進行長寬轉換（reshape）
合并數據集（merge）
刪除（drop）或保留（keep）觀測值或變量
導出excel，CSV文件（export）或dta文件（save）
復習stata中添加注釋的方法：stata中綠色的文字都是不運行的
三種添加注釋的方法，見do file 26
do clear all? 注意 // 與all之間需要有空格，否則會報錯
Stata IC最多只有2,048個變量; SE:32,767; MP: 120,000
盡管是MP，但一開始打開stata時，默認5000給變量，因此：如果你的變量數>5000，需要設置。例如: set maxvar 8000
正式數據分析開始之前，我們可以打開一個log文件，這樣可以將屏幕上的東西都保留在log文件里，便于與co-author其他作者進行分享
若沒有指定，自動保存在工作路徑下
use "child.dta", clear 該command沒有指定去哪兒找child.dta 因此會去pwd工作路徑下尋找
use child.dta, clear 雙引號可以不加，但如果文件名中有空格，就必須加“”雙引號，因此建議都加雙引號“”
import delimited "child.csv", clear 報錯了！！！ 以前用的方法 在Index of /RePEc/bocode/i (bc.edu) 下載import command 也沒能解決 因此有了下面的鏈接的方法
stata17運行遇到“Java installation not found” - 知乎 (zhihu.com)
安裝好stata17后遇到Java installation not found問題解決方案，親測可行
1、下載java17以上的版本：https://www.oracle.com/java/technologies/javase/jdk17-archive-downloads.html
（注意：下載壓縮包版本，windows64位的不建議下載jbk.17.0.6版本，因為很慢）
2、將壓縮包解壓后將文件放到stata17安裝的路徑下
3、在stata中輸入命令（路徑地址根據自己java保存的地址更改）
set java_home "C:\Users\ramfa\Desktop\jdk-17_windows-x64_bin\jdk-17.0.2\"
4、指令運行結束后，關閉stata，再打開就好了
方法來源于簡書作者陸陸柒柒https://www.jianshu.com/p/95fe650fb
上述文件已下載在：D:\stata17\外部命令\jdk-17.0.6_windows-x64_bin
import delimited "child.csv", clear
此處stata很聰明的把變量名導入了，沒有把它當作第一行的observation，可以看到變量名都是小寫的
在導入CSV文件時，Stata自己判斷是否把CSV文件的第一行作為變量名導入（但是有時會出錯 因此更保險的方法如下：指定stata不把csv的第一行作為變量名導入 import delimited "child.csv", varnames(nonames) clear
我們也可以讓Stata必須把CSV文件的第一行作為變量名導入 import delimited "child.csv", varnames(1) clear
導入delimited file時，stata會自動把變量名都變成小寫；我們可以讓Stata導入的時候把變量名全變成大寫 import delimited "child.csv", varnames(1) case(upper) clear
其他內容見Lecture 26.do 文件
以上代碼很多 記不住怎么辦？可以窗口菜單操作
??
窗口菜單 運行的命令，如果想要保存至log中，下次直接用：import delimited "D:\stata17\外部命令\數據導入-stata教程-醫咖會\child.csv", delimiter(comma) varnames(1) case(upper) clear ?
log.file中一行太長不好看，換行的話，可以用 ///，記得///和command之間需要空格
在導入Excel文件時，Stata默認第一行不是變量名…余下見lecture26.do文件
import excel "child.xlsx", clear
import excel "child.xlsx", firstrow clear
也可以使用窗口菜單導入excel文件
?
其他形式數據怎么導入？窗口菜單（本節課學了dta csv excel）
還有其他不支持的格式 可以用 stat/transfer? 或者 R語言進行數據轉換，一般轉化成csv文件是最穩妥的，因為不管什么統計學軟件都可以讀取csv
第二十七講 數據整理2 intro to data management 2
依舊使用lecture26
import excel "child.xlsx", firstrow case(preserve) clear
以CHILD_SEX變量為例，進行數據的清洗和整理
codebook CHILD_SEX 字符型變量 因為“1” “2” “M”有雙引號
option：replace? 在destring時 replace原有變量或generate新變量
給變量和變量的數值加標簽
注意：都有空格
生成新變量 ?變成二分類變量
* 特別注意：在Stata中，缺失被定為為無限大
觀測值排序
sort 從小到大 gsort + 從小到大 gsort -從大到小
第二十八講 數據整理3 intro to data management 3
依舊用lecture26.do
打開browser 可以直接點擊 也可以browse
變量排序
合并數據集
child.dta 每個孩子一條觀測值，也就是每個孩子只有一行，現在我們想合并anthro.dta，為一個孩子的各種數據的dta，每個孩子有好幾條觀測值，好幾行數據，分別是不同的時間點進行的測量，是一個長的long dataset
打開一個數據集，打開的就叫master dataset
master dataset（child）?? using dataset（anthro）
one-to-one???? 一行 ?????????????????一行
many-to-one??? 多行 ?????????????????一行
one-to-many??? 一行 ?????????????????多行（本次就是第三種情況）
many-to-many? 多行（未知） ??????????多行（未知）（自己不知道多少行，讓stata自己決定，耗時長，并且不佳，知道的話盡量還是對應用以上三種之一）
3580未merge，并且來自master，說明anthro是child的子集，而child是初始baseline，因此有孩子失訪了
自動生成了一個新變量_merge，master only即只在master中，不在using dataset中
只想保留match成功的（既有maternal information又有infant anthropology measurement的孩子）：keep if _merge == 3 ??drop _merge? 完成數據集合并
reshape 有時我們有一個長數據集，想變成寬數據集，反之亦然
help reshape 理解：i是id，j是時間點，stub為測量內容（如身高），長數據集：時間點是分開的，寬數據集：時間點跟在stub后面
child數據集中 j對應visit（各個時間點）
reshape long stub, i ( i ) j ( j )? （reshape (to) long變成長數據集）
reshape wide CHILD_ADJ_AGE WEIGHT-HEAD_CIRC, i(CHILD_PIDX) j(VISIT)
-意為從x1變量到xn變量
long → wide 變量名增加了，visit被drop掉了
8096為2024的四倍
因為：每個孩子的visit不一定都有四行，long→ wide→long，還原的時候stata自己補充了
egen（egenerate）生成統計量變量
by CHILD_PIDX: egen max_visit = max(VISIT)
egen height_mean = rowmean(HEIGHT6 HEIGHT12 HEIGHT24 HEIGHT36)? 四個時間點的身高平均值
egen height_miss = rowmiss(HEIGHT6 HEIGHT12 HEIGHT24 HEIGHT36) 看這四個變量中，每個孩子到底有幾個缺失值
有沒有確實少一點的？tab height_miss
egen height_nonmiss = rownonmiss(HEIGHT6 HEIGHT12 HEIGHT24 HEIGHT36) 看一下沒缺失的有幾個 height_miss 和 height_nonmiss 互補
help egen 自己探索！很強大的功能 ?多看看下面的examples
export excel "Sample", replace //有問題！
沒有保存變量名字（第一行）！
export excel "Sample.xlsx", firstrow(variables) replace //“變量名”作為第一行導出
excel打開csv文件，亂碼，但其實數據是完好的，因此最好把變量和變量標簽都保存成數值格式，界面也盡量使用英文的
窗口菜單操作
可以勾選是否保留第一行名字等
?
第二十九講 高級繪圖1 advanced graphing 1
第8行 保留變量keep 保留觀測值keep if ；保留了所有matched觀測值child
9… 去掉 變量_merge 多余的行，因為這一行的作用就是篩選下_merge==3的child
13…
14…
15…不知道藍/紅色是誰，待會教加標簽
22…
換不起行？因為///前沒有空格，stata把///當成command了
√
缺點：都排在一列上了，想讓他們更分散均勻點
option：jitter ( ) 抖點-讓點散開
twoway (scatter WEIGHT MOM_AGE if CHILD_SEX==1, jitter(2)) ///
(scatter WEIGHT MOM_AGE if CHILD_SEX==2, jitter(2))
缺點：單個點太大了，顏色可以接受
?twoway (scatter WEIGHT MOM_AGE if CHILD_SEX==1, jitter(2) msize(tiny)) ///
(scatter WEIGHT MOM_AGE if CHILD_SEX==2, jitter(2)? msize(tiny))
twoway (scatter WEIGHT MOM_AGE if CHILD_SEX==1, jitter(2) msize(tiny)) ///
(scatter WEIGHT MOM_AGE if CHILD_SEX==2, jitter(2)? msize(tiny)) ///
(lfit WEIGHT MOM_AGE if CHILD_SEX==1) ///
(lfit WEIGHT MOM_AGE if CHILD_SEX==2)
缺點：顏色有黃有綠，更不知道誰是誰了
缺點：置信區間有點細，想加粗 clwidth(thick)
觀察到x軸在15之后，45之前出現點，想改一改x軸
?xlabel(15(5)45)
y軸想延長到25 ylabel(5(5)25)
?xtick(15(1)45) ytick(5(1)25) 間隔1 畫tick，就是坐標軸的小尺度
做好一個圖后，更改細節就很快很方便了
第三十講 高級繪圖2 advanced graphing 2
學習一個新的command：coefplot
如果之后的文章使用了coefplot，最好引用一下這篇文章
ssc install coefplot 報錯了！
cannot write in directory C:\Users\���\ado\plus\c
解決辦法：https://bbs.pinggu.org/thread-10685955-1-1.html
每次記得運行一下這個do.file 就可以使用ssc install 命令了
第9行
圖上給出了這幾個的β estimate 以及95% CI
但在實際中，我們很少關注截距是怎樣的，而且截距的范圍很大（-5000，20000）。導致mpg trunk等置信區間被壓縮成一個點
第13行 coefplot, drop(_cons)
第15行 coefplot, keep(trunk turn)
第32行 coefplot D F, drop(_cons) xline(0) 缺點：legend不清楚D、F代表什么 給它們加label
灰線和黑線均勻分布在縱坐標對應刻度的兩邊，兩條線的距離可以更改
總結
以上是生活随笔為你收集整理的医咖会stata 笔记（自己能看懂版的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 幼儿园小班上计算机课 作业内容是手口一致
- 下一篇： 《Did I Buy the Wrong