【序言】
上周我的小文章,無意中引來了許多師長和學友的關注,這是我始料未及並受寵若驚的。平心而論,我的初衷已經實現了:那就是向學友們說明,量化研究在實踐中存在著極複雜的數據清理,使得嚴格意義上的「可複製性」只存在於理論可能性。要真正實踐這一要求,原作者對操作過程的詳細說明(可用附錄形式)或公開代碼是必要的。
我的文章重點想證實一個觀察:在目前的學術條件下,由於量化研究對數據清理的複雜操作,可複製性已經在事實上不可能。我身體力行,選擇了一篇論文實證之。然而,因為兩個因素,我上次操作的效果打了折扣。第一,由於精力與水平所限,我的複製出現了兩個不該出現的錯誤。雖然這兩個錯誤對複製失敗的責任是有限的(下面詳說),但觀者難免覺得我的實踐難以令人信服。第二,原作者不可避免地認為我在質疑他的研究,要同我一辯是非,過多精力放在了探討作者的研究設計方面,這不是我的關切所在。
對於第二點,我(再次)澄清:我做上篇文章的初衷,並不是要證明原作者的研究存在問題;而是藉以表現複製一項量化研究,是何其困難。本文與原作者的結論差異,只說明複製研究非常困難,不能說明作者的研究存在問題。本來在最初的複製文章開頭我加了這句話,但由於騰訊公眾號糟糕的同步設置,這句話沒有出現在最終版本中,待我發現時已追之不及。我的邏輯前提,就是假設原研究沒有錯誤,而複製者卻仍然無法加以復現。在我複製的過程中,確實發現了作者少量的(疑似)誤操作,但這僅僅是本次實驗的副產品。
對於第一點,我的兩處重大錯誤是毋庸諱言的。一,在2002年的數據處理中,對「住房類型」的分類採納了錯誤的解釋;二,在2014年的「住房建築面積」變量中,沒有合併前一期數據的相關變量,導致這一變量缺失值太多。
但我要為自己做兩點辯護。首先,雖然錯誤很明顯,但這兩個錯誤不是複製失敗的主要原因。第一個錯誤會導致表4中2002年結果異常,但對表5和表6不會有影響;第二個錯誤導致我放棄了2014年數據中「住房建築面積」這個變量。但這樣做相當於2014年沒有測量住房建築面積。但理論上說,這不會對我的其他結論產生影響。社會世界不是量子王國,「住房資本」或「住房套數」也不是薛丁格的貓,不會因為測量或不測量某個變量而產生變化。
第二,我犯這兩個錯誤,本身就是複製研究之難的一個體現。第一個錯誤是編碼本錯漏造成的:CHIP2002的編碼本沒有解釋住房類型1-6的具體含義。我雖然誤信網上信息,但也是由於網上的解釋與前一期的類似變量的分類一致,我才予以採信。第二個錯誤當然主要在我自己犯傻。但是,複製研究作為義務勞動,一定是在短時間內完成的。我在高強度的數據處理過程中,雖然看到了相關變量,但一念之差,沒有深究。對於CFPS這樣結構複雜的數據來說,維護方提供的編碼本在用戶體驗上差強人意。學術活動是受社會生活和物質環境制約的,研究者的時間精力、研究工具的便利與否、乃至於研究環境的物理特徵(比如防火牆對網際網路的限制),都有可能對研究本身產生影響。
我的觀點是,之所以複製研究很難,是因為當代社會科學的數據處理已經成為了一個高度複雜的系統。但與理工科實驗相比,這項工作的自動化程度還很低,要求的人工主觀判斷還很多。對於類似吳老師這樣的研究來說,研究者在數據清理過程中可能需要做出上百個判斷。這些判斷有些固然是深思熟慮的後果,有些卻只是習慣使然、一念之間。即使有一定規則制約,多少也因人而異。不同研究者要在這些問題上獨立做出完全一致的判斷,其概率微乎其微。受期刊篇幅的限制,變量說明不可能如下文一般長篇累牘。不要說吳老師原文不到一頁的變量說明,即使是四五頁的詳細說明,也難免會有涉及不到的地方。
吳老師在回覆中再三強調量化研究要「嚴謹」。嚴謹當然是對的,做什麼研究都要嚴謹。但是我想說,對於這樣一個超複雜結構,並非「嚴謹」就可以實現研究的可複製。藉此我想說明,量化社會科學正在面臨的可複製危機,一部分是由于越來越複雜的變量操作造成的,只有非常詳細的變量說明或公開代碼,才有可能逐步克服這一問題。吳老師暗示我對量化研究有「偏見」。我如果真的對量化研究有偏見,也就不會做這件事了。但我確實希望說明:在實證主義的外衣之下,社會科學的量化研究是一個充斥研究者主觀性、受限於研究者能動性和研究環境的活生生的過程。承認這一點是真正的科學態度的體現。我因而把我的工作理解為對量化社會學研究的知識社會學探索。
我既將我的行為看做行為藝術,也看做一次知識社會學的實驗。作為行為藝術,我想向讀者展示:在數據處理的過程中,一個研究者需要進行多少細小的主觀判定。作為實驗,我希望檢驗:在有了正常標準的操作說明後(作者原文解釋略簡單,加上對拙文的回覆可算正常標準的操作說明),在眾多讀者甚或作者本人的監督和指導下、在極盡我個人之所能之「嚴謹」後,究竟能不能復現這一研究。為此,我將主要的代碼都展現出來,感興趣的讀者只需要複製粘貼,即可在自己的電腦上重現我的操作,以此來監督我做這項實驗,也限制我個人水平這個外生變量。我完全相信在這一次的複製中,我仍然犯了錯誤、甚至是重大錯誤。我歡迎任何人向我不吝指出,這不是對我的冒犯,反而恰恰支持了我的論點。
最後,請吳老師千萬不要認為我在針對您個人。固然,我對吳老師的部分操作有所保留,但這只是本次實驗的副產品。我對吳老師個人非常尊重,對這次給他造成的不便也深感歉疚。他日相逢滬上,當浮一大白謝罪。
是為序。
【說明】
1. 正文部分由以下部分組成
【0. R環境說明】
【1. CFPS2010的數據清理】
【2. CFPS2010的樣本選擇和研究分析單位問題】
【3. CFPS2014的數據清理】
【4. 基尼係數】
【5. 泰爾係數】
【6. 住房套數不平等】
【7. 住房面積不平等】
【8. 住房資產不平等】
【9. 評估複製效果】
2. 提要
a. 在得到原作者較為詳細的操作說明,並修正了先前若干錯誤之後,得到了比較接近原文的描述統計。但在如何處理退休人員、如何選取個體樣本等問題上,仍與原文存在較大差別。
b. 仍然沒有復現作者關於2014年住房不平等超過收入不平等的論斷。
c. 作者的泰爾係數計算可能存在問題。復現結果與原文相反。
d. 在模型部分,復現了作者29個主要觀察中的20個。復現了作者關於住房套數不平等與面積不平等的主要結論,但是沒有復現作者關於資產不平等在金融化時期擴大的結論。3. 關於工具
吳老師已經在回覆中說明他是用STATA完成的這項研究。照理說我也應該使用STATA。但重頭做起確實太花時間,並且R作為面向對象的語言,無疑更適合處理這類需要大量合併和拆分的數據。更適合也意味著更少的出錯可能。除非特殊說明,本次複製的代碼為R語言。本文完整代碼見:
https://github.com/duplicate-soc/replicator/blob/master/Replicate%20Wu2019
4. 關於1988年和2002年的數據
在我的初次複製中,1988年數據的結果與原文比較接近,吳老師似乎也比較認可。但其中有一處小錯誤,我略作修改,沒有對數據產生太大影響。2002年的數據除了在住房套數上出現錯誤,似乎也沒有被過多批評,我檢查之後也決定不做其他修改。下文不再贅述這兩期數據的處理。
【0. R環境說明】
本文使用以下命令包。
library(pacman)p_load(tidyverse, gtools, psych,stringr, gmodels, AER, readstata13, ineq, foreign, VGAM,DescTools)
【1. CFPS2010的數據清理】
CFPS2010數據由五個部分組成:家庭關係庫、家庭問卷庫、成人問卷庫、少兒問卷庫、社區問卷庫。我認為原文使用了其中的三個部分:家庭問卷庫、成人問卷庫、社區問卷庫。原文對如何選擇城市樣本沒有解釋,但作者在回復拙文時說:
「在城鄉樣本選擇方面,CFPS有兩個選項,一個是統計局的城鄉分類(變量名稱為urban),另一個是訪問員記錄的社區分類(變量名為cz7)。由於很多劃入城市的區域原來是農村,這些區域很多住房還是宅基地房。為了保證分析的是城市的住房樣本,我剔除了CFPS中的農村、郊區和村改居社區樣本,以城市家庭作為研究對象。研究住房的學者都知道,中國的城鄉住房在產權、價格等方面存在巨大差異,城市樣本的選擇對研究結果會有顯著影響。」
讀入數據使用的是readstata13的包。
# read individual-level datadata2010a <- read.dta13("ecfps2010adult_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")# read household-level datadata2010h <- read.dta13("ecfps2010famecon_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")# read community-level datadata2010c <- read.dta13("ecfps2010comm_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")
錄入後有33600條成年人樣本,14798條家庭樣本。
接下來,上述引文說對樣本的選擇是基於urban這個變量選擇城市(urban =1),然後基於cz7去掉農村(cz7=3)和郊區(cz7=4),最後在剩餘樣本中通過cz701這個變量去掉村改居社區(cz701)。
# select the urban cases# three variables: "urban", "cz7","cz701"# urban == 1, 城市;# cz7 == 1 | 2 城市或城鎮,去掉農村和郊區# cz701 != 9 城市不是村改居# cz702 != 9 城鎮不是村改居data2010c %>% select(cid, cz7, cz701, cz702) -> communitydata2010a <- merge(data2010a, community, all.x= TRUE, by = "cid")data2010a %>% filter(urban == 1) %>% filter(cz7 == 1 | cz7 ==2 ) %>% filter(cz701 != 9) %>% filter(cz702 !=9) -> data2010a
這一輪篩選之後,留下9670個成年人樣本,分屬4478個家庭。
本項研究在個體層面上涉及以下變量:性別(變量名:gender)、世代(qa1y_best)、戶籍(qa2)、黨籍(qa7_s_1)、學歷(cfps2010edu_best)、婚姻(qe1)、單位(qg305)、職業(qg307code)、退休(qg301)、省份(provcd)。這項數據對成年人的標準似乎是16周歲,我覺得這也可以接受,何況數量也不多(小於18周歲226個人),我沒有做處理。我將這些變量中拒絕回答或不知道的選項改為缺失值。
# individual-level cleandata2010a$qe1[data2010a$qe1 < 0] = NA # qe1 is for marriagedata2010a$income[data2010a$income < 0] = NAdata2010a$q1y_best[data2010a$qa1y_best < 0] = NA# QG305 您現在主要是在哪個機構工作,codebook 231頁# delete if qg305 == -1 don't know; == -2 refuse; == -8 unapplicable# == 17 unable to judgedata2010a$qg305[data2010a$qg305 < 0] = NAdata2010a$qg305[data2010a$qg305 == 17] = NA# QG307CODE 職業事後編碼data2010a$qg307code[data2010a$qg307code < 0] = NAdata2010a$qg307code[data2010a$qg307code > 69999] = NAdata2010a$qa1age[data2010a$qa1age < 0] = NA
在生成適合分析的變量的過程中,比較麻煩的職業和單位兩個變量。在原文表1中,可以看到工作單位中的「離退休」與職業中的「退休人員」是不一樣的,數值上「退休人員」時而比「離退休」要高一些(1995;2002;2012),時而低一些(1988;2010;2014)。但是,涉及這個問題的只有一個變量(qg301:您是否已經離/退休(包含病退、內退等非正式退休))。
我是這麼理解的:考慮單位時,離退休人員待遇也會隨單位而有所差別,因此優先將這群人劃入退休前單位,如果沒有寫單位,則單獨列為「離退休人員」;考慮職業時則相反,退休人員優先劃為「退休人員」。但這無法解釋為什麼「退休人員」時而高於「離退休」,時而又偏低。我暫時如此處理,只能等待作者解釋。
data2010a %>% mutate(SEX = gender) %>% mutate(GENERATION = cut(qa1y_best, breaks = c(-Inf,1939.5, 1949.5, 1959.5, 1969.5, 1979.5, Inf), labels = c("Before 1939", "1940-9", "1950-9","1960-9","1970-9","1980"))) %>% # Generation mutate(URESIDENCE = (qa2==3)) %>% # Urban hukou mutate(MARRIAGE = (qe1 == 2)) %>% # Marriage mutate(PARTY = (qa7_s_1 == 1)) %>% # PARTY MEMBER mutate(COLLEGE = (cfps2010edu_best >=5)) %>% # 大專及以上 mutate(INCOME = income) %>% mutate(OCCUPATION = cut(qg307code, breaks = c(9999, 19999, 29999, 39999, 69999), labels = c("Cadre","Professional","Clerk","Worker"))) %>% mutate(PROVINCE = as.factor(provcd)) %>% select(pid, fid, indno, SEX, URESIDENCE, GENERATION, MARRIAGE, PARTY, COLLEGE, INCOME, OCCUPATION, PROVINCE, qg3, qg301, qg305) -> data2010
data2010 %>% mutate(OCCUPATION = as.character(OCCUPATION)) %>% mutate(OCCUPATION = ifelse(qg3==0, "Worker", OCCUPATION)) %>% # "Worker" is short for Worker and Unemployed mutate(OCCUPATION = ifelse(qg301==1, "Retired", OCCUPATION)) -> data2010
# if the case is not working at all, he/she is "non-welfare", which means # being out of the welfare state.# but if that is because he/she has been retired, then he/she is attributed as retired# but if the retired person has reported his/her former workplace, then he/she is # classified into the workplace.data2010$WORKPLACE[data2010$qg3==0] <- "Non-welfare"data2010$WORKPLACE[data2010$qg301==1] <- "Retired"data2010$WORKPLACE[data2010$qg305 == 1 | data2010$qg305 == 13 | data2010$qg305 == 14 | data2010$qg305 == 2] <- "Government"data2010$WORKPLACE[data2010$qg305 == 10 | data2010$qg305 == 11 | data2010$qg305 == 12| data2010$qg305 == 5 | data2010$qg305 == 6 | data2010$qg305 == 7 | data2010$qg305 == 77 | data2010$qg305 == 8 | data2010$qg305 == 9 ] <- "Non-welfare"data2010$WORKPLACE[data2010$qg305 == 3 | data2010$qg305 == 4] <- "SOE or Collective"data2010 %>% select(-c(qg3,qg301,qg305)) -> data2010到這裡,對2010年個體數據的處理就結束了。
下面對2010年家庭數據加以處理。
家庭級別的數據主要涉及到如下變量。家庭人口數(familysize)、當前住房建築面積(fd2)、其他住房建築面積(fd702)、當前住房估值(fd4)、其他房產估值(fd703)、房貸或其他購房借款(mortage,報告CFPS官方:拼錯了)、房屋產權(fd1)、是否有其他房產(fd7)。
這裡(又)有個棘手的問題,住房估值以萬元為單位,房貸卻以元為單位。住房估值中有一些異常的高值,顯然是錯把萬元當做了元。經過觀察(比對住房估價出現高值時的住房面積),我認為可以1000為一條界限,超過1000的都是弄錯單位的。對「其他房產估值」我加一條限定條件:其他住房面積小於200平米。
#household-level data clean# building area of the housedata2010h$fd2[data2010h$fd2 < 0] <- NA # -1 unknown; -2 Refuse# the current value of your house in 10,000 yuan# If house assset is more than 5000, I concern it as in wrong unit of measurementdata2010h$fd4[is.na(data2010h$fd4)] <- -8 # replace NAsdata2010h$fd4[data2010h$fd4 >=1000] <- data2010h$fd4 / 10000 data2010h$fd703[data2010h$fd703 >=1000 & data2010h$fd702 < 200] <- data2010h$fd703 / 10000describe(data2010h$fd703)data2010h$fd4[data2010h$fd4 < 0] <- NAdata2010h$fd702[data2010h$fd702==-8] <- 0 data2010h$fd702[data2010h$fd702<0] <- NA # Area of other placesdata2010h$fd703[data2010h$fd703==-8] <- 0data2010h$fd703[data2010h$fd703<0] <- NA # Value of other places
另一個棘手的問題是產權類型比較複雜。
我採用最嚴格的定義,僅將「完全自有」視為佔有房屋。這當然也是我個人的一個小決定。
data2010h %>% mutate(NHH = familysize) %>% mutate(AREA = (fd2+fd702) / NHH) %>% mutate(ASSET = (fd4 + fd703 - mortage/10000) / NHH) %>% mutate(HOUSEOWN = (fd1 == 1)) %>% # owned solely by family mutate(MULTI = (fd7 == 1)) %>% select(fid, NHH, AREA, ASSET, HOUSEOWN, MULTI) -> data2010household
合併家庭數據與個體數據。
data2010 <- merge(data2010, data2010household, all.x = TRUE, by = "fid")
計算家庭人均收入。
這裡需要指出,吳老師在回復我的過程中說了這樣一段話:
「由於購房是家庭行為,我們的收入採用了家庭人均收入,而不是個體的收入。布文上說他的操作是「選擇每一戶中收入最高的個體,因為他/她的情況最能反映該戶購買住房的水平。」 布文這種方法會存在一些問題,一是容易造成樣本量丟失;二是從實際來看,年輕人的收入可能更高,但購房人可能是他的父母。對個體樣本選擇的不同,是原文與布文在研究結論上存在差異的另一個重要原因。」
吳老師在這裡有意無意地混淆了兩個問題。在計算收入時,我採用的同樣是「家庭人均收入」。引文中說的這個操作,是在處理另一個問題,下節會專門討論。
# aggregate household incomeaggregate(INCOME~fid, data2010, sum, na.action = na.omit) -> householdindata2010 <- merge(data2010, householdin, all.x= TRUE, by = "fid")data2010$INCOME_A <- data2010$INCOME.y/data2010$NHH # Average Income Within a Household
最後增加一個表示房產套數的變量,1表示一套,2表示兩套及以上。
至此,對CFPS2010的數據清理結束。新生成的數據包含個體樣本9670人,分屬4478個家庭。在處理過程中沒有損失變量。
【2. CFPS2010的樣本選擇和研究分析單位問題】
真正棘手的問題在於,原文是如何從上述9670個個體(假設作者操作與我一樣)中選出最終納入模型的3325個的。
這裡涉及研究設計中一個複雜的問題,就是這一研究的分析單位既是家庭(人口數、人均收入),也是個體(世代、職業、單位、性別、黨籍)。這種情況下我會傾向於使用多層模型(multilevel model)。但既然作者沒有用,這個問題暫時不作展開。
我們看到,作者事實上是以個體為分析單位,而家庭的作用是展現在個體中的。這樣的話,就有必要選出一個具有代表性的家庭成員。比如說為了展現富豪家庭的收入情況,非要以個體來考慮的話,不應該選擇一位不工作的主婦。
岔開說一句,我覺得作者以個體為分析單位的做法有生態謬誤的風險,會將家庭層面的特徵與個體特徵混淆起來。比如,我們都在大城市見過一類夫妻:一方在體制內爭取福利、解決戶口;另一方在體制外努力掙錢、創造流動性。這類家庭在住房等問題上的優勢其實是兩者組合的後果,即家庭後果,不是某一方的特徵造成的。
拋開這點不論,我最初的想法是,從家庭中選出收入最高的那位成員作為代表。這麼做,也與前三年CHIP數據中以男性戶主為主要分析單位的情況保持了一致。
作者在回復我的過程中雖然誤解了我的意思,但至少否認了自己是如此操作的。那麼,我只能想出另一種方法:選擇每戶家庭的第一位受訪人。這種處理不能保持與前三年的一致,也無法客觀反映該家庭的收入情況。
我將兩種方式都實踐了一下。
# select the family members with highest income as the headdata2010 <- data2010[order(data2010$fid, -data2010$INCOME.x),]data2010 %>% group_by(fid) %>% slice(1) %>% ungroup() -> df2010_income# select the family members ranked as the firstdata2010 <- data2010[order(data2010$fid, data2010$pid),]data2010 %>% group_by(fid) %>% slice(1) %>% ungroup() -> df2010_order# drop cases with nasdf2010_income <- na.omit(df2010_income)df2010_order <- na.omit(df2010_order)
生成了兩個新數據,按我的第一種思路形成的df2010_income有3145個樣本,按第二種思路形成的df2010_order有3138個樣本。都與原文的3325個樣本接近(雖然還是不同的)。
這一次比之前接近很多,可見原作者在回覆中提供的信息還是很重要。我們處理退休人員的手法還是有差異。從男性比例來看,我傾向於認為作者是使用了類似第二種思路的手法,但不排除作者有第三種選擇的可能。
【3. CFPS2014的數據清理】
這一部分基本與前面一樣。但是感謝原作者在回復拙文時的指正,部分變量需要與前期數據合併。
CFPS2014的數據組成與CFPS2010一樣。首先我們需要讀入成人資料庫、家庭資料庫、社區數據區,以及供匹配所用的CFPS2012家庭資料庫。
data2014a <- read.dta13("ecfps2014adult_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")
# read household-level datadata2014h <- read.dta13("ecfps2014famecon_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")
# read community-level datadata2014c <- read.dta13("ecfps2014comm_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")
# read cfps2012 household-level datadata2012h <- read.dta13("ecfps2012famecon_201906.dta",convert.factors = TRUE, generate.factors = FALSE, encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE, missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE, add.rownames = FALSE, nonint.factors = FALSE, select.rows = NULL, select.cols = NULL, strlexport = FALSE, strlpath = ".")然後選擇按照cz7和cz701這兩個變量選擇樣本,步驟與上一節基本類似。
# Select the urban sample, confirmed by the authordata2014c %>% select(cid14, cz7, cz701) -> communitydata2014a <- merge(data2014a, community, all.x = TRUE, by = "cid14")data2014a %>% filter(urban14 == 1) %>% filter(cz7 ==1 | cz7 == 2) %>% filter(cz701 !=9) -> data2014a
一共6495個個人,分屬2758個家庭。
下面清理個人層面的數據。
data2014a$p_income[data2014a$p_income<0] <- NA # Set missing value for household income as missing value
data2014a %>% mutate(SEX = cfps_gender) %>% mutate(BORN = cfps_birthy) %>% mutate(GENERATION = cut(BORN, breaks = c(-Inf,1939.5, 1949.5, 1959.5, 1969.5, 1979.5, Inf), labels = c("Before 1939", "1940-9", "1950-9","1960-9","1970-9","1980"))) %>% # Generation mutate(URESIDENCE = (qa301 == 3)) %>% mutate(MARRIAGE = ((qea0 == 2))) %>% # at Marriage mutate(PARTY = (cfps_party == 1)) %>% # PARTY mutate(COLLEGE = (cfps2014sch>=5 & cfps2014sch<=8)) %>% # college or not mutate(INCOME = p_income) %>% mutate(OCCUPATION = cut(qg303code, breaks = c(9999, 19999, 29999, 39999, 69999), labels = c("Cadre","Professional","Clerk","Worker"))) %>% mutate(PROVINCE = as.factor(as.character(provcd14))) %>% select(pid, fid10, SEX, URESIDENCE, GENERATION, MARRIAGE, PARTY, COLLEGE, INCOME, OCCUPATION, PROVINCE, employ2014, qg2, qi101) -> data2014
data2014 %>% mutate(OCCUPATION = as.character(OCCUPATION)) %>% mutate(OCCUPATION = ifelse(employ2014 == 0 | employ2014 == 3, "Worker", OCCUPATION)) %>% mutate(OCCUPATION = ifelse(qi101==1, "Retired", OCCUPATION))-> data2014
data2014$WORKPLACE[data2014$employ2014 == 0 | data2014$employ2014 == 3] <- "Non-welfare"data2014$WORKPLACE[data2014$qi101==1] <- "Retired"data2014$WORKPLACE[data2014$qg2==1] <- "Government"data2014$WORKPLACE[data2014$qg2== 2 | data2014$qg2 == 3] <- "SOE or Collective"data2014$WORKPLACE[data2014$qg2>=4] <- "Non-welfare"
select(data2014, -c(employ2014,qg2,qi101)) -> data2014個人數據清理完畢。
下面處理家庭層面的數據。
上次我犯了錯誤,在2014年的住房面積中需要併入2012年的數據。這次改過來。
CFPS2012中最好的建築面積變量,我認為是fq701_best;關於房產的信息,相關變量是產權(fq1;1、2)、其他房產(fr101);關於價值的變量時現住房市價(houseprice1_best)、其他房產的總價值(fr2test)
CFPS2014數據中,住房面積變量取決於兩個變量,先問住房面積是否變化(FQ8),如果是再確認新的住房面積(FQ801)。對於房產信息和房產價值卻不是這麼處理的,房產信息只問了是否還有房產,沒有問「房產產權是否變化」,涉及現住房產權的變量是fq2,涉及是否有其他房產的變量時fr1。對於房產價值,現房價值修正值變量時fq6_best,其他房產市價是fr2_best。所以,在CFPS2012中只需要涉及建築面積的變量,即fq701_best。
首先,從CFPS2012中提取房產面積變量,加以清理之後合併入2014年的家庭數據中。生成一個表示房屋面積的變量,如果房屋面積沒有變化,那麼面積就等於2012年時的面積;如果有變化,就等於2014年時的面積。
# select the cases from data2012h for constructive areadata2012h %>% select(fid10, fq701_best) -> area2012colnames(area2012) <- c("fid10","area2012")area2012$area2012[area2012$area2012 <0] = NAdata2014h <- merge(data2014h, area2012, all.x = TRUE, by = "fid10")
data2014h$area2014 <- ifelse(data2014h$fq8==1, data2014h$area2012, data2014h$fq801)CFPS2014對房產價值的數據比較整齊,我決定不對1000萬以上的房價做修正。這一年的房屋所有權問題上出現了新提法,家庭成員擁有「完全產權」或「部分產權」。為了與前文保持一致,我只計算擁有「完全產權」的情況。
data2014h$area2014[data2014h$area2014 < 0] <- NA # the difference between refusal and unknown is no longer importantdata2014h$fq6_best[data2014h$fq6_best==-8] <- 0data2014h$fq6_best[data2014h$fq6_best==-1] <- NAdata2014h$fr2_best[data2014h$fr2_best==-8] <- 0data2014h$fr2_best[data2014h$fr2_best==-1] <- NA
data2014h %>% mutate(NHH = fml2014num) %>% mutate(AREA = area2014 / NHH) %>% mutate(ASSET = (fq6_best + fr2_best -mortage/10000) / NHH) %>% mutate(HOUSEOWN = (fq2 ==1)) %>% mutate(MULTI = (fr101>=1)) %>% select(fid10, NHH, AREA, ASSET, HOUSEOWN, MULTI) -> data2014household
data2014 <- merge(data2014, data2014household, by = "fid10")照舊產生一個表示家庭人均收入的分組變量和一個表示房屋產權數量的變量。
# aggregate household income and divided into four groupsaggregate(INCOME~fid10, data2014, sum, na.action = na.omit) -> householdindata2014 <- merge(data2014, householdin, by = "fid10")data2014$INCOME_A <- data2014$INCOME.y/data2014$NHHdata2014 %>% mutate(INCOMEGROUP = quantcut(INCOME_A,4, na.rm=TRUE)) %>% mutate(INCOMEGROUP = factor(INCOMEGROUP, labels = c("1st","2nd","3rd","4th"))) -> data2014
# generate a variable for house ownershipdata2014$HOUSE <- "0"data2014$HOUSE[data2014$HOUSEOWN == TRUE] <- "1"data2014$HOUSE[data2014$MULTI == TRUE] <- "2"
# change province into factordata2014$PROVINCE <- as.factor(data2014$PROVINCE)
最後,仍然以上一節說的兩種策略生成兩個資料庫。
# select cases by incomedata2014 <- data2014[order(data2014$fid, -data2014$INCOME.x),]data2014 %>% group_by(fid10) %>% slice(1) %>% ungroup() -> df2014_income# select cases by orderdata2014 <- data2014[order(data2014$fid, data2014$pid),]data2014 %>% group_by(fid10) %>% slice(1) %>% ungroup() -> df2014_order
按收入選擇的數據只有1797個樣本,按順序選擇的數據有1875個樣本。與原文的2086個相比差距還是不大的。缺失值主要出現在職業、工作兩個變量上。這是我目前與原文差別最大的地方,但我目前無法使結果更接近。
比較兩組數據與原表數據,差別比2010年要大一些。主要差異仍然在處理體制外工人、離退休人員等問題上。但相比而言,使用自然應答順序得到的數據(df2010_order)還是更接近原表數據。
以上完成了2014年的數據清理。
總體而言,目前四期數據的描述統計已經與原文非常接近。
【4. 基尼係數】
下面展示從四組數據中計算的基尼係數。原作者在回復我時曾說,
「布文2010年的基尼係數為0.469,2014年就上升到0.631,4年時間基尼係數上升了0.162。從變化規律來看,沒有任何國家的基尼係數能在4年時間內上升0.162,這只能是數據處理或計算上出了差錯」
我覺得「只能是數據處理或計算上除了差錯」的說法不嚴謹。因為2010年和2014年兩期數據沒有進行匹配,應該視為兩個樣本(不然原文在解釋上又有新的問題了)。既然是兩個樣本,兩者的差距就有可能是樣本誤差造成的。基尼係數既然基於樣本數據,當然也存在標準誤和置信區間,出現較大提升並不奇怪。所以,我在下面僅展示置信度95%的置信區間。
基尼係數的計算,我用的是DescTools這個包。
library(DescTools)Gini(df1988$INCOME.y/df1988$NHH, conf.level = 0.95, type = "basic")Gini(df2002$INCOME.y/df2002$NHH, conf.level = 0.95, type = "basic")Gini(df2010_income$INCOME.y/df2010_income$NHH, conf.level = 0.95, type = "basic")Gini(df2014_income$INCOME.y/df2014_income$NHH, conf.level = 0.95, type = "basic")Gini(df2010_order$INCOME.y/df2010_order$NHH, conf.level = 0.95, type = "basic")Gini(df2014_order$INCOME.y/df2014_order$NHH, conf.level = 0.95, type = "basic")
Gini(df1988$AREA, conf.level = 0.95, type = "basic")Gini(df2002$AREA, conf.level = 0.95, type = "basic")Gini(df2010_income$AREA, conf.level = 0.95, type = "basic")Gini(df2014_income$AREA, conf.level = 0.95, type = "basic")Gini(df2010_order$AREA, conf.level = 0.95, type = "basic")Gini(df2014_order$AREA, conf.level = 0.95, type = "basic")
對於收入和面積的基尼係數比較,我在下表中列出來了。家庭人均收入的基尼係數在我的兩個數據裡,增速都比原文要快很多。
我的上一篇複製報告忘記指出,基尼係數也是要求非負的。我當時之所以能夠計算,是因為刪除了資產中的負值。作者在回復我時說:
「在數據處理的時候發現,只有兩個樣本的住房價格(值)減去住房貸款後的數值為負數,可以忽略不計。為了嚴謹起見,我在數據分析中刪除了住房淨資產為負數的樣本。」
但是,在我複製的時候,資產為複製的不止兩個樣本。比如2002年的數據中有70個家庭資產為負,這裡面有42個是因為報告住房資產為0(這是原始數據,不是我的重新編碼),確有不實的可能性,但也有28個家庭是實實在在報告資產為負的。比如編號3212830031的家庭,報告自家房產價值8萬,購房借貸15萬;編號3403000078的家庭,報告自家房產價值1.5萬,借貸4萬;等等。作者固然可以說,這可能是受訪者對自己房產估值有誤。但是,這裡每個數據都是基於受訪者的自己估值。如果懷疑這一點,是否就應該對整個數據處理加以修改?不能說資產為正,估值就準確;資產為負,估值就不準確吧?
報告負資產家庭的代碼如下,以2002年為例。
df2002%>% filter(df2002$ASSET < 0) %>% select(pcode)-> negative2002data200202 <- read.table(file = '21741-0002-Data.tsv', sep = '\t', header = TRUE)data200202 %>% select(pcode, h414, h417) -> asset2002negative2002 <- merge(negative2002, asset2002, all.x = TRUE, all.y = FALSE, by = "pcode")summary(negative2002$h414==0)
CFPS的幾期數據中負值確實不多,2010年有兩例,2014年有四例。
因此我退一步,接受作者的說法,刪除這些資產為負的樣本。
與作者的結果挺接近的,但還是不能說住房資產不平等加劇得比收入不平等更嚴重。
【5. 泰爾係數】
在處理泰爾係數的問題上,原文有一處確定的錯誤和一處潛在的錯誤。
確定的這處錯誤是完全可以理解的。作者在表3中報告人均住房面積組間不平等的佔比時,將「世代」和「收入組」搞反了。看起來,這兩行只是順序在上下兩部分中顛倒了一下,但仔細看,「收入組佔比」那一行應該是「世代佔比」,反之亦然,看數字的大小關係就可以看出來。不過他在行文解釋時沒有說反。我明白這近乎吹毛求疵。我只是想說,有些錯誤並不是作者嚴謹就可以避免的,也不是責任編輯和匿名評審細心就可以發現的,量化研究的細節之多,不犯錯真的是不可能的。
第二個潛在的錯誤就比較嚴重了,但我僅僅是猜測。
我在之前的複製文章中寫過泰爾係數的公式:
這個公式是我從維基百科抄來的,作者在97頁上的公式(2)與這個是一樣的。再抄一個分解公式:
這兩個公式都要求xi是正數。因為這是對數運算,所以要求要高於基尼係數(基尼係數要求非負數)。
我在R上嘗試了兩種計算泰爾係數的函數:ineq包中的Theil函數和IC2包中的decompGEI函數。decompGEI函數不允許變量取值為0,但Theil函數卻允許變量取值為0。為什麼呢?因為Theil函數默認刪去取值為0的樣本。
作者在回應我時,只說房產淨值為負數怎麼辦,卻沒有說明房產淨值為零怎麼處理。在我重製的2002年數據中,資產為0的樣本有955個;2010年15個;2014年320個。作者的描述統計中也有相當數量沒有房產的樣本。所以房產淨值為0是一個需要重視的問題。
對於Theil函數(或類似的算法)來說,有沒有零是不會對結果產生影響的。下面我生成兩組數據,一組包括0,一組不包括0,用Theil函數算,結果是一樣的。
df2002[df2002$ASSET >= 0,] -> df2002tdf2010_income[df2010_income$ASSET >= 0,] -> df2010tdf2014_income[df2014_income$ASSET >= 0,] -> df2014t
df2002[df2002$ASSET > 0,] -> df2002pdf2010_income[df2010_income$ASSET > 0,] -> df2010pdf2014_income[df2014_income$ASSET > 0,] -> df2014p
Theil(df2002t$ASSET) # 0.3777506Theil(df2010t$ASSET) # 0.8321001Theil(df2014t$ASSET) # 0.6915311
Theil(df2002p$ASSET) # 0.3777506Theil(df2010p$ASSET) # 0.8321001Theil(df2014p$ASSET) # 0.6915311我為什麼說原文有可能用了類似的算法呢?因為這個結果與作者的結果很接近。2002年時0.378 vs 0.377;2010年時0.832 vs. 0.826; 2014年時0.692 vs. 0.698。
在討論不平等時,房產為0的群體當然應該考慮。如何考慮呢?可以把房產為0改成房產為近乎0的正數(如0.0001)。使用以下代碼來分解泰爾係數。其他泰爾係數也以相似方法處理。
df2002t$ASSET[df2002t$ASSET==0] <- 0.0001df2002t$COLLEGE <- as.factor(df2002t$COLLEGE)decompGEI(x=df2002t$ASSET, z=df2002t$COLLEGE, w = NULL, alpha = 1, ELMO = TRUE)decompGEI(x=df2002t$ASSET, z=df2002t$WORKPLACE, w = NULL, alpha = 1, ELMO = TRUE)decompGEI(x=df2002t$ASSET, z=df2002t$OCCUPATION, w = NULL, alpha = 1, ELMO = TRUE)decompGEI(x=df2002t$ASSET, z=df2002t$GENERATION, w = NULL, alpha = 1, ELMO = TRUE)decompGEI(x=df2002t$ASSET, z=df2002t$INCOMEGROUP, w = NULL, alpha = 1, ELMO = TRUE)decompGEI(x=df2002t$ASSET, z=df2002t$PROVINCE, w = NULL, alpha = 1, ELMO = TRUE)
這樣,我們在人均住房面積方面的結論是一致的,但是在人均住房資產方面,我得到了與原文完全不同的結果。
作者說:
「在住房市場戶進程中,不同省份、收入和職業群體間住房資產泰爾係數增長明顯,組間不平等也呈上升趨勢。」(101頁)
在我的重製數據中,不同單位、職業和收入組之間的組間不平等恰恰呈現出下降趨勢。省際不平等確實是顯而易見。
【6. 住房套數不平等】
下面使用四個數據重製表4的廣義定序Logit模型。
作者在回復時質疑我使用的STATA命令。oglm和作者使用的gologit2隻是同一作者(Richard Williams)開發的兩個不同版本。具體可參看作者的以下說明:https://www3.nd.edu/~rwilliam/gologit2/。
這一模型相對少見。我吃不準相關的R命令包是否採取了完全相同的解法,因此與作者一樣使用了STATA的gologit2命令。
對於兩期CFPS,我只使用按收入選取的樣本(df_income),因為我認為這樣操作相對更合理。
先在R上使用類似下面的命令將R數據整理並存儲為STATA數據。
# produce stata datadf2014_income$HOUSE <- as.factor(as.character(df2014_income$HOUSE))df2014_income$WORKPLACE <- as.factor(as.character(df2014_income$WORKPLACE))df2014_income$OCCUPATION <- as.factor(as.character(df2014_income$OCCUPATION))df2014_income %>% mutate(GENERATION = factor(GENERATION, levels = c("1960-9","Before 1939", "1940-9","1950-9","1970-9", "1980"))) %>% mutate(WORKPLACE = factor(WORKPLACE, levels = c("Non-welfare", "Government", "SOE or Collective", "Retired"))) %>% mutate(OCCUPATION = factor(OCCUPATION, levels = c("Worker", "Cadre","Professional","Clerk","Retired"))) -> df2014_income
df2014_income <- na.omit(df2014_income)write.dta(df2014_income, "df2014_income.dta")然後在STATA上運行類似下面的命令。
// install gologit2findit gologit2// for CHIP2002 CFPS2010 CFPS2014xi: gologit2 HOUSE SEX i.GENERATION URESIDENCE PARTY COLLEGE MARRIAGE NHH i.INCOMEGROUP i.WORKPLACE i.OCCUPATION i.PROVINCE// for CHIP 1988xi: gologit2 HOUSE SEX i.GENERATION PARTY COLLEGE NHH i.INCOMEGROUP i.WORKPLACE i.OCCUPATION i.PROVINCE
結果展示如下表。
在這張表中,我只展示我自己複製的結果。但用不同顏色標明與原文結果的相似程度。綠色表示接近。
作者在102頁《住房套數不平等機制分析》中的主要結論:
(1.1)產權化和產業化早期(1980-2002),高學歷者住房優勢增強。(1.2)產權化改革(1980-1998)時期,公有制從業者住房優勢增加。
(1.3)金融化階段(2009-),高學歷、高收入、高職業地位群體和黨員有優勢。
(1.4)婚姻狀況、家庭人口數、家庭收入的影響越來越顯著。
我的複製數據支持了結論(1.1)和結論(1.2)。
對於結論(1.3),支持「行政企事業單位管理者」和「黨員」有優勢;不支持「高學歷」的優勢。高收入的優勢由顯著變為不顯著。
對於結論(1.4),家庭人口數的影響一直顯著,但不確定效應是否一直增加(「越來越顯著」應該是指效應增加)。婚姻狀況和家庭收入的影響均由顯著變為不顯著。
【7. 住房面積不平等】
下面對表5的住房面積重新進行OLS回歸分析。
原作者對面積進行標準化,以便於橫向比較係數大小。
這部分的代碼入下。
# Table 5: OLS Regression for Arealibrary(jtools)df1988$AREA.S <- standardize(df1988$AREA)model1988 <- lm(AREA.S ~ SEX + GENERATION + PARTY + COLLEGE + NHH + INCOMEGROUP + WORKPLACE + OCCUPATION + PROVINCE, data = df1988)summ(model1988, digits = getOption("jtools-digits", 3))
df2002$AREA.S <- standardize(df2002$AREA)model2002 <- lm(AREA.S ~ SEX + GENERATION + URESIDENCE + PARTY + COLLEGE + MARRIAGE + NHH + INCOMEGROUP + WORKPLACE + OCCUPATION + PROVINCE, data = df2002)summ(model2002, digits = getOption("jtools-digits", 3))
df2010_income$AREA.S <- standardize(df2010_income$AREA)model2010 <- lm(AREA.S ~ SEX + GENERATION + URESIDENCE + PARTY + COLLEGE + MARRIAGE + NHH + INCOMEGROUP + WORKPLACE + OCCUPATION + PROVINCE, data = df2010_income)summ(model2010, digits = getOption("jtools-digits", 3))
df2014_income$AREA.S <- standardize(df2014_income$AREA)model2014 <- lm(AREA.S ~ SEX + GENERATION + URESIDENCE + PARTY + COLLEGE + MARRIAGE + NHH + INCOMEGROUP + WORKPLACE + OCCUPATION + PROVINCE, data = df2014_income)summ(model2014, digits = getOption("jtools-digits", 3))
作者結論見原文102頁和106頁。
複製結果見下表。
作者主要結論如下:
(2.1) 1988-2002年國有機構對住房條件的影響由不顯著而顯著。
(2.2) 管理精英和專業精英在1988-2008年期間都有顯著優勢。
(2.3)辦事人員在1988-1998年期間有優勢,到1999年之後沒有優勢。
(2.4)高收入群體在1988-2002年期間優勢增加。
(2.5)1940-1949世代在1995和2002年間有優勢。
(2.6)黨員的優勢在1988-2002年期間增加。
(2.7)在2009年之後,學歷、世代、收入、職業造成的不平等持續存在。(2.8) 其中,學歷優勢在擴大。
我的複製結果支持了大部分上述結論(2.1-2.6)。
但是2009年之後學歷的優勢不顯著,更不擴大。職業造成的不平等僅限於普通工人或失業者與「管理精英」之間的差別,確實存在,但有所縮小。
【8. 住房資產不平等】
最後是對住房資產的廣義Tobit模型分析。
Tobit模型中,作者預設了住房資產非負(左刪失)。我對此保留一定意見,但尊重作者的選擇。與作者(可能)不同,我沒有刪去資產為負的樣本,而是將這些樣本設置為0。
在模型實現方面,我使用了VGAM包。並且直接將人民幣換算為2014年的幣值。代碼如下:
# Table 6: generalized Tobit regression# package use: VGAM
# transfer CNY into 2014CNYdf2002$ASSET2014 <- df2002$ASSET * (1.0113*1.0382*1.0178*1.0165*1.0482*1.0593*.9927*1.0318*1.0555*1.0262*1.0263*1.0192)df2010$ASSET2014 <- df2010$ASSET * (1.0555*1.0262*1.0263*1.0192)
df2002$ASSET2014 <- ifelse(df2002$ASSET2014 < 0, 0, df2002$ASSET2014)summary(tobit2002 <- vglm(ASSET2014 ~ SEX + GENERATION + URESIDENCE + PARTY + COLLEGE + MARRIAGE + NHH + INCOMEGROUP + OCCUPATION + WORKPLACE + PROVINCE, tobit(Lower = 0), data = df2002))
df2010$ASSET2014 <- ifelse(df2010$ASSET2014 < 0, 0, df2010$ASSET2014)summary(tobit2010 <- vglm(ASSET2014 ~ SEX + GENERATION + URESIDENCE + PARTY + COLLEGE + MARRIAGE + NHH + INCOMEGROUP + OCCUPATION + WORKPLACE + PROVINCE, tobit(Lower = 0), data = df2010))
df2014$ASSET2014 <- ifelse(df2014$ASSET < 0, 0, df2014$ASSET)summary(tobit2014 <- vglm(ASSET2014 ~ SEX + GENERATION + URESIDENCE + PARTY + COLLEGE + MARRIAGE + NHH + INCOMEGROUP + OCCUPATION + WORKPLACE + PROVINCE, tobit(Lower = 0), data = df2014))順便說,在這裡上次我犯了一個非常尷尬的錯誤。因為太尷尬了我就不詳述了。要不犯錯真的難。
回歸的結果非常接近原文,要是每一個模型都可以這樣就好了。(3.1)「房改售房以後,職工擁有了住房產權……體制內單位職工的係數由負轉正」;(3.3)1998-2008年期間,高收入的資產優勢擴大;(3.6)2009年以後,不同職業和收入群體的住房資產差異顯著;(3.7)同時期,最高25%收入階層與其他階層的住房資產分化呈擴大趨勢。複製數據支持3.2-3.5。對於3.2,我無法進行1995年和2002年的比較,但支持2002年管理精英係數顯著這一點。對於3.1,我曾在前一篇文章中批評作者有時候將結論「寫成一些比較籠統的、幾乎不會出錯的論斷,而不是清晰的、可操作的經驗性」論斷。作者似乎沒有理解我這句批評的所指。比如像下面的這段話,很難聯繫到表格內的經驗數據:「在福利住房體系下,職工不擁有居住公房的產權,住房的資產價值無法顯現,體制內單位和高學歷職工的住房資產積累處於『劣勢』。房改售房以後,職工擁有了住房產權,住房資產價值顯現,體制內單位職工的係數由負轉正、管理精英的係數由不顯著變為顯著。」「福利住房體系」應該是指1995年,這一年「體制內單位」係數確實是負的,說「處於劣勢」沒有問題;「房改售房以後,職工擁有了住房產權,住房資產價值顯現,體制內單位職工的係數由負轉正」,可能是指2002年「國有機構」的係數轉為正數,但是這個係數不顯著啊,而且後面三年都是負效應,無法說明房改提高了制度內職工的優勢。既然這些數字沒有說明任何變化,作者為什麼要專門點出來呢?我之前指摘的是這裡,但作者在回復時卻要求我讀110-111頁《結論與討論》部分的觀點,有點奇怪。「在住房金融化階段……不同職業和收入群體的住房資產差異顯著。」這一句話對應的經驗證據是「行政企事業單位管理者」在2010年和2012年顯著,「專業和科研人員」在2010年顯著。但是「管理者」在2014年又不顯著了,而「專業和科研人員」在2012年和2014年都不顯著;辦事人員作為高於「普通職工和失業」的群體,係數一直不顯著。這些問題是不是都值得討論呢?細究起來,是不顯著的情況多、顯著的情況少,為什麼以一句「不同職業和收入群體的住房資產差異顯著」就概括了呢?
對於3.7,從係數看,最高收入組的優勢在2002-2010年間擴大了,但隨後又有所縮小。對於這個問題,作者自己數據的證據就不是很強。我之前提過,2010、2012和2014年應該視作三個獨立的樣本。我們對比原文表6中這三年「最高25%」收入組的係數在95%置信度下的置信區間。## discuss the coefficients of "top 25%"coefficient <- c(13.772, 14.785, 16.084)se <- c(1.946, 2.343, 2.704)year <- c(2010, 2012, 2014)errorbar <- as.data.frame(cbind(coefficient,se, year))library(sfsmisc)errbar(year, coefficient, coefficient+1.96*se, coefficient-1.96*se)原文中這三年的係數,只是「剛好」在遞增狀態下而已。但要說最高收入組在2010年之後優勢比2002年時要大,這沒問題。我的數據也支持這一點(至少2010年是支持的)。這裡我要說幾句批評原作者的話。我認為作者的寫作沒有達到量化研究的清晰性標準。作者再三說自己的基本觀點是住房不平等的變化是階段性的。那麼,在最後總結部分,作者應該告訴讀者,在住房產權化階段(1980-1998),住房產權在職業、收入、單位、世代之間分別是如何變化的,住房面積和住房資產又是如何變化的;在產業化階段(1999-2008)和金融化(2009)階段,這些群體之間的不平等又是如何變化的,各階段的特徵是什麼、規律是什麼、說明了什麼問題。
這段話前面說不平等「具有顯著的階段性特徵」,後面卻給出了不加限定的整體性特徵,使人困惑。
「在市場化各個階段,住房不平等具有顯著的階段性特徵。…… 住房不平等在產權化階段表現為住房產權不平等,在產業化階段表現為住房面積差異性,在金融化階段轉為住房資產不平等。……在住房套數方面,高收入、高職業地位和高學歷者優勢積累效應明顯,年輕世代處於明顯劣勢。在住房面積方面,住房市場化減少了體制內外單位以及不同職業群體間的住房差異……在住房資產方面,住房市場化加劇了不同地區、收入、職業和學歷群體間的住房資產差異,高收入、高職業地位和高學歷者的優勢凸顯,不同地區住房資產差異顯著。」(111頁)
雖然我不足為作者法,但為了說明我所理解的清晰性,我幫作者做了一個表總結他的主要觀察。綠色表示得到了複製數據支持的觀察。作者一共做出了29個觀察(其實有更多其他次要的,這裡不納入了),有20個觀察得到了複製數據的支持。通過觀察該表的每一列,大致能夠認識到為什麼作者說:「住房不平等在產權化階段表現為住房產權不平等,在產業化階段表現為住房面積差異性,在金融化階段轉為住房資產不平等。」前面只是從數量上評估複製的效果。從質上看,作者所說的「住房套數」的不平等主要發生在產權化時期得到了支持;面積不平等主要發生在產業化時期,這一說法本身不是很準確,應該說發生在產權化與產業化早期(主要是1995-2002年),可以說得到了支持。但資產不平等主要發生在金融化時期的說法沒有得到支持。如果竟然有人讀到了這裡……我覺得還有許多延伸的問題可談,甚至可以做一個改進的研究設計。但是這篇文章僅關於對這項研究的複製。我覺得我已經窮盡了我個人對複製一篇論文的時間和精力預算。複製成功率大概是多少呢?我想也就是打六十分吧。