Чтение данных
Я скачал данные NES2012 здесь , затем прочитал NES2012.sav это с foreign pacakge.
library(foreign)
data <- read.spss(file = "NES2012.sav",
to.data.frame = T,
use.value.labels = T)
Убедитесь, чтоиспользуйте метки значений и считайте их как data.frame
Таблица сопряженности
t1<-table(data$gender,data$voted2012)
print(t1)
Вывод на консоль:
Did not vote Voted
-1. Inapplicable 0 0
Male 537 2137
Female 572 2267
У нас есть уровень для пропущенных значений впеременная пола мы можем подмножество / перекодировать набор данных, но, в этом случае, мы можем просто подмножество таблицы
t1 <- t1[2:3,] # pick the second to the third row and al the columns
print(t1)
Консольный вывод:
Did not vote Voted
Male 537 2137
Female 572 2267
Теперь процент, мы используем проп.таблицу и умножить, затем округлить
pt <- round(prop.table(t1) * 100,2)
Вы можете добавить поля с помощью apply, cumsum или rowum, но addmargins - самый простой.
ptm <- addmargins(pt)
Теперь измените имена
colnames(ptm)[3] <- "Total" # only the third one
rownames(ptm)[3] <- "Total" # same thing
print(ptm)
Вывод на консоль:
Did not vote Voted Total
Male 9.74 38.76 48.50
Female 10.38 41.12 51.50
Total 20.12 79.88 100.00
Среднее сравнение
Я предполагаю, что тест, который вы используете, является ANOVA, из-за уровня вопросов, которые мне кажутсякак вводный курс статистики или количественных методов в SSCC.Но, не зная шкалу измерения «Рейтинга термометров федерального правительства», я бы не рекомендовал этот метод.Просто для выполнения задания давайте идти с ним!С ANOVA мы работаем под нулевым гипотезом (H0 :), что нет никакой разницы между средними значениями.
Сначала давайте возьмем простую и простую таблицу с tapply, для этого мы должны подготовить очистку и перекодировать переменные.Ftgr_fedgov должен быть числовым, а pid_3 - множителем с 3 уровнями.Я перекодирую уровни, которые не являются числовыми "1" ... "99" (пропущенные значения и т. Д.), В NA, которые в моем случае являются первыми 8 уровнями, возвращаемыми функцией levels().Обратите внимание, что у меня есть коэффициент, потому что я указал, что при чтении файла spps .sav в первых строках.
levels(data$ftgr_fedgov)[1:8] <- NA
принудительно преобразовать в числовые значения
data$ftgt_fedgov <- as.numeric(data$ftgr_fedgov)
использовать анонимную функцию счтобы избавиться от NA с помощью na.rm = T
the.means <-tapply(data$ftgt_fedgov, data$pid_3, function(x) mean(x,na.rm = T))
meantab <- cbind(the.means)
colnames(meantab) <- "Mean Federal Gov rating"
print(meantab)
Вывод на консоль:
Mean Federal Gov rating
Dem 32.81931
Ind 23.73232
Rep 19.51978
Теперь давайте перейдем к ANOVA:
myanova <-anova(lm(data = data, ftgt_fedgov ~ pid_3))
print(myanova)
Вывод на консоль:
Analysis of Variance Table
Response: ftgt_fedgov
Df Sum Sq Mean Sq F value Pr(>F)
pid_3 2 164624 82312 342.58 < 2.2e-16 ***
Residuals 5426 1303697 240
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Итак, у нас достаточно доказательств, чтобы отвергнуть ноль, теперь мы можем посмотреть, какая группа отличается, есть много способов сделать это, самый распространенный - Тьюки.
myanova <- aov(data$ftgt_fedgov ~ data$pid_3)
TukeyHSD(myanova)
Вывод на консоль:
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = data$ftgt_fedgov ~ data$pid_3)
$`data$pid_3`
diff lwr upr p adj
Ind-Dem -9.086995 -10.217986 -7.956004 0
Rep-Dem -13.299528 -14.576857 -12.022198 0
Rep-Ind -4.212533 -5.515013 -2.910053 0
В этом случае все три группы означают разные (p прил.)
, вы можете увидеть это с помощью:
par(mfcol = c(1,2))
boxplot(data$ftgt_fedgov ~ data$pid_3, main = "Federal Government Thermometer Rating",
col = c("lightblue","lightgoldenrodyellow","pink"))
plot(TukeyHSD(myanova))
Кроме того, вам нужно будет проверить модель с учетом невязок, нормальности переменных подмножеств и т. Д., Но мне кажется, что для назначения это не требуется.
Дополнительное задание
Моя гипотеза состоит в том, что люди, которые посещают более частоРелигиозные службы наименее благоприятны для легализации марихуаны.
DV = pot_legal3
IV = религия_attend
Обе переменные являются факторами, поэтому я использую chisq.test() со своей таблицей.
t1 <- table(data$relig_attend,data$pot_legal3)[,1:3] # removed weird column with [,1:3]
ht <- chisq.test()
ht
Вывод на консоль:
Pearson s Chi-squared test
data: table(data$relig_attend, data$pot_legal3)[, 1:3]
X-squared = 470.58, df = 8, p-value < 2.2e-16
Значение p меньше 0,001, у нас достаточно доказательств, чтобы отклонить нулевой гипотезы независимости.
Теперь давайте проверимстандартизированные остатки для проверки направления зависимости.с участком
library(ggplot2)
# chi & table into data.frame
df1 <- as.data.frame(t1)
df1$res <- as.vector(round(ht$stdres,2))
df1$exp <- as.vector(round(ht$expected,0))
df1$obv <- as.vector(ht$observed)
ggplot(df1,aes(Var1,Var2,fill = res, label = paste(df1$obv," / ",df1$exp,"\n",df1$res,sep = "")))+
geom_raster()+
scale_fill_gradient2(low = "blue", mid = "white",high = "red")+
geom_text()+
guides( fill = guide_colorbar(title.position = "top"))+
scale_x_discrete(position = "top")+
theme_minimal()+
theme(axis.text.y = element_text(angle = 90, size = 14, hjust = 0.5),
axis.text.x = element_text(size = 12),
plot.caption = element_text(face = "italic"),
legend.position = "bottom",
legend.direction = "horizontal",
legend.key.height = unit(1,"mm"),
legend.key.width = unit(14,"mm"),
legend.title = element_text(size = 10),
legend.title.align = 0.5,
plot.title = element_text(hjust = 0.5,face = "bold.italic", size = 20, colour = "darkseagreen" )) +
labs(title ="Pot & God",
y = "Marihuana legalization" ,
x = "Religious service attendance",
fill = "Standarized residuals",
caption = "Data: SAGE NES2012 N: 5440",
subtitle = "Observer / expected")
