GLMER graph info

gistfile1.txt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Generalized Linear Models
with binary categorical data
with binomial error structure
and logit link function

Categorical data is heteroskedastic. Residuals are not normally distributed. It violates assumptions of linear models. It can be approximated as a logistic curve (s-like-curve). This curve represents probability values: probability of observing 1 or 0 at y when steps in x change. (Ideal for perception data.)

model = glm(response ~ step in continuum or Hz, data=study, family="binomial") # regular generalized (logistic) linear model

Estimates are log-odds. If estimate of contrast is negative and estimate of intercept is positive with an increase in factor you have a decrease in odds. You need to calculate the inverse of the logit transform to obtain a probability (0-1).

logit.inv(estimate of intercept)
logit.inv(estimate of intercept + estimate of contrast * 0.3)
"if you are interested in 0.3 (30% of x)"
logit.inv = function(x){exp(x) / 1+exp(x)}

Generalized Linear Mixed-Effects models with LME4

glmer(y ~ x + z + (1|subject), family=binomial)

For "categorical perception" studies

lmer(response ~ step + group + (1+step|group), family = binomial)
unique(subject) # ¿Cuántos niveles tiene "subject"?
unique(step) # ¿Cuáles son los márgenes o rangos de "step"?

model = glmer(response ~ 1 + condition + (1+condition|subject), data=mydata, family="binomial")
Estimate the intercept and the slope of "condition"

model.null = glmer(response ~ 1 + (1+condition|subject), data=mydata, family="binomial")
This is a simpler model, a model in which there is no slope for "condition". We test the simple model (without slopes for "condition") against the model in which it was specified.

anova(model.null, model)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Plotting logit functions

response = as.numeric(as.character(response)) # backtransform the response into numbers

# windows(width=8, height=6)

quartz(width=8,height=6)
plot(unique(Hz), rep(c(0,1),5), type="n" #fit the margins of the plot, do not plot the points
        xaxt = "n", yaxt="n", xlab="", ylab="")

axis(side=1, at=unique(Hz), labels=unique(Hz))  # side refers to the axis of the square
                                                                                                # side=1 is bottom and then clock-wise
axis(side=2, at=c(0,1), labels=c("ba","pa"), las=2)     # side=2 is left side, add the axis information
                                                                                                        # las=2 labels are pependicular to axis
mtext(side=1, "F1 (Hz) in continuum", line=2.5, cex=1.5) # cex=1.5 makes fontsize bigger
points(jitter(F1_Hz, factor=0.2),jitter(response))      # plot x and y columns in Table
                                                                                                        # Add random jitter for visibility, otherwise everything is plotted on top of other stuff

xpoints = seq(0,360, by=0.001) # generate a sequence of points from 0 to 360 in small steps
fixef(model) # fixef() is to get the coefficients of the GLMER model
ypoints = intercept estimate + contrast estimate * xpoints # o, lo que es lo mismo...
ypoints = fixef(model)[1] + fixef(model)[2] * xpoints
# BUT ypoints are NOT probability values! Everything is still in logit (transformed in log odds space). Need to transform back.
ypoints = logit.inv(ypoints) # Now you have values bound between 0 and 1
points(xpoints,ypoints,type="l", lwd=2) # type="l" this is a line
points(..., col=) # argument for transparent points