-
Notifications
You must be signed in to change notification settings - Fork 0
/
LCA_profileplot_function.R
141 lines (113 loc) · 3.87 KB
/
LCA_profileplot_function.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
profile_plot = function(data){
f = with(data, cbind(tax, religion, free_election, state_aid, civil_rights, women)~1)
min_bic <- 1000000
for(i in 1:7){
lc <- poLCA(f, data, nclass=i, maxiter=3000,
tol=1e-5, na.rm=FALSE,
nrep=10, verbose=TRUE, calc.se=TRUE)
if(lc$bic < min_bic){
min_bic <- lc$bic
LCA_best_model <- lc
}
}
probs = LCA_best_model$probs
n_class = length(LCA_best_model$P)
profile_tb = data.frame(
tax = replicate(n_class, NA),
religion = replicate(n_class, NA),
free_election = replicate(n_class, NA),
state_aid = replicate(n_class, NA),
civil_rights = replicate(n_class, NA),
women = replicate(n_class, NA))
for (i in 1:6) {
if (length(probs[[i]][1,]) < 10) {
probs[[i]] = cbind(probs[[i]], matrix(0, nrow = nrow(probs[[i]]), ncol = 10 - ncol(probs[[i]])))
}
profile_tb[, i] = probs[[i]] %*% 1:10
}
rownames(profile_tb) = paste(rep("class", n_class), seq(1, n_class, 1), sep = "_")
profile_tb = rownames_to_column(profile_tb)
colnames(profile_tb)[1] = "class"
profile_long = reshape::melt(profile_tb, id.vars = "class")
p = ggplot(profile_long, aes(x = variable, y = value, group = class, color = class)) +
geom_point(size = 2.25, aes(shape = class))+
geom_line(size = 1.00) +
labs(x = NULL, y = "Mean value of the response") +
theme_bw(base_size = 14)+
ggtitle(paste(paste("class", 1:length(LCA_best_model$P), sep = "_"),
round(LCA_best_model$P, 3), collapse = ", "))+
theme(plot.margin=unit(c(1.5,1.5,1.5,1.2),"cm"))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 13))+
theme(plot.title = element_text(hjust = 0.5, size = 10))
print(p)
return(p)
# library(plotly)
#
# plotly_p = ggplotly(p, tooltip = "all") %>%
# layout(legend = list(orientation = "h", y = 1.2))
# print(plotly_p)
# return(plotly_p)
}
gen_country = function(country_name){
a = nonnegtive %>%
filter(country == country_name) %>%
dplyr::select(c("tax", "religion", "free_election", "state_aid",
"civil_rights", "women"))
return(a)
}
country_list = unique(nonnegtive$country)
for (country in country_list) {
pdf(paste(country,"LCA.pdf", sep = "_"), width = 9)
profile_plot(gen_country(country))
dev.off()
}
#######-----------Prototying code
# f = with(data, cbind(tax, religion, free_election, state_aid, civil_rights, women)~1)
# min_bic <- 100000
#
#
# for(i in 2:7){
# lc <- poLCA(f, data, nclass=i, maxiter=3000,
# tol=1e-5, na.rm=FALSE,
# nrep=10, verbose=TRUE, calc.se=TRUE)
# if(lc$bic < min_bic){
# min_bic <- lc$bic
# LCA_best_model<-lc
# }
# }
#
# probs = LCA_best_model$probs
# n_class = length(LCA_best_model$P)
#
# profile_tb = data.frame(
# tax = replicate(n_class, NA),
# religion = replicate(n_class, NA),
# free_election = replicate(n_class, NA),
# state_aid = replicate(n_class, NA),
# civil_rights = replicate(n_class, NA),
# women = replicate(n_class, NA))
#
# for (i in 1:6) {
# if (length(probs[[i]][1,]) < 10) {
# probs[[i]] = cbind(probs[[i]], matrix(0, nrow = nrow(probs[[i]]), ncol = 10 - ncol(probs[[i]])))
# }
# profile_tb[, i] = probs[[i]] %*% 1:10
# }
#
#
# rownames(profile_tb) = paste(rep("class", n_class), seq(1, n_class, 1), sep = "_")
# profile_tb = rownames_to_column(profile_tb)
# colnames(profile_tb)[1] = "class"
# profile_long = reshape::melt(profile_tb, id.vars = "class")
#
#
# p = ggplot(profile_long, aes(x = variable, y = value, group = class, color = class)) +
# geom_point(size = 2.25)+
# geom_line(size = 1.25) +
# labs(x = NULL, y = "Mean value of the response", main = "Profile plot") +
# theme_bw(base_size = 14)
#
# library(plotly)
#
# ggplotly(p, tooltip = "all") %>%
# layout(legend = list(orientation = "h", y = 1.2))