2
$\begingroup$

So I have data from a randomized blind trial of 1mg of nicotine gum on dual n-back working memory scores; I analyzed them as usual with a t-test and found a small increase in means but a large increase in standard deviations on a f-test! Strange. I also have data for each day on mood/productivity that day on a 1-5 scale.

I wondered: is nicotine following an inverse U-curve, where it causes higher scores on the worser days (1-3) and lower scores on the better days (3-5)? I look around and it seems I want a multinomial logistic regression comparing the placebo & active days.

I enter the data & load mlogit:

nicotine <- read.table(stdin(),header=TRUE)
day      active mp score
20120824 1      3  35.2
20120827 0      5  37.2
20120828 0      3  37.6
20120830 1      3  37.75
20120831 1      2  37.75
20120902 0      2  36.0
20120905 0      5  36.0
20120906 1      5  37.25
20120910 0      5  49.2
20120911 1      3  36.8
20120912 0      3  44.6
20120913 0      5  38.4
20120915 0      5  43.8
20120916 0      2  39.6
20120918 0      3  49.6
20120919 0      4  38.4
20120923 0      5  36.2
20120924 0      5  45.4
20120925 1      3  43.8
20120926 0      4  36.4
20120929 1      3  43.8
20120930 1      3  36.0
20121001 1      3  46.0
20121002 0      4  45.0
20121008 0      2  34.6
20121009 1      3  45.2
20121012 0      5  37.8
20121013 0      4  37.2
20121016 0      4  40.2
20121020 1      3  39.0
20121021 0      3  41.2
20121022 0      3  42.2
20121024 0      5  40.4
20121029 1      2  41.4
20121031 1      3  38.4
20121101 1      5  43.8
20121102 0      3  48.2
20121103 1      5  40.6

library(mlogit)
Nicotine <- mlogit.data(nicotine,shape="wide", choice="mp")
mlogit(score ~ (active + mp)^2, Nicotine)
Error in solve.default(H, g[!fixed]) : 
  Lapack routine dgesv: system is exactly singular
Calls: mlogit ... mlogit.optim -> as.vector -> solve -> solve.default

The error also happens even with the simplest call I can think of:

mlogit(score ~ active, Nicotine)
Error in solve.default(H, g[!fixed]) : 
  Lapack routine dgesv: system is exactly singular
Calls: mlogit ... mlogit.optim -> as.vector -> solve -> solve.default

Reading the documentation for mlogit didn't much help, and look at the other questions having the same error, they're different enough I can't tell whether they apply or not.

Thank you for your assistance.

$\endgroup$
  • 1
    $\begingroup$ Take a look at table(nicotine$mp, nicotine$active). It appears you have separation: mp=4 is always associated with active=0. This is likely the problem. $\endgroup$ – Jason Morgan Nov 25 '12 at 3:05
  • 2
    $\begingroup$ Maybe I'm missing something, but could you explain why you choose multinomial regression to model the dependent variable "memory score"? To me, it seems "memory score" is not categorical but can be measured at a fine scale. $\endgroup$ – caracal Nov 25 '12 at 9:56
  • $\begingroup$ Your variables are alternative-specific, meaning that they do not vary by alternative (what your alternatives are, I cannot tell). For this model, you need to call it with mlogit(score ~ 1 | active, Nicotine). The call is mlogit(Depvar ~ generic | altspec) $\endgroup$ – gregmacfarlane Nov 25 '12 at 13:11
  • $\begingroup$ @JasonMorgan Yes, that's some separation there alright! When I collapse the 4/5 days into just a 4, I can look at stuff like summary(mlogit(mp ~ active, Nicotine)) which produce no relationship as expected. $\endgroup$ – gwern Nov 26 '12 at 0:49
1
$\begingroup$

You don't want multinomial logit as your dependent variable is a score that is nearly continuous. I would start by plotting the data e.g. with

with(nicotine, stripchart(jitter(active)~score, vertical = TRUE))

which doesn't reveal any obvious pattern.

Then you could look at a linear model:

m1 <- with(nicotine, lm(score~as.factor(active)))
summary(m1)

(this is equivalent to the t-test you ran) which shows a small and nonsignificant difference. Plotting m1 doesn't reveal anything particularly interesting to my eyes, either. You say you found large differences in variances but

with(nicotine, sd(score[active == 1]))
with(nicotine, sd(score[active == 0]))

shows the difference to be not all that large (and the stripchart shown above agrees).

Then you could add mp to the model:

m2 <- with(nicotine, lm(score~as.factor(active) + mp))
summary(m2)

which also shows only very small differences and a miniscule $R^2$

There's probably other things you could do, but it looks like there is not much to find here.

$\endgroup$
  • $\begingroup$ Thanks for your thorough look! I agree, it looks like a but. I'll keep those tricks in mind next time. $\endgroup$ – gwern Nov 26 '12 at 0:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.