社会科学のためのデータ分析入門 章末問題解答(2章-2) Rコード
これまでの章の解答はこちら
章末問題解答(1章-1)と(1章-2)のRコードの記事はこちらこちら(1章-1)(1章-2)から確認できます。
はじめに (Textbook Solution: Quantitative Social Science: An Introduction )
Rを使った統計学の日本語のテキストとして非常に定評のある社会科学のためのデータ分析入門
欠点なのかは分かりませんが、こちらのテキストには章末問題の解答がついていません。そして日本語でも英語でもwebで公開されていません（2018年冬ごろの時点では）。2018年冬に私が上巻の章末問題を解いたのですが、一度公開してみようと思ったので複数の記事に分けて投稿していこうと思います。誰かの役に立てればとも思っているのですが、私のコードにミスがあった場合に指摘していただけると嬉しいです。
I would highly appreciate if you could point out mistakes.
また同じ変数に関するプロットをする場合でも複数の方法を使ったりもしています。
2章-1 (Chapter2 - Section 1)
スクリプトをベタ張りしています。
## Chapter 2 Causality ## Exercise Solution setwd("~/qss/CAUSALITY") # ご自身のディレクトリを選択 ## ----------------------------------------------------- ## The author of this script uses Japanese-Version QSS. ## ----------------------------------------------------- ## ----------------------------------------------------- ## Section 2 ## Q1 ## Canvasser: ?K???? gay <- read.csv("gay.csv") head(gay) dim(gay) ## Create subset for wave = 1 & study = 1. wave1 <- subset(gay, gay$study == 1 & gay$wave == 1) ## wave1 <- subset(gay, subset = (study == 1) & (wave == 1)) ## same as above sum(gay$ssm) ## no missing value wave1a <- tapply(wave1$ssm, wave1$treatment, mean) ## Q2 ## Create subset for wave = 2 & study = 1. wave2 <- subset(gay, gay$study == 1 & gay$wave == 2) ## Calculate mean for each of smm. wave2a <- tapply(wave2$ssm, wave2$treatment, mean) ## Convert wave2a into data.frame for simplicity's sake. wave1b <- as.data.frame(wave1a) wave2b <- as.data.frame(wave2a) ## Sample Average Treatment Effect for the Treated for Gay & Straight ## Difference-in-Difference Estimators DIDgay1 <- (wave2b[4,] - wave1b[4,]) - (wave2b[1,] - wave1b[1,]) DIDstraight1 <- (wave2b[5,] - wave1b[5,]) - (wave2b[1,] - wave1b[1,]) ## Maybe this part (wave2b[1,] - wave1b[1,]) is not need? ## Please give me feedback on it :D. ## DID DIDgay1 - DIDstraight1 ## Q3 DIDgay2 <- (wave2b[2,] - wave1b[2,]) - (wave2b[1,] - wave1b[1,]) DIDstraight2 <- (wave2b[3,] - wave1b[3,]) - (wave2b[1,] - wave1b[1,]) ## Interpretation: there is no big difference b/w the two ## that means there is a possibility that the bias does not exist. DIDgay2 - DIDstraight2 ## Q4 ## Create subsets for wave = 3:7 & study = 1 wave3 <- subset(gay, gay$study == 1 & gay$wave == 3) wave4 <- subset(gay, gay$study == 1 & gay$wave == 4) wave5 <- subset(gay, gay$study == 1 & gay$wave == 5) wave6 <- subset(gay, gay$study == 1 & gay$wave == 6) wave7 <- subset(gay, gay$study == 1 & gay$wave == 7) ## Calculate mean for each of smm. wave3a <- tapply(wave3$ssm, wave3$treatment, mean) wave4a <- tapply(wave4$ssm, wave4$treatment, mean) wave5a <- tapply(wave5$ssm, wave5$treatment, mean) wave6a <- tapply(wave6$ssm, wave6$treatment, mean) wave7a <- tapply(wave7$ssm, wave7$treatment, mean) ## Convert waves into data.frame for simplicity's sake. wave3b <- as.data.frame(wave3a) wave4b <- as.data.frame(wave4a) wave5b <- as.data.frame(wave5a) wave6b <- as.data.frame(wave6a) wave7b <- as.data.frame(wave7a) ## ATE DIDgay3 <- (wave3b[4,] - wave1b[4,]) - (wave3b[1,] - wave1b[1,]) DIDstraight3 <- (wave3b[5,] - wave1b[5,]) - (wave3b[1,] - wave1b[1,]) DIDgay4 <- (wave4b[4,] - wave1b[4,]) - (wave4b[1,] - wave1b[1,]) DIDstraight4 <- (wave4b[5,] - wave1b[5,]) - (wave4b[1,] - wave1b[1,]) DIDgay5 <- (wave5b[4,] - wave1b[4,]) - (wave5b[1,] - wave1b[1,]) DIDstraight5 <- (wave5b[5,] - wave1b[5,]) - (wave5b[1,] - wave1b[1,]) DIDgay6 <- (wave6b[4,] - wave1b[4,]) - (wave6b[1,] - wave1b[1,]) DIDstraight6 <- (wave6b[5,] - wave1b[5,]) - (wave6b[1,] - wave1b[1,]) DIDgay7 <- (wave7b[4,] - wave1b[4,]) - (wave7b[1,] - wave1b[1,]) DIDstraight7 <- (wave7b[5,] - wave1b[5,]) - (wave7b[1,] - wave1b[1,]) DIDgay3; DIDga4; DIDgay5; DIDga6; DIDgay7 DIDstraight3; DIDstraight4; DIDstraight5; DIDstraight6; DIDstraight7 ## DID DIDgay1 - DIDstraight1 DIDgay3 - DIDstraight3 DIDgay4 - DIDstraight4 DIDgay5 - DIDstraight5 DIDgay6 - DIDstraight6 DIDgay7 - DIDstraight7 ## The difference b/w gay & straight canvasser is still positive ## even in wave 7, and effects seems to be lasted. ## Q5 ## Create subset & calculate means study2 <- subset(gay, subset = (study == 2) & (wave == 1)) study2a <- tapply(study2$ssm, study2$treatment, mean) study2b <- as.data.frame(study2a) ## Randomization seems to be done prorerly ## (because of almost no difference b/w two means). ## Q6 ## Create subset & calculate means study2.2 <- subset(gay, subset = (study == 2) & (wave == 2)) study2.2a <- tapply(study2.2$ssm, study2.2$treatment, mean) study2.2b <- as.data.frame(study2.2a) DIDstudy2 <- (study2.2b[4,] - study2b[4, ]) - (study2.2b[1, ] - study2b[1, ]) ## Compare the results of study 1 & study 2. DIDstudy2; DIDgay1 ## The difference seems not so big. We can conclude that ## the result of study 2 of wave 2 is consistent with study 1. ## Q7 ## Create subsets of different waves study3 <- subset(gay, subset = (study == 2) & (wave == 3)) study4 <- subset(gay, subset = (study == 2) & (wave == 4)) study7 <- subset(gay, subset = (study == 2) & (wave == 7)) ## Create means of different waves study3a <- tapply(study3$ssm, study3$treatment, mean) study4a <- tapply(study4$ssm, study4$treatment, mean) study7a <- tapply(study7$ssm, study7$treatment, mean) ## Convert 3 studies into data.frame for simplicity's sake. study3b <- as.data.frame(study3a) study4b <- as.data.frame(study4a) study7b <- as.data.frame(study7a) diff3 <- (study3b[4, ] - study2b[4, ]) - (study3b[1, ] - study2b[1, ]) diff4 <- (study4b[4, ] - study2b[4, ]) - (study4b[1, ] - study2b[1, ]) diff7 <- (study7b[4, ] - study2b[4, ]) - (study7b[1, ] - study2b[1, ]) diff3; diff4; diff7 ## Study 2 also has the positive effects of asking by gay canvasser on ## marriage script. Overall, if gay canvasser asks about marriage, ## rate of suppor for same-sex marriage become higher.
章末問題解答(1章-1) Rコード
https://www.econ-stat-grad.com/entry/statistics/qss/solution/ch1-1www.econ-stat-grad.com
章末問題解答(1章-2) Rコード
https://www.econ-stat-grad.com/entry/statistics/qss/solution/ch1-2www.econ-stat-grad.com
ggplot2::geom_pointで人口に関する地図を書く（その1）
はじめに
この記事では、以下のような関西地域の人口に関する地図をRで作成する方法を紹介する。海外新聞ぽいフォントとかを使いお洒落に仕上げることが目標。出来るだけ丁寧に書いているつもり。全てのコードを載せている。ステップバイステップでやるのを心がけた。バリエーションがないため、レンジを区切って離散クラスにすればもっと綺麗に描画できるが、疲れたので次回以降にする。（間違いなどあればコメントください。）
1. データソース
国土数値情報のウェブサイトから「1kmメッシュ別将来推計人口（H30国政局推計）（shape形式版）」の関西二府四県のデータをダウンロードする。ファイルサイズが大きいので注意。
ダウンロードしたshapeファイルをQGIS上で開くとこんな感じになる（EPSG:6674）。
2. データ読み込み & 確認
# Setup library(ggplot2) library(dplyr) library(sf) # library(ggthemes) library(viridis) library(extrafont) # Set your wd setwd() # Load the data shiga <- read_sf("1km_mesh_suikei_2018_shape_25/1km_mesh_2018_25.shp") kyoto <- read_sf("1km_mesh_suikei_2018_shape_26/1km_mesh_2018_26.shp") osaka <- read_sf("1km_mesh_suikei_2018_shape_27/1km_mesh_2018_27.shp") hyogo <- read_sf("1km_mesh_suikei_2018_shape_28/1km_mesh_2018_28.shp") nara <- read_sf("1km_mesh_suikei_2018_shape_29/1km_mesh_2018_29.shp") wakayama <- read_sf("1km_mesh_suikei_2018_shape_30/1km_mesh_2018_30.shp")
試しに先ほどの図と同じようなものを書いてみる。
# Draw Maps kansai <- ggplot() + geom_sf(data = shiga, color = "red") + geom_sf(data = kyoto, color = "green") + geom_sf(data = osaka, color = "blue") + geom_sf(data = hyogo, color = "cyan") + geom_sf(data = nara, color = "magenta") + geom_sf(data = wakayama, color = "yellow") + theme_minimal() + theme(panel.grid = element_blank()) + labs(x = "Longitude", y = "Latitude") + theme(axis.text = element_text(size = 10), axis.title.x = element_text(size = 10), axis.title.y = element_text(size = 10)) kansai ggsave(kansai, filename = "kansai.png", width = 6, height = 8)
こんな感じになる。色は適当。後で描画するものを映えさせるためにあえてダサくした。
3. データの変換etc
上で見たように、今読み込んでいるデータは1kmメッシュのデータである。これをポイント（点）に変換したい。重心を取る。
sf::st_centroidを用いる。
# Calculate centroid of polygon using sf::st_centroid shiga$centroids <- st_transform(shiga, 6674) %>% st_centroid() %>% st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>% st_geometry() kyoto$centroids <- st_transform(kyoto, 6674) %>% st_centroid() %>% st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>% st_geometry() osaka$centroids <- st_transform(osaka, 6674) %>% st_centroid() %>% st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>% st_geometry() hyogo$centroids <- st_transform(hyogo, 6674) %>% st_centroid() %>% st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>% st_geometry() nara$centroids <- st_transform(nara, 6674) %>% st_centroid() %>% st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>% st_geometry() wakayama$centroids <- st_transform(wakayama, 6674) %>% st_centroid() %>% st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>% st_geometry()
データフレームに変換する。
# Convert to data frame shiga <- as.data.frame(shiga, xy = T) kyoto <- as.data.frame(kyoto, xy = T) osaka <- as.data.frame(osaka, xy = T) hyogo <- as.data.frame(hyogo, xy = T) nara <- as.data.frame(nara, xy = T) wakayama <- as.data.frame(wakayama, xy = T)
longitudeとlatitudeを抽出し、それぞれを独立の列に変数として加える。
# Extract and split lat long coordinates from the point # refer: https://stackoverflow.com/questions/38637070/extract-and-split-lat-long-coordinates-from-wkt-point-data-in-r shiga$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", shiga$centroids)) shiga$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", shiga$centroids)) kyoto$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", kyoto$centroids)) kyoto$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", kyoto$centroids)) osaka$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", osaka$centroids)) osaka$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", osaka$centroids)) hyogo$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", hyogo$centroids)) hyogo$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", hyogo$centroids)) nara$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", nara$centroids)) nara$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", nara$centroids)) wakayama$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", wakayama$centroids)) wakayama$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", wakayama$centroids))
必要のない変数を削除する。everything()で順番を変更する。これは好みの問題。Stataのorderコマンドと同じ。
shiga <- shiga %>% dplyr::select(long, lat, PTN_2015) %>% dplyr::select(long, lat, PTN_2015, everything()) %>% dplyr::rename(Population = PTN_2015) kyoto <- kyoto %>% dplyr::select(long, lat, PTN_2015) %>% dplyr::select(long, lat, PTN_2015, everything()) %>% dplyr::rename(Population = PTN_2015) osaka <- osaka %>% dplyr::select(long, lat, PTN_2015) %>% dplyr::select(long, lat, PTN_2015, everything()) %>% dplyr::rename(Population = PTN_2015) hyogo <- hyogo %>% dplyr::select(long, lat, PTN_2015) %>% dplyr::select(long, lat, PTN_2015, everything()) %>% dplyr::rename(Population = PTN_2015) nara <- nara %>% dplyr::select(long, lat, PTN_2015) %>% dplyr::select(long, lat, PTN_2015, everything()) %>% dplyr::rename(Population = PTN_2015) wakayama <- wakayama %>% dplyr::select(long, lat, PTN_2015) %>% dplyr::select(long, lat, PTN_2015, everything()) %>% dplyr::rename(Population = PTN_2015) # Append the prefectures' data kansai_pop <- dplyr::bind_rows( shiga, kyoto, osaka, hyogo, nara, wakayama )
4. Maps−地図を書く−
冒頭で掲載した地図を作成する。
# Draw a map (1) pop_map1 <- ggplot() + geom_point(data= kansai_pop, aes(x = long, y = lat, color = Population), alpha = 0.5) + scale_color_gradient(low = "#edd6d4", high = "#7A0018") + labs(x = NULL, y = NULL, title = "Kansai's Regional Demographics", subtitle = "Population in Kansai Municipalities, 2015", caption = "Source: National Spatial Planning and Regional Policy Bureau") + theme_void() + theme( plot.title = element_text(vjust = -2, hjust = 0.1), plot.subtitle = element_text(vjust = -2, hjust = 0.1), plot.caption = element_text(vjust = 3, hjust = 0.8), text = element_text(family = "Optima"), plot.background = element_rect(fill = "#f5f5f2"), legend.position = c(0.13, 0.18)) pop_map1 ggsave(pop_map1, filename = "pop_map1.png", width = 6, height = 8)
せっかくgeom_pointを使っているので、ポイントのサイズを人口に合わせて変えてみる。
# Draw a map (1) pop_map2 <- ggplot() + geom_point(data= kansai_pop, aes(x = long, y = lat, size = Population, color = Population), alpha = 0.3) + scale_color_gradient(low = "#edd6d4", high = "#7A0018") + # Show a legend for a color guides(size = FALSE) + scale_size_area(max_size = 15) + labs(x = NULL, y = NULL, title = "Kansai's Regional Demographics", subtitle = "Population in Kansai Municipalities, 2015", caption = "Source: National Spatial Planning and Regional Policy Bureau") + theme_void() + theme( plot.title = element_text(vjust = -2, hjust = 0.1), plot.subtitle = element_text(vjust = -2, hjust = 0.1), plot.caption = element_text(vjust = 3, hjust = 0.8), text = element_text(family = "Optima"), plot.background = element_rect(fill = "#f5f5f2"), legend.position = c(0.13, 0.18)) pop_map2 ggsave(pop_map2, filename = "pop_map2.png", width = 6, height = 8)
参考にしたサイト timogrossenbacher.ch
【つぶやき】米研究、危うい中国排除 中国は「独立」へ着々：日本経済新聞
【感想】走ることについて語るときに僕の語ること（What I Talk About When I Talk About Running）
村上春樹さんの走ることについて語るときに僕の語ること
この本は小説家でありランナーでもある村上さんのエッセイ集（村上さん本人はメモワール（回顧録）と呼んでいる）である。彼は33歳でランニングを始め、30年以上毎日平均10km走る生活を続けいてる。これまで何十回も世界各地の大会でフルマラソンを完走しており、トライアスロンや100kmウルトラマラソンも何度も完走してきている。
「小説をしっかり書くために身体能力を整え、向上させる」ことが走ることの第一目的であると言う。彼の仕事として小説を書き続ける力の源泉を、「走ること」と「彼の人生観」と絡めながら語ってくれる。話は村上さんがどのようにしてランニングを始めたのか、これまで参加してきた大会やそこでの完走までの感想、そして30年以上小説家（ランナー）として走り続けてきた経験から学んだことが書かれており、創造をする上で規則正しく健康的な生活することがどれほど重要か教えてくれる。
これ以上の説明は不要である。私の感想としては、修士課程入学前に出会っておきたい本だった。早く読んでおきたかった習慣に関する本と言えば、できる研究者の論文生産術 どうすれば「たくさん」書けるのか
前書きで引用されている、村上さんがパリで読んだ雑誌に書かれていたランナーのコメントが印象的だった。"Pain in inevitable. Suffering is optional." これはそのランナーのマントラ（呪文）であり、マラソンのレース中に自らを叱咤激励するために頭の中で唱えている言葉だそうだ。訳せば「痛みは避けがたいが、苦しみはオプショナル（こちら次第）」のようになる。「きつい」と感じることは避けようのない事実だが、「もう駄目」かどうかは自分次第ということだ。これはマラソンやランニング以外の生活で経験することにも通ずるのではないかと思う。カッコいい言葉なので、さっき紙に書いて部屋の壁に貼り付けた。
ちなみにこの本のタイトルは、レイモンド・カーヴァーの短編集のタイトル愛について語るときに我々の語ること（What We Talk About When We Talk About Love）
前回適当に書いて投稿した読書感想が以外に人気だったらしく、少し続けてみようと思った。さて、私は晩ご飯を食べてから（今日も）ランニングをしてくる。
【感想】デカルトの憂鬱 マイナスの感情を確実に乗り越える方法
今週末は論文の改訂を行い、読みかけだった本一冊を含め、二冊読み終えることができた。その内の一冊が結構よかったので簡単に紹介する。
その本はこちら、『デカルトの憂鬱 マイナスの感情を確実に乗り越える方法
この本を二言で表現するならば、「デカルトの解説書」、「デカルトから学ぶ現代の人生指南書」となるだろう。この本は21章から構成されており、全373ページと結構長い。特徴として、21の動詞が章題になっているが、「辞書」や「参考書」のように読み進めることも可能なので、あまり興味がないところは流し読みをしてもいいし、読まなくてもいい。私は、二周目以降の読書を楽しむ傾向にある。もちろん一周目はそれはそれで刺激があり楽しいが、いい作品というのは、無添加天日干しのスルメのように、噛めば噛むほど味が出る。深みに到達できる。これは音楽でも映画でも同じだと思う。
少し私の好みについて書いてしまったが、二周目以降簡単に味わい深いところを噛めるように、人生指南章をメモランダムとしてここに書いておく（既に裏表紙に書き込んでいるが）。 ◎がVery Good、◯がGoodという評価。もちろん、章全体では何を言っているか分からなくても、素敵な考えや文章は本全体に散りばめられていた。
【◎】4章『デカルトは 意外と「休む」』、14章『デカルトは 少しずつ「慰める」』、15章『デカルトは 穏やかに「暮らす」』
【○】7章『デカルトは 自分のなかを「旅する」』、8章『デカルトは ミツバチのように「本を読む」』、11章『デカルトは 魂を耕すために「学ぶ」』、17章『デカルトは 遠大に「準備する」』
全ての人が楽しめる内容ではあると思うが、知的生産を本業としている人、挑戦・挫折・リセットを繰り返している人、精神的に強くなりたい人にはかなりオススメできる本である。おそらくこのブログの読者である、研究活動を行っている皆さんにもオススメしたい。もちろん、哲学を学びたい、デカルトが好きなどの目的も十分に達成できると思う。哲学素人なので分からないけど。
さてさて、今日は天気も微妙なのでスタバに行かずに部屋に籠ることにしている。夕飯までには、少し残している論文の改訂を終えたい。次回の読書感想は、『愛について語るときに我々の語ること
【つぶやき】自分の価値を自問する：日本経済新聞
昨日、2020年9月11日の朝刊一面の記事。
【パクスなき世界】自分の価値を自問する：日本経済新聞
www.nikkei.com
最後の3段落にはこう書かれている。
企業に所属するだけで豊かさの恩恵を被ることができる時代は終わった。断絶をあおるのでなく、変化が不可逆的だと理解し、時代に合った価値観へのシフトを促すのは政治の役割だ。
突然のウイルスは特定の組織や価値観との結びつきに依存する社会のもろさを露呈させた。一方で自身の強度を高め、身近な人たちとの新たなつながりを求めるなど個のポートフォリオの組み替えも始まる。
ペストの流行が個性を花咲かせたルネサンスにつながったように、不透明な時代を主体的に生きる強靱（きょうじん）さをどう身につけるか。自分が生き抜くよりどころを増やすためにも、新たな挑戦を始めるしかない。
昨今、この第1段落に書かれているような意見をよく聞く。特にコロナウイルスのパンデミック以降、3日に1回はなんらかの記事で目にしている気がする。自分がどのように生きていくのか。1年先のことすら見えない時代だからこそ、自分が生き抜くよりどころを増やすという考えもある。自分はどうしよう。
【2020年9月7日】台風10号
今回の台風10号で被災された方々にお悔やみとお見舞いを申し上げます。
台風を含め、近年、自然災害が猛威を振るっていますね。被災状況の映像や写真を見て唖然とするしかありません。いつ自分の住む地域が選ばれるのか。いろいろ考えさせられます。
下の一つ目の記事は少し読む価値がありそうなので載っけておきます。
【名曲を紹介します】to U - Bank Band with Salyu
名曲を紹介します。今年自分で7万回くらい再生した曲です。
今年2020年。ありがとうございます。 www.youtube.com
2012年 www.youtube.com
2012年東北みちのく www.youtube.com
2010年 www.youtube.com
2009年 www.youtube.com
2008年 www.youtube.com
2007年、コブクロ優勝。 www.youtube.com
2007年ミスチル HOME TOUR www.youtube.com
おまけ
エムステ www.youtube.com
これもむっちゃ好き www.youtube.com
