Hi Colin, I have never heard of a binomial distribution z statistic with (or without for that matter) a continuity correction, but I am not a statistician. Other's may have some ideas there. As for other ways to analyze the data, I skimmed through the article and brought the data and played around with some different analyses and graphs. I attached a file headache.txt with all the R script (including the data in an Rish format). It is really a script file (i.e., .R) but for the listservs sake I saved it as a txt. There are quite a few different things I tried in there so hopefully it gives you some ideas. Regardless of the analysis type used and whether one considers proportion that "significantly improved" or the raw frequency or intensity scores, I would say that concluding the treatment was effective is a good conclusion. The only real concern could be that people would naturally get better on their own (a control group would be needed to bolster the causal inference drawn from a pre/post measurement). However, given at least what I know about migraines, it is often a fairly chronic condition so over a relatively short time period, it seems implausible to conclude that as many people would be improving as this study reported. Cheers, Josh
On Sat, Nov 19, 2011 at 7:43 PM, Colstat <colstat at gmail.com> wrote:
hey, Joshua I was reading this paper, in attachment, and reproducing the results. I was really confused when he said in the paper "The results were then statistically analyzed using binomial distribution z statistics with continuity correction."? The data is binomial?? To me, this is a paired t-test. What command should I use to get those results (the first paragraph in Results section)?? Basically, it's a pre and post treatment problem. What other graphical analysis do you think is appropriate? reshape package? lattice package, namely conditional graph? I know this might be too much, but I do really appreciate it if you do take a look at it. Thanks, Colin On Sat, Nov 19, 2011 at 10:15 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
Hi, I am not clear what your goal is. ?There is a variety of data there. You could look at t-test differences in preIntensity broken down by sex, you could use regression looking at postIntensity controlling for preIntensity and explained by age, you could.... Why are you analyzing data from an article? ?What did the article do? What you mention---some sort of z statistic (what exactly this was of and how it should be calculated did not seem like was clear even to you), histogram, t-test, lattice, are all very different things that help answer different questions, show different things, and in one is a piece of software. Without a clearer question and goal, my best advice is here are a number of different functions some of which may be useful to you: ls(pos = "package:stats") Cheers, Josh On Sat, Nov 19, 2011 at 3:01 PM, Colstat <colstat at gmail.com> wrote:
Dear R experts, I am trying to analyze data from an article, the data looks like this Patient Age Sex Aura preCSM preFreq preIntensity postFreq postIntensity postOutcome 1 47 F A 4 6 9 2 8 SD 2 40 F A/N 5 8 9 0 0 E 3 49 M N 5 8 9 2 6 SD 4 40 F A 5 3 10 0 0 E 5 42 F N 5 4 9 0 0 E 6 35 F N 5 8 9 12 7 NR 7 38 F A 5 NA 10 2 9 SD 8 44 M A 4 4 10 0 0 E 9 47 M A 4 5 8 2 7 SD 10 53 F A 5 3 10 0 0 E 11 41 F N 5 6 7 0 0 E 12 49 F A 4 6 8 0 0 E 13 48 F A 5 4 8 0 0 E 14 63 M N 4 6 9 15 9 NR 15 58 M N 5 9 7 2 8 SD 16 53 F A 4 3 9 0 0 E 17 47 F N 5 4 8 1 4 SD 18 34 F A NA ?5 9 0 0 E 19 53 F N 5 4 9 5 7 NR 20 45 F N 5 5 8 5 4 SD 21 30 F A 5 3 8 0 0 E 22 29 F A 4 5 9 0 0 E 23 49 F N 5 9 10 0 0 E 24 24 F A 5 5 9 0 0 E 25 63 F N 4 19 7 10 7 NR 26 62 F A 5 8 9 11 9 NR 27 44 F A 5 3 10 0 0 E 28 38 F N 4 8 10 1 3 SD 29 38 F N 5 3 10 0 0 E How do I do a binomial distribution z statistics with continuity correction? basically normal approximation. Could anyone give me some suggestions what I (or R) can do with these data? I have tried tried histogram, maybe t-test? or even lattice? ?what else can I(or can R) do? help please, thanks so much. ? ? ? ?[[alternative HTML version deleted]]
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
-- Joshua Wiley Ph.D. Student, Health Psychology Programmer Analyst II, ATS Statistical Consulting Group University of California, Los Angeles https://joshuawiley.com/
Joshua Wiley Ph.D. Student, Health Psychology Programmer Analyst II, ATS Statistical Consulting Group University of California, Los Angeles https://joshuawiley.com/ -------------- next part -------------- #### Data (provided via dput) dat <- structure(list(patient = 1:29, age = c(47L, 40L, 49L, 40L, 42L, 35L, 38L, 44L, 47L, 53L, 41L, 49L, 48L, 63L, 58L, 53L, 47L, 34L, 53L, 45L, 30L, 29L, 49L, 24L, 63L, 62L, 44L, 38L, 38L), sex = structure(c(1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Female", "Male"), class = "factor"), nonaura = structure(c(1L, 2L, 3L, 1L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 3L, 3L, 1L, 3L, 1L, 3L, 3L, 1L, 1L, 3L, 1L, 3L, 1L, 1L, 3L, 3L), .Label = c("A", "A/N", "N"), class = "factor"), csm = c(4L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 4L, 5L, 5L, 4L, 5L, 4L, 5L, 4L, 5L, NA, 5L, 5L, 5L, 4L, 5L, 5L, 4L, 5L, 5L, 4L, 5L), freq0 = c(6L, 8L, 8L, 3L, 4L, 8L, NA, 4L, 5L, 3L, 6L, 6L, 4L, 6L, 9L, 3L, 4L, 5L, 4L, 5L, 3L, 5L, 9L, 5L, 19L, 8L, 3L, 8L, 3L), inten0 = c(9L, 9L, 9L, 10L, 9L, 9L, 10L, 10L, 8L, 10L, 7L, 8L, 8L, 9L, 7L, 9L, 8L, 9L, 9L, 8L, 8L, 9L, 10L, 9L, 7L, 9L, 10L, 10L, 10L), freq1 = c(2L, 0L, 2L, 0L, 0L, 12L, 2L, 0L, 2L, 0L, 0L, 0L, 0L, 15L, 2L, 0L, 1L, 0L, 5L, 5L, 0L, 0L, 0L, 0L, 10L, 11L, 0L, 1L, 0L), inten1 = c(8L, 0L, 6L, 0L, 0L, 7L, 9L, 0L, 7L, 0L, 0L, 0L, 0L, 9L, 8L, 0L, 4L, 0L, 7L, 4L, 0L, 0L, 0L, 0L, 7L, 9L, 0L, 3L, 0L), btxa = structure(c(2L, 1L, 2L, 1L, 1L, 3L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 2L, 1L, 2L, 1L, 3L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 2L, 1L), .Label = c("E", "SD", "NR"), class = "factor")), .Names = c("patient", "age", "sex", "nonaura", "csm", "freq0", "inten0", "freq1", "inten1", "btxa"), row.names = c(NA, -29L), class = "data.frame") table(dat$btxa) ## note p = .1 is 10% (much higher than the 1% they said they used) ## eliminated + significantly reduced binom.test(x = 16 + 8, n = 16 + 8 + 5, p = .1) ## eliminated only binom.test(x = 16, n = 16 + 8 + 5, p = .1) ## unclear what they compare the 8 SDs against ## perhaps *just* the no responses?? binom.test(x = 8, n = 8 + 5, p = .1) ## even against everyone (terribly unfair) ## do not approach p < .04 binom.test(x = 8, n = 16 + 8 + 5, p = .1) #### Examination of baseline characteristics t.test(freq0 ~ sex, data = dat) t.test(inten0 ~ sex, data = dat) t.test(csm ~ sex, data = dat) t.test(age ~ sex, data = dat) ## note how you treat the single A/N case ## will affect these results t.test(freq0 ~ I(nonaura == "N"), data = dat) t.test(inten0 ~ I(nonaura == "N"), data = dat) t.test(csm ~ I(nonaura == "N"), data = dat) t.test(age ~ I(nonaura == "N"), data = dat) ###### Tests of pre post frequency ###### ## (paired) wilcox test for frequency (with and without correction) with(dat, wilcox.test(x = freq0, y = freq1, paired = TRUE, conf.int = TRUE, correct = TRUE, exact = FALSE)) with(dat, wilcox.test(x = freq0, y = freq1, paired = TRUE, conf.int = TRUE, correct = FALSE, exact = FALSE)) ## paired t-test with(dat, t.test(x = freq0, y = freq1, paired = TRUE)) ## see if age predicts post controlling for pre ## center freq0 and age so intercept is the expected ## post frequency (freq1) for an individual with average ## pre frequency and age dat$cfreq0 <- with(dat, freq0 - mean(freq0, na.rm = TRUE)) dat$cage <- with(dat, age - mean(age, na.rm = TRUE)) summary(lm(freq1 ~ cfreq0 + cage, data = dat)) ## age appears to be related to frequency of headaches ###### Tests of pre post intensity ###### ## wilcox test for intensity (with and without correction) with(dat, wilcox.test(x = inten0, y = inten1, paired = TRUE, conf.int = TRUE, correct = TRUE, exact = FALSE)) with(dat, wilcox.test(x = inten0, y = inten1, paired = TRUE, conf.int = TRUE, correct = FALSE, exact = FALSE)) with(dat, t.test(x = inten0, y = inten1, paired = TRUE)) ## see if age predicts post controlling for pre ## center inten0 so intercept is the expected ## post intensity (inten) for an individual with average ## pre intensity and age dat$cinten0 <- with(dat, inten0 - mean(inten0, na.rm = TRUE)) summary(lm(inten1 ~ cinten0 + cage, data = dat)) ## age appears to be related to intensity of headaches ## see if pre frequency predicts post intensity ## or if pre intensity predicts post frequency summary(lm(inten1 ~ cinten0 + cage + cfreq0, data = dat)) summary(lm(freq1 ~ cinten0 + cage + cfreq0, data = dat)) ########### Some data visualization ################# require(reshape) require(ggplot2) ## categorical variables catvars <- c("patient", "sex", "nonaura", "btxa") ## continuous variables (ordered, and ignoring the centered ones we made) ## (ordered so the plots put freq0/1 and inten0/1 next to each other contvars <- c("age", "csm", "freq0", "freq1", "inten0", "inten1") ## scatter plot matrix of continuous variables plot(dat[contvars]) ## correlation matrix cor(dat[contvars], use = "pairwise.complete.obs") mdat <- melt(dat[, c(catvars, contvars)], id = c("patient", "sex", "nonaura", "btxa")) ggplot(data = mdat, aes(x = variable, y = value)) + geom_boxplot() + geom_jitter() ## age overwhelms things a bit ## we could create standardized variables ## some care is required because we want to be able to compare ## how frequency and intensity change, so they must be standardized together ## of course we could also make it a long data set with a time (pre/post) indicator ## but that makes it trickier to plot the nonchanging variables ## get mean/sd for overall freq and inten desc <- with(dat, c(mf = mean(c(freq0, freq1), na.rm = TRUE), sdf = sd(c(freq0, freq1), na.rm = TRUE), mi = mean(c(inten0, inten1)), sdi = sd(c(inten0, inten1)))) ## standardize stdat <- within(dat, { age <- scale(age) csm <- scale(csm) freq0 <- (freq0 - desc["mf"])/desc["sdf"] freq1 <- (freq1 - desc["mf"])/desc["sdf"] inten0 <- (inten0 - desc["mi"])/desc["sdi"] inten1 <- (inten1 - desc["mi"])/desc["sdi"] }) ## melt data for plotting stdmdat <- melt(stdat[, c(catvars, contvars)], id = c("patient", "sex", "nonaura", "btxa")) ## colour jittered points based on aura ggplot(data = stdmdat, aes(x = variable, y = value)) + geom_boxplot() + geom_jitter(aes(colour = nonaura)) ## it looks a bit like highest post frequencies are ## over represented in nonauras summary(lm(freq1 ~ I(nonaura == "N"), data = dat)) summary(lm(freq1 ~ I(nonaura == "N") + freq0, data = dat)) ## nonsignificant after account for pre frequency ## ## create long data just for freq tmp <- melt(dat[, c("patient", "nonaura", "freq0", "freq1")], id = c("patient", "nonaura")) ## fit linear mixed effects model ## random intercept by patient ## value is freq0 and freq1 ## variable indicates freq1 or freq0 ## nonaura seems to predict overall frequency ## the negative effect of variable indicates that ## frequency decreases from freq0 to freq1 require(lme4) summary(lmer(value ~ variable + I(nonaura=="N") + (1|patient), data = tmp))