New answers tagged

0

This example used poisson distribution, which defaults to link log, right? I'd expect the random effect, even if ALL groups had one observation, to estimate the variance in the log, whereas the residual error would then measure the error in the prediction including the random effect. Example, 10 observations, first without random effect: summary(pp <- ...


0

Certainly not; both have their place. As already said in the comments, IRT has no way to account for continuous responses. Further, when the sample size is small, CTT may be your only option as IRT (and other latent variable measurement models such as CFA for that matter) is a large sample method (e.g., a common rule of thumb is to have a sample > 1000).


0

This line in the RSS function requires that parameters is of length 3. names(parameters) <- c("beta_value", "gamma_value","delta_value") But in your optimization optim(c(0.5, 0.5), you are giving only a vector of length 2 as starting point.


1

It looks as though you are interested in extracting (i.e., observing) latent variable scores for prediction purposes (i.e., not necessarily to make inferences). Given this, I would not rule out PCA either (while dually noting its similarity to EFA; see link below for more details). Also, since your goal prediction, I would not suggest using CFA, as the most ...


0

This depends whether your goal is only to test the hypothesis you tested, or if you are also interested in estimating and criticizing your regression model. First a simple example: If you are interested in a simple regression model $y=\beta_0 + \beta_1 x+\epsilon$ with possible values for $x$ being the interval $[a, b]$, then the optimal design (at least for ...


2

I'm answering only to the citation at the beginning of the question. I did not consider the analysis in R provided in the question. I admit that the author of the first passage may have made some confusion in terminology. Let us define properly. $v$ is the eigenvector's values. It is the eigenvector from eigencecomposition of the covariance matrix of the ...


0

This can be achieved by computing the $\Sigma$ first from the prcomp() function to make it numerically equivalent to that of the svd()and it turned out after inspecting the source code of prcomp() this can be achieved by the following: sigma <- prcomp.pca$sdev * sqrt(max(1, nrow(S) - 1)) round(sigma, 2) == round(svd.pca$d, 2) # TRUE diag(sigma) %*% ...


3

When the online references become convoluted, it is time to invest in some good books on the topics of mixed effects modelling or longitudinal data analysis. You can pose a question on this forum asking people to recommend such books to you. Going back to your question, the first thing you should note is that your terminology needs to be more nuanced. In ...


0

You can use "emmeans" for planned comparisons - it's just a matter of setting up your own custom contrasts (or comparisons). Note that some people don't adjust for multiplicity when considering multiple planned comparisons, while others do. You can refer to this excellent post by Very statisticious on Custom contrasts in emmeans for details on how ...


3

The only benefit of na.exclude over na.omit is that the former will retain the original number of rows in the data. This may be useful where you need to retain the original size of the dataset - for example it is useful when you want to compare predicted values to original values. With na.omit you will end up with fewer rows so you won't as easily be able to ...


2

I don't know of a way that you can impose a constraint on the variances. However, I can think of two ways you could achieve something similar: instead of fitting (1|A) + (1|B) + (1|C), using (1|A:B:C) combine the three variables A, B, and C into a new variable, D which represents all the unique combinations of the three and fit (1|D) Both approaches ...


2

As Jeremy pointed out, EFA, CFA, and IRT models scores would usually be in close agreement. This is especially true in the case of unidimensional scales, or 2nd-order factor models (since this will take you back to almost the same configuration when working on the higher-order factor). Moreover, PCA, which does not take into account measurement error but is ...


3

most important: is it appropriate to report the results as written above in the log-link format or do I have to convert the Estimates and Std. Errors by exp()? Whether it is appropriate or not depends on your audience. Some people may be happy with the estimates being on the log count scale, but for other it may be better to exponentiate them so that they ...


0

Your title says "with lme4 or nlme package", but your text says How can I achieve this using R packages, perhaps with nlme or lme4? I know that ASREML can do it but I do not have hold and I love R for being robust as well as free. This approach is not based on these two packages, but it is open source and very flexible. GBLUP with arbitrary ...


-1

This sounds like the use case for a correlated random-effects model: https://journals.sagepub.com/doi/pdf/10.1177/1536867X1301300105


3

Model.2 doesn't make much sense because you are asking the software to estimate random intercepts for Subject three times. That is because: (1|Subject) + (X1|Subject)+ (X2|Subject) is the same as: (1|Subject) + (1 + X1|Subject) + (1 + X2|Subject) where the 1 in each term specifies that you want random intercepts for the term on the right side of the |. The ...


2

Although it is difficult to determine general behavior from one dataset, it looks like these data would be amenable to a simple, robust analysis that looks for sudden drops followed a little later by a sudden rise. (I adopt the language of time series analysis in which row numbers are viewed as succeeding points in time.) To find a drop, predict the data in ...


1

Thanks for including the link to the data. This makes it clear that you have observations that cover a variable number of days. I'm going to assume that your 52 week data are 52 normal, 7 day weeks. One approach might be to convert the external variable data to a daily rate. For example, if you have 18 cases in 9 days, the daily rate is 2 for each of those 9 ...


2

Recall that logistic regression with a single predictor $x$ models the probability for $y$ to be TRUE (or 1) as $$ P(y=1) = \frac{1}{1+\exp(-\beta_0+\beta_1x)}. $$ We can now approach your question in two ways. Either you have your estimated coefficients $\hat{\beta}_0$ and $\hat{\beta}_1$. Or you have a fitted model object from a previous call to glm(). Let'...


2

@chl is right. I added it to to the function output to remind persons what hypothesis they are testing because it is often not clear what the alternative ($H_A$) and the null hypothesis ($H_0$) is. So it just tells you what the null hypothesis is and nothing about the acutal result. p < 0.05 means that $H_0$ can be rejected. So in your case the parallel ...


0

Because of this structure, I don’t think I can just do a usual Cox Proportional Hazard model (like coxph()). That's not necessarily true, and if it is true it might not be for the reason you think. The fact that, for technical reasons, many individuals have the same values for areas and concentrations doesn't by itself rule out a simple coxph() analysis. ...


1

Answer in two steps. First, data transform testing. Second, model fit testing. Concerning data transformation, from prior work for the type of data the OP has offered, neither model is appropriate. Body scaling is not linear, so linear models are not as useful as logarithm transformed data and variables, which leads to power function formulas. For power ...


0

The fable package replaces the hts package and produces prediction intervals. It is also much easier to handle the aggregation structure. Here is some code using the same example as in your question (updated to include multiple models). library(tsibble) library(feasts) library(fable) library(dplyr) df <- as_tsibble(hts::htseg1$bts) %>% mutate( ...


0

The tests of the coefficients are for the effect of that variable over and above all the other variables in the model. The ANOVA table tests each variable added one by one, so size by itself, new over and above size, taxes over and above size and new. Notice that $\sqrt{30.558} = 5.528$ but that does not apply for any of the other F and t values. I do not ...


1

The first step, if you are not using unsupervised learning (data reduction), is to show evidence that there is a predictive signal. Fit the full model without penalty and get the likelihood ratio $\chi^2$ test with say 33 degrees of freedom. This corrects for having 33 chances to find something. If this $\chi^2$ is not large you do not have a basis for ...


2

Assume model is: $$y=X\beta+\epsilon$$ where $\epsilon \sim N(0,\sigma^2I)$ The residuals, $\hat{\epsilon_i}$ are estimates of true errors, $\epsilon_i$ The are related to each other by the residual-maker matrix, $M \equiv I-P=I-X(X'X)^{-1}X'$ Therefore, $$Var(\hat{\epsilon})=Var(My)=MVar(y)M'=\sigma^2MM'$$ Since $M$ is idempotent and symmetric this becomes: ...


1

We found the solution Needed to first specify the ordered levels: pgsi_cat_levels <- c("non_prob", "low_risk", "mod_risk", "prob_gamb") And then use these levels: mutate(pgsi_label = factor(t2_pgsi_cat, levels = ...


0

In Python, the Anomaly Detection Toolkit (ADTK) provides really a nice interface and suit of functions. This talk from 2019 provides a walkthrough of the features, but essentially the same material can be found in the examples in the docs. The package provides 13 built-in methods for detection ranging from the very simple, e.g. thresholds, to the complex e.g....


2

ACME is an acronym for "Average Causal Mediation Effects" ADE is an acronym for "Average Direct Effects" Total Effect is the sum of ACME and ADE So the "Indirect Effect" that you are seeking is simply the row forACME


4

Is a Mixed Model Appropriate for Repeated Measures of Multiple Covariates? Yes, you have repeated measures within subjects, and you are not interested in specific subject effects, so a mixed model is appropriate for modelling these data. There are a few things to note here. The structure of the random effects in the two models are not equivalent. The lme ...


2

If the probability of getting a particular result on one test is $0.03$, then when you have two tests, you have two chances to get that result, so it's $0.03*2$ (technically, it's $0.03*2-0.03^2$, but that rounds to $0.06$). We decide whether to reject the null based on whether the p-value is greater than $\alpha$. So when we adjust the test, we can either ...


0

First, make sure you coded your dummies correctly. For example, a categorial variable X with 2 possible values A and B, a dummy should be created for either A or B, not both: D = 1 if X = A and 0 otherwise. Or D = 1 if X = B and 0 otherwise. If 2 dummies are created, one for each value, then the dummies are linearly dependent and OLS just won't run. ...


8

The question is why, not what does the manual say? The answer is that alpha and the P value are different things, and they work in opposite ways. When you get a P value, you are asking "at what alpha value would I reject the null?". In other words, we are going to view our observed |t| as the critical value of the test, and the P value is the alpha ...


5

From the "Details" section of the p.adjust help page (?p.adjust): Details The adjustment methods include the Bonferroni correction ("bonferroni") in which the p-values are multiplied by the number of comparisons. So the output is 0.06 because 2 x 0.03 = 0.06: number of comparisons is multiplied, not divided.


0

Ah. I realized that I missed the point of cross-validation. Cross-validation, in the context with which I'm using it, is just to estimate out-of-sample prediction error. But with the model fit to the entire dataset. What I actually wanted were bootstrapped coefficients, which use resampling to estimate bias of the coefficients. Related, but different. This ...


4

Nice answer from Thomas! I think we may need some more information though before finalizing a modelling approach. For example, the currently proposed model ignores the fact that the response variable was collected on different days. One way I would think about this modelling exercise is like this: We start out with 2 rounds of experiments, so Round can be ...


3

Adding (1|round) as a random effect to prevent pseudoreplication is right, a good article about this subject can be found here. To ensure your model assumes everything correct you could compare the df’s of your model summaries and check if they take account for your pseudoreplicates and are therefore lower in the model with round as a random effect. The ...


4

Yes, AIC (or AICc) are still "generally appropriate". We make the assumption of a Gaussian error (see below). Avoid computing AIC manually if both models are of the same type, that said: Comparing across model types requires attention to detail to make sure that parameters are counted using similar rules, and that additive constants are ...


2

An alternative approach (I agree it is devious, bit it is also interesting) is to transform your function $$y=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3x_2^2 + \epsilon$$ into $$y=\beta_0 + \beta_1 x_1 + \beta_2 \frac{x_2}{3} + \beta_3(x_2-2)(x_2-5) + \epsilon$$ This is the same quadratic polynomial but now you have $\hat{y}_{x_2=5} - \hat{y}_{x_2=2} = \...


3

I can recommend the new shiny based packages to everyone, it makes data visualisation and inspection interactive and thus easier than writing code in R espacially in the beginning. A good example would be ggplotgui


0

This is an interesting question that has different answers according to the context. I agree with what have answered you before, so here I will focus more on the context. It is normal that you have been confused when looking for McNemar's test interpretation on the Internet. The main reason is that both the historical origin of the test (Genetics), as well ...


0

Even with only 100 observations, then assuming that the data are missing at random or missing completely at random, it is likely that a pricipled approach to missing data such as multiple imputation will provide much better resuts that removing rows/columns or any kind of single value imputation. The general approach to multiple imputation is: create ...


2

To use matching, you need a binary variable that is the purported cause. You do not have that here. length is the purported cause, and it is not binary. You cannot use matching for this type of analysis.


0

When a model with breakpoints fits the data well, then the choice of h should not matter too much, as long as it is large enough to fit the model reliably on every subset. With a model like this (linear trend + strong autocorrelation), you will at least 25 observations per subset, maybe more. This is also what I get when I set.seed(1) before running your ...


1

I'm on my phone so can't check any details, but my best guess is that the upper bound of the envelope (hi) is NA when the shading stops. You can use as.data.frame to convert the envelope to a data.frame and inspect the values. This may be a mathematical problem with the estimator and not a computational problem / bug, but I can't investigate right now.


0

You can use ConsReg package. cran.r-project.org/web/packages/ConsReg/index.html It's very easy to use


0

When writing it up, you should interpret the model as follows. Your variable "FIRST_cat" was significant, you can reject the null hypothesis that it is equal to zero. If that variable is "high" it has a positive effect, otherwise a negative effect. Furthermore, that means that the probability of "InsomniaT3" being "high&...


0

The following books on Statistical modeling using R might be helpful Title : An Introduction to Statistical Learning: with Applications in R ISBN: 9781461471370 This book presents some of the most important modeling and prediction techniques, along with relevant applications. Topics include linear regression, classification, resampling methods, support ...


2

If you check the warning posted by R after fitting your lmer model, Mod1, you'll notice the dreaded singular fit message. You can google that error message to learn more about it but suffice it to say that you should not trust a model which comes with this warning. It seems that there's virtually no variation among the random intercepts in your model so this ...


1

Take note of the comments by @DimitriyV.Masterov. Difference-in-differences works well with repeated cross-sections data. Your equation would look something like the following: $$ y_{ist} = \alpha + \gamma S_s + \lambda A_{t} + \delta (S_s \times A_{t}) + \theta X_{ist} + \epsilon_{ist} $$ where $y_{ist}$ denotes movie $i$ in society $s$ at time period $t$. ...


Top 50 recent answers are included