New answers tagged

0

The scale() function doesn't translate correctly outside of the fitted model, because it depends on the values of the variable being scaled, which are not the same in the original data frame as in the reference grid. To get meaningful results, you need to scale the variable ahead of time. Your results: > lstrends(m1s, ~Type, var="Var", adjust = "tukey", ...


0

You could use a two-part mixed effects model for semi-continuous data. This combines a logistic regression for the dichotomous indicator that the outcome is zero or not with a log-normal model for the continuous part. Two-part mixed effects models are available in the GLMMadaptive package I’ve written - you can find a sample analysis here.


1

Rather than transforming data to fit a model, I suggest using a model that fits the data. Quantile regression does not assume homoscedasticity and quantile regression for multilevel (or mixed) models now exists. There is an R package lqmm available from CRAN that should work.


2

When you include the term (trialtype | subject) in your model, and given that trialtype is a categorical variable with four levels, you include four separate correlated random effects per subject (i.e., you postulate an unstructured covariance matrix for your random effects). This is a complex model with many variance components, and it's no surprise that it ...


1

You don't explain what statistical model you are fitting. I am going to assume that you are fitting a Poisson or negative binomial glm with a log-link to the number of Type A utterances. I assume that you want to compare counts between instruction types and possibly between other experimental conditions as well. As always, the appropriate analysis depends ...


1

For the first part of your question, reporting p values by themselves provides little useful information. This thread is a good introduction to the limits, even the dangers, of relying on p values. A reader is going to want to know details about the size of the sample on which the results are based, the magnitudes of individual effects, of differences among ...


2

You have indeed a complex type of design and a rich dataset to explore it. A couple of considerations: To explore the following (nested) structure, you will need to put the data in the long format with measurements over time in the same book and type put one underneath the other, and also create a time variable. It would be expected that sales of the same ...


0

your auto.arima model is probably WAY over-modelled as the sum of the ar coefficients is approximately 1. suggesting non-invertibility perhaps due to the near cancellation as a result of the ma(1) coefficent . GLS appears to be slightly smarter as it CORRECTLY flagged the model . Post your data and I will give you more definitive corrections. Your data ...


0

I used the same strategy and I had exactly the same issue. It seems that parameters > 1 are possible for explosive time series. For the parameter that should be > 1 from your auto.arima results, try to set it to 0.99 in your GLS model without using the argument fixed = TRUE. The model will use the value that you provide as starting value and should finally ...


0

You can get the proportions of sum of squares using aov, but this assumes you have the correct model. This is almost the same as using least squares. Define all fixed effects as factors. n=200 set.seed(123) var1=rnorm(n,100,5) var2=rnorm(n,300,2) week=1:n%%7 area=1:n%%4 salesvolume=round(100+var1*20+var2*35+ifelse(week>3,150,0)+ifelse(area>2,100,0)+...


1

If you standardize your variables, then the fixed effects coefficients you will tell you about the relative importance of your variables compared on the same grounds. I.e., one standard deviation increase of var1 increases salesvolume by X, whereas the same increase in var2 increases salesvolume by Y.


1

In GLMMs you do not have an analogue of multivariate error terms for which you can define such a correlation structure. A potential way to achieve something like this in GLMMs would be to use observation-level random effects, and define such a correlation structure for their variance-covariance matrix. I think this should be provided by the glmmTMB package.


8

A couple of points: Even there are exceptions, the general rule involving interactions is that you need to include the lower order terms. That is, in a model in which you want to include a 3-way interaction, you also need to include the main effects, and all 2-way interactions. In that regard, model mod2 does that. In general interactions are complex terms ...


2

The model you are fitting postulates that per SCHOOL you have a different random intercept per RACE (BTW perhaps it would be more easier to interpret if you exclude the intercept in the random-effects part, i.e., lmer(TEST_SCORE ~ PRIOR_TEST_SCORE + (0 + RACE | SCHOOL))). However, the formulation you have used also postulates that the random intercepts for ...


2

Ben is correct. But if interested, the emmeans way is emm <- emmeans(model, “bodypart”) contrast(emm, “eff”) This compares each EMM with the grand mean. If you have more than one fixed factor in the model, this’d be the way to go.


1

Yes emmeans can do this via emmeans(fitted_model, ~ bodypart). However, in this particular case you can get the same results by fitting the model weight_change ~ 0 + bodypart + (1|subject) (a -1 would work just as well in place of 0); this instructs R to suppress the model intercept, and in this case would correspond to the model $\Delta w_{ij} = \beta_i ...


0

Exactly, as Dimitris mentioned, RNA-Seq tends to follows a Poisson distribution or a reparameterized Negative Binomial, where the NB now models "over-dispersed Poisson" (variance > mean) data rather than the traditional "trials until failure" model. There are a few different R packages in Bioconductor for RNA-Seq analyses that might be worth checking out. ...


0

Indeed, you can simply include these covariates in the formula of lmer() as fit <- lmer(y ~ x + z + w + (1 | i) + (1 | o), data = dta) and you will get the estimates for their effects. In particular, the model that you will fit is $$y_{io} = \beta_0 + \beta_1 X_i + \beta_2 Z_o + \beta_3 W_{io} + u_i + v_o+ \varepsilon_{io},$$ with $u_i \sim \mathcal N(...


2

RNA-seq data are (normalized) count data exhibiting a mean-variance relationship. The linear mixed model works for normal data/error terms that don’t have this relationship. You could work with a Poisson or Negative Binomial mixed model instead.


0

Posting an edited version of my old comment as an answer: I second this comment by @Kodiologist "...I think that anything useful to the analyst as a modeling diagnostic will be useful to the reader, too, to help the reader decide if you made good modeling decisions". Withholding potentially useful information because readers may not be statistically ...


0

A short answer to the updated question: your model is too complex. You will need to either simplify your model (typically one removes more complex random effects) or to switch to a Bayesian approach to fitting your model. A more thorough explanation can be found here.


0

Dimitris comments are useful, but I think the fundamental issue here is the assumption that you have to simplify your model. If your initial hypothesis was x ~ y * z, fit this model, check residuals, consider if you might be overfitting (i.e. do you have enough data to fit this model), and if everything is fine report the results and move on. Removing ...


1

A couple of points: You have not stated what test exactly are you performing when you run ANOVA. I would suggest that you best fit your linear mixed model using restricted maximum likelihood (REML), which is the default in most software, and then you do an F-test for the interaction using the Satterthwaite's degrees of freedom. In R you can do this using ...


4

To better check for over-dispersion you can use the simulated residuals provided by the DHARMa package. If you want to account for over-dispersion, you can use a negative binomial mixed effects model. This is provided by glmer.nb() in the lme4 package and also by mixed_model() of the GLMMadaptive package. The glmer.nb() fits the model using the Laplace ...


2

"But when I ran the code, I realized that model2 gave the exact same output as model1." If you check the estimate of variance of error term and random intercept, two models give you different results. Also the variance of fixed effect are different, degree of freedom are not the same. For the fixed effects (intercept and slope), fitting two models should ...


1

The ANOVA is Type II sum of squares and the method you used below is Type III sum of squares. In general, they don't agree. Models for independent data give the same results. For mixed models, they differ. Type II and Type III sum of squares handle the random effect differently. The Type III method assumes, (incorrectly in most cases) that if the "true" data ...


0

First, the P values are adjusted for the fact that you are simultaneously testing six comparisons. Second, you’re comparing apples with oranges — tests of nested models aren’t the same as tests if specific linear functions if predictions. For reporting, just give the P values. Declaring things “significant” or not based on a .05 threshold is poor ...


2

A couple of points: The fixed-effects part of the model specifies the mean structure. Based on the design of your study and your research question, you would probably want to assume that there is a difference in the average longitudinal evolutions between the two groups. To achieve this you would need to include the interaction term between your follow-up ...


3

The translation from glmer() to glmmPQL() is correct. The error you receive indicates that the optimization algorithm behind lme() that is internally used in glmmPQL() did not converge successfully. You could try setting as optimization algorithm the optim() function instead of the nlminb() (the default). But note that the glmmPQL() algorithm (penalized ...


4

A couple of points: Yes, you can compare these two models because model md.logistic1 is nested within model md.mm. With the anova() function you do a likelihood ratio test to compare these two models. In particular, the null hypothesis you are testing is that the variance for the random effect you have included for the grouping factor trainId is zero versus ...


2

You can also obtain coefficients from a generalized linear mixed model that have the desired marginal interpretation. These coefficients are provided by function marginal_coefs() of the GLMMadaptive package. For an example, check here.


Top 50 recent answers are included