Reproduce the results of analysis reported in the mental impairment and SES example by adapting the R script in the lecture note and from the rating experiment example.
The mental health status of a sample of 1,660 young New York residents in midtown Manhattan was classified by their parents’ socioeconomic status (A=High, F=Low).
Source: Srole, L. & Fisher, A.K. (1978). Mental health in the metropolis. The midtown Manhattan study. Reported in Agresti, A. (1989). Tutorial on modeling ordered categorical response data. Psychological Bulletin, 105, 290-301.
Data
Column 1: Subject ID
Column 2: SES, numerical code A=6, … , F=1
Column 3: SES, reference code for A
Column 4: SES, reference code for B
Column 5: SES, reference code for C
Column 6: SES, reference code for D
Column 7: SES, reference code for E
Column 9: Mental impairment, well=1, mild=2, moderate=3,impaired=4
#setwd("C:/Users/Ching-Fang Wu/Desktop/glm")
require(pacman)## Loading required package: pacmanpacman::p_load(tidyverse, MASS, HH)# load data
dta <- read.table("mentalSES.txt", h = T)head(dta,8)## Id ses A B C D E mental
## 1 1 6 1 0 0 0 0 1
## 2 2 6 1 0 0 0 0 1
## 3 3 6 1 0 0 0 0 1
## 4 4 6 1 0 0 0 0 1
## 5 5 6 1 0 0 0 0 1
## 6 6 6 1 0 0 0 0 1
## 7 7 6 1 0 0 0 0 1
## 8 8 6 1 0 0 0 0 1str(dta)## 'data.frame': 1660 obs. of 8 variables:
## $ Id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ ses : int 6 6 6 6 6 6 6 6 6 6 ...
## $ A : int 1 1 1 1 1 1 1 1 1 1 ...
## $ B : int 0 0 0 0 0 0 0 0 0 0 ...
## $ C : int 0 0 0 0 0 0 0 0 0 0 ...
## $ D : int 0 0 0 0 0 0 0 0 0 0 ...
## $ E : int 0 0 0 0 0 0 0 0 0 0 ...
## $ mental: int 1 1 1 1 1 1 1 1 1 1 ...Id、ses和mental的資料屬性為數值，應該改為類別變數。
#轉換資料屬性並新增變數
dta$fId<-factor(dta$Id)
dta$fses<-factor(dta$ses)
dta$fmental<-factor(dta$mental)#修改factor內部名稱
dta$fses<-factor(dta$fses, levels = c(1,2,3,4,5,6), labels = c("F","E","D","C","B","A"))
dta$fmental<-factor(dta$mental, levels = c(1,2,3,4),labels = c("well","mild","moderate","impaired"))with(dta, tapply(fId, list(fses, fmental), FUN = function(x) length(unique(x))))## well mild moderate impaired
## F 21 71 54 71
## E 36 97 54 78
## D 72 141 77 94
## C 57 105 65 60
## B 57 94 54 40
## A 64 94 58 46左欄ABCDEF的順序反了，該如何反轉?
a<-as.data.frame(with(dta, tapply(fId, list(fses, fmental), FUN = function(x) length(unique(x)))))#新增Total欄位，並加總數值
a$Total<- a$well + a$mild + a$moderate +a$impairedknitr::kable(apply(with(a,table(dta$fses,dta$fmental)), 2, cumsum))| well | mild | moderate | impaired | |
|---|---|---|---|---|
| F | 21 | 71 | 54 | 71 |
| E | 57 | 168 | 108 | 149 |
| D | 129 | 309 | 185 | 243 |
| C | 186 | 414 | 250 | 303 |
| B | 243 | 508 | 304 | 343 |
| A | 307 | 602 | 362 | 389 |
覺得這裡應該可用accumulate，但是還沒想出來怎麼做….
##likert scale plotting
#m<- as.data.frame(matrix(dta$Count, 4, 2))
#a<-as.data.frame(with(dta,ftable(dta$sesf,dta$mentalf)))#names(a) <- c("well", "Mild","Moderate","Impaired")#rownames(a) <- c("F", "E", "D", "C","B","A")#print(m)#likert(t(m), as.percent = FALSE, main = "", ylab = "Stimulus")
#likert(t(m), as.percent = TRUE, main = "", ylab = "Stimulus")#
# reorder the levels of mother height to S, M, T (rather than M S T)
#dta$Rating <- factor(dta$Rating, levels(dta$Rating)[c(1,2,4,3)])#summary(m0 <- polr(Rating ~ Stimulus, data = dta, weight = Count, method = "probit", Hess = TRUE))#confint(m0)##
# probabilities for each category
#1 - pnorm(2.0138 - 1.009)
#pnorm(2.0138 - 1.009) - pnorm(1.0104 - 1.009)
#pnorm(1.0104-1.009) - pnorm(-0.9848 - 1.009)
#pnorm(-0.9848-1.009)#x <- seq(-4, 4, by = .05)
#plot(x, pnorm(x), type = "n",
# xlab = "Perceptual dimension", ylab = "Probability")
#lines(x, pnorm(x), lty = 2, col = "steelblue")
#lines(x, pnorm(x - 1.009), lty = 2, col = "blue")
#grid()
#points(c(-0.9848, 1.0104, 2.0138), c(0.16, 0.845, 0.98),
# col = "magenta", pch = 16)
#points(c(-0.9848, 1.0104, 2.0138), c(0.025, 0.5, 0.84),
# col = "cyan", pch = 16)
#text(1.4, 0.3, "Noise+Signal", cex = 0.7)
#text(0, 0.7, "Noise", cex = 0.7)