裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

Julia で,散布図と回帰直線の信頼限界帯を描画--その 2

2022年03月14日 | Julia

https://data-viz-lab.com/ggplot2

であるが,色々あって,最終的には

が成果物として挙げられている。

所要のプログラムは 51 行とのこと(キータッチ数は結構あるなあ)。
結構細かい指定が一杯ある。
とても,「簡単にきれいな図が得られます」とは言えないような。

# First draw a basic plot with main arguments
plot <- ggplot(data = data_iris,             #set the data
               mapping = aes(x=Sepal.Length, #set x axis
                             y=Petal.Length, #set y axis
                             colour=Species) #set color by species
) +
  geom_point() +              #set a scatter plot
  geom_smooth(method = "lm")  #set a linear line

# Then make it nicer
plot <- plot + 
  theme(axis.line = element_line(colour = "black",  #set color of axis line
                                 size = 0.7,        #size
                                 linetype = "solid"), #and line type
        legend.justification = c(1, 1), #set legend inside plot
        legend.position= c(1, 0.5),     #set position by x and y axis
        legend.title = element_text(size = 12, color = "black", 
                                    face = "bold"),
        legend.text = element_text(size = 10, color = "black", 
                                   face = c("italic")),
        legend.background = element_blank(),
        legend.key = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank()
  ) +
  labs(title="Association between the length of petal and setal", 
       subtitle="From IRIS dataset") +
  scale_x_continuous(name = "Lenght of Sepal (cm)", 
                     breaks = seq(4,8,1),
                     limits = c(4,8) 
  ) +
  scale_y_continuous(name = "Lenght of Petal (cm)",
                     breaks = seq(0, 8, 1),
                     labels = seq(0, 8, 1),
                     limits = c(0, 8)) +
  theme(axis.title.x = element_text(size = 12, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y=element_text(size=12, vjust = 0.8),
        plot.title = element_text(size = 12,  #font size and adjust
                                  hjust = 0.5,#adjust
                                  face = "bold"),
        plot.subtitle = element_text(size = 10,  #font size and adjust
                                     hjust = 0.5,#adjust
                                     face = "bold")
  ) +
  scale_colour_viridis(discrete = TRUE, alpha=0.6)

plot

Julia では,前回のプログラムを少し修正して,複数のデータフレームでの描画を蓄積して最終結果を表示するという戦略を立てた。

まずは,scatterplot() のバージョンアップ版

function scatterplot(df; xlabel="", ylabel="", grouplabel="", col=:black, conflevel=0.95,
                     add=false) # append
    gr( tick_direction=:out,
        grid=false,
        markerstrokewidth=0,
        alpha=0.4,
        fontfamily="serif",
        guidefontfamily="times",
        guidefontsize=8,
        titlefontfamily="times", # append
        titlefontsize=10,        # append
        label="")
    df2 = DataFrame(x=df[:,1], y=df[:,2])
    model = lm(@formula(y ~ x), df2);
    minx, maxx = extrema(df2.x)
    pred = DataFrame(x=range(minx, maxx, length=500));
    pr = predict(model, pred, interval=:confidence, level=conflevel);
    add || (p = plot(xlabel=xlabel, ylabel=ylabel)) # append
    p = scatter!(df2.x,
                df2.y,
                color=col);#,
                #xlabel=xlabel, # delete
                #yaxis=ylabel); # delete
    plot!(      pred.x,
                pr.prediction,
                linecolor=:red,
                linewidth=2,
                color=:gray,
                ribbon=(pr.prediction .- pr.lower,
                        pr.upper .- pr.prediction))
    for gl in grouplabel
        annotate!(gl[1],
                gl[2],
                text(gl[3],
                "times",
                8,
                gl[4]))
    end
    p
end

これを使うにあたっての事前準備。まあ,単にデータフレームの分割であるが。

gd = groupby(iris, :Species);

続いては,分割されたデータフレームについて scatterplot() を呼び出す。
ただし,2 つ目以降は,グラフの追加(add=true)。

for i in 1:length(gd)
    scatterplot(gd[i][:,[:SepalLength, :PetalLength]],
        xlabel="Length of Sepal(cm)",
        ylabel="Length of Petal(cm)",
        grouplabel=[[(5.5, 2.2, "Setosa", :red),
          (6.5, 3.8, "Versicolor", :green),
          (6.0, 6.4, "Virginica", :blue)][i]],
        col=[:red, :green, :blue][i],
        add=i != 1)
end

最後に,タイトルを付けておきましょう。

title!("Association between the length of petal and sepal\nFrom IRIS dataset")

出来上がりは,以下の通り。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村