New answers tagged

0

I realize this is at the study-level and your first column looks like the number of subjects in each study. Do you have another column that contains subject-years for each study? This will be helpful for calculating the exposure-adjusted mean since it is likely not every subject was followed for the same amount of time. My suggestion would be to feed this ...


1

Multicollinearity is not a problem. In the world of mixed effects modeling, it's thought of as within, vs between, cluster variation, and is assessed using the intra-class correlation (ICC). However, higher ICC means that you have less power than if the ICC was lower. This is why, if all else is equal, you would rather have fewer ratings from more ...


2

What you need is a way of checking your mathematical analysis by means of a simulation that relies as little as possible on that analysis. Otherwise--if you base the code on your calculations--the code will give results that are only as correct as the calculations themselves. It wouldn't be much of a check at all. In this instance, rejection sampling is an ...


9

Regression analysis accommodates situations where the explanatory variables/covariates can have any values, including zero. There is no particular model needed for this --- it can be implemented in almost any model. If there are lots of zeros for a particular covariate, the only issue this creates is that it affects the leverage of the data points and one ...


0

One possible approach for the requested graphical results: Start by simulating points in the rectangle with vertices $(0,0)$ and $(1,1)$ so that $X\sim\mathsf{UNIF}(0,1)$ and, independently, $Y\sim\mathsf{Beta}(2,1).$ Then throw away points not in $A.$ Then plot histograms of the remaining $X$'s and $Y$'s, which will be negatively correlated. Simulation ...


3

Your value for the CDF of $X$ is incorrect, but perhaps this is just bad notation on your part. What you have written is in fact the PDF of $X$. Simulate $U_1$, a standard uniform random variable. Plug this value into the inverse CDF $F^{-1}_X(U_1)$ of the marginal density of $X$. (This inverse may not be tractable, although it exists, and the Acceptance-...


4

The two commands you use estimate two different models. To illustrate, I use the Grunfeld data set; you can think of its firm variable as your Country variable and its year variable as your Year variable. Your 1st model: You effectively estimate a two-ways fixed effects model where the time fixed effect are explicitly modelled via dummies (the part +factor(...


1

If the structural break only affect a subset of the model parameters, you could use the whole time series to model all other parameters (carrying over info) and only use post-break info to infer the new parameters. The r pacakge mcp can do this. As a simple example, here's a model with an intercept change but with no change in 2nd-order autoregressive ...


0

Ideally this analysis would include the observation time when a subject is eligible to fill 'scrips at a particular pharmacy. The linear model for the response would be: $$ \log \left( \frac{\text{Events}}{\text{Obs time}}\right) = \beta_0 + \sum_{i=1}^p \beta_i X_i$$ where $X_i$ is the vector of covariates adjusting the intensity $\lambda$. In a Poisson ...


0

What you are describing doesn't really make sense. To reduce the confounding effects of a variable, you need to find units within each treatment group (i.e., each level of Var1) that have similar values of the variables (other than glucose). The outcome (glucose) should have no role in this part of the analysis. It should only be used in the final step of ...


1

I'm not 100% sure I understand the question, but I will do my best by interpreting what you might be interested in. Question 1: How do I explore the multicollinearity in my data? Firstly, figure out if your variables are multicollinear (compute VIF or effective rank). If they indeed are multicollinear, it may be of interest to know the structure of ...


1

Regarding exponential smoothing yielding autocorrelated residuals, it depends on how strong the autocorrelation is (magnitude, not only statistical significance). If it is strong, consider using another model or modelling the residuals and adjusting your forecasts accordingly. Thank you very much dear Mr. Hardy. I managed to fix this by playing with $p$, $d$...


0

I prefer to equivalently view the AFT model in terms of a generalized linear model like logistic regression or Poisson regression. In these models there is no "error term." There is simply a likelihood and a parameterization of the mean that is estimated through maximum likelihood. Most documentation describe an AFT model in terms of a log scale ...


1

how can you interpret the intercept of the output of the survreg model other than as the log of the scale parameter of the distribution? Does it tell us anything useful about the data from a biological perspective...? This difficulty in interpreting the intercept is inherent in AFT models, as you note. For symmetric error terms like that in a log-normal ...


2

The title question "How many variables are needed for a good multiple regression model?" doesn't have a general answer. In principle good models can be found with any number of variables, and for sure any model becomes worse if variables that are in fact important are removed. A major problem with a large number of variables is that estimation of ...


1

If you keep the parameters $n$ and $p$ of the distribution fixed at $n=14$ and $p=\frac 1 4$ then of course the variance will also be constant. You have to vary on of the parameters. I assume that you want $p$ to vary. Then you can plot this variance using the curve function in R. See also its help page (?curve) for more details. curve(14 * x * (1 - x), ...


0

For left figure, this can be accomplished by using predict.gam/predict.bam using type = 'terms'. You just need to set up a new data for each pm10 and lag. For example, for lag0 pm10 <- seq(min(pm10, na.rm = T), max(pm10, na.rm = T), length = 1000) l <- rep(0, 1000) newdata <- data.frame(pm10 = pm10, lag = l) If there are multiple smoothing terms, ...


2

It is a bit suprising this question is asked as (almost?) all textbooks on panel data econometrics introduce the Hausman test (Hausman (1978)) as the standard test to discriminate between the fixed effects and random effects specification early on. Also the plm package's first vignette [1] introduces this test; it is implemented in phtest. Picking up the ...


0

You probably mean bootstrap aggregation (a.k.a. bagging) combined with time series techniques such as ARIMA or exponential smoothing. There is a paper* and another one** about ets + bagging and how it works very well on M3 competition data. See this section of the "Forecasting Principles and Practice" textbook for a brief overview with R code and ...


3

With data. First, it is important to know how the var.test in R works, if you do have the (normal) data. Consider data as follows: set.seed(2021) x1 = rnorm(20, 19, 5); x2 = rnorm(30, 10, 8) stripchart(list(x1,x2), pch="|", ylim=c(.5,2.5)) With respective sample variances $S_1^2, S_2^2.$ The test statistic for testing $H_0: \sigma_1^2/\sigma_2^2 ...


1

I agree with the accepted answer. Just a comment on your code: In your trainControl function you specify method="repeatedcv", number=k.folds, repeats=3. Repeatedcv will automatically create 3 different partitions of 7 folds. However, you also specify index=df1.folds As far as I understandd, this conflicts with the above since it provides your ...


2

Your data is noisy, there is no clear pattern. The wide intervals tell you about the uncertainty, the intervals for the predictions are wide because the model is uncertain. If you had a clear, repetitive pattern that could be predicted, they would be narrower. For example, here you can find an example of SARIMA results for such time-series from the great ...


1

Assuming that you wish to use a logarithmic link function (which is the default for what the glm function is implementing here), your proposed linking equation is: $$\mu_i = c(\boldsymbol{\beta}) \exp(\boldsymbol{\beta}^\text{T} \mathbf{x}_{(i)}),$$ which can also be written as: $$\log(\mu_i) = h(\boldsymbol{\beta}) + \boldsymbol{\beta}^\text{T} \mathbf{x}_{(...


0

Your question is likely one of the following two: Is there significant difference in the numerical column across days, controlling for the groups A, B, C, D? Is there significant difference in each group's occurrence across days, controlling for the numerical column? In the first case, you can run an ANOVA to test the hypothesis. In the second case, you ...


1

From the the package documentation glht extracts the number of degrees of freedom for models of class lm (via modelparm) and the exact multivariate t distribution is evaluated. For all other models, results rely on the normal approximation. Alternatively, the degrees of freedom to be used for the evaluation of multivariate t distributions can be given by ...


0

This isn't really different from the issue in any regression model with interaction terms. The only difference is whether baseline predictor values are included in the baseline hazard (Cox) or contribute to an intercept term (other regression models). When you write in R with its default treatment coding of predictors: outcome ~ Control * Gender where ...


1

I intend to get one metric per test that I can plot for all the 100 models to show that the models from which I am fetching variable importance are good enough. Getting a single metric per test will require some thought and care, as illustrated in the following excerpt from your predicted survival curves, the 8th and 10th rows of your predicted survival ...


1

In this instance, and in general, what we are doing is shuffling the data into a random ordering. This has the effect of assigning each observation to one of the treatment groups at random, or, for a continuous variable, randomly reordering the observations with respect to the continuous variable breaks the relationship between the response matrix and the ...


5

There's a single expected value in the chi-squared that's below 5 (it's about 4.6). You can see it if you use chisq.test(gt)$expected. R has implemented a common rule of thumb to check that the asymptotic chi-squared approximation will be reasonable -- one that is sometimes overly conservative (and, rather less frequently, not cautious enough). There are a ...


8

Note that in the model summaries, the regression coefficients and standard errors are the same. In the comparisons outputs, the estimates and SEs are the identical. Also the t ratios in the first are identical to the z ratios in the second. The difference is in the degrees of freedom, which are taken to be 27 in the first model and infinite in the second. ...


0

Multivariate linear regression can be handled directly by the standard lm() function in R. As Fox and Weisberg say (page 3): Multivariate linear models are fit in R with the lm() function. The procedure is the essence of simplicity: The left-hand side of the model is a matrix of responses, with each column representing a response variable and each row a ...


4

Disclosure: I wrote the samc R package used in this answer This answer is more of a supplement to Stephan Kolassa's answer In it, he showed how to construct a transition matrix representing the problem: Image credit: Stephan Kolassa's answer Now, as indicated in comments/another answer, this matrix can be simplified a variety of different ways, but I'm ...


1

Seems very close to Shannon - related theorems. If you posit "HTH" as your "end of message" string, you want to estimate the chance that "HTH" shows up in random data. I suspect a little digging into his work will provide the equations/ formulas of interest. And because I can't resist, "HTH"


0

I assume that this is because you're not specifying the fold IDs for your observations and cv.glmnet, which is used by the model you are fitting, is randomly assigning observations to folds and a different number of these calls to cv.glmnet is happening during the optimisation/fitting of the mixture of GLMNET models. The help for FLXMRglmnet() does tell you ...


1

There is no penalization applied by default to the parametric terms of a GAM fitted by {mgcv}. If you tried to fit: gam(y ~ x + z) you would get back the equivalent (up to details of the actual algorithm and implementation of course) of glm(y ~ x + z) because these terms are not subject to any penalization. However, what you say is correct; that is a ridge ...


1

Your base model is incorrectly specified; factor by smooths must have the by factor included as a parametric categorical term in the model, hence you need: gam(resp ~ species + s(x1, by = species) + s(x2, by = species) + s(x3, by = species) + s(location, bs = "re") This allows for the mean of the response in each level of ...


0

You can try Monte Carlo simulation: Generate a large number of datasets with the properties that you know about (initially normal, truncated 20%, anything else you know or can reasonably assume?) Use some goodness-of-fit measure to compare these datasets to the one given (you can start with simply plotting and eyeballing) Focus on the things you know about ...


1

If you can safely assume that your underlying data is normally distributed, then as Otto Kässi writes, you have a truncated normal distribution. If you know where it was truncated, this is good (and with 800 data points below the point of truncation, simply using the maximum observation will likely be a sufficiently good estimate of it, and any uncertainty ...


2

Here is part of what ?wilcox.test will give you: if both ‘x’ and ‘y’ are given and ‘paired’ is ‘FALSE’, a Wilcoxon rank sum test (equivalent to the Mann-Whitney test: see the Note) is carried out. In this case, the null hypothesis is that the distributions of ‘x’ and ‘y’ differ by a location shift of ‘mu’ and the alternative is that they differ by some ...


1

Mixture distributions make the most sense if we have grounds to suspect that the data belong to multiple different subpopulations, as per the textbook you cite. Usually, we don't know which of the subpopulation a given instance belongs to - and we may not even know how many such subpopulations exist. As such, we could even consider straightforward ANOVA as a ...


2

The chance of A being $>2\sigma$ from the mean is $2.28\%$. pnorm(12, 10, 1, lower.tail = FALSE) Thus the chance of A & B (2 independent events both occurring) is equal to $A*B$ or $0.052\%$ of the time. Therefore expected frequency is once every 1932 times (agrees with the simulation results reasonably well). The medium and mode is a bit more ...


12

There is a fun way to answer this problem using martingales, and in particular using https://en.wikipedia.org/wiki/Optional_stopping_theorem. I first saw this trick in the book A First Look at Rigorous Probability Theory by Jeffrey S. Rosenthal, in the martingale chapter. (I don't have the book in front of me at the moment but I'll edit and add a page or ...


2

You need to distinguish between the instantaneous hazard as a function of time, $h(t)$, and the corresponding cumulative hazard $H(t)$. The latter is the integral of the former up through time $t$, $H(t)= \int_0^t h(\tau) d \tau$. As instantaneous hazard is necessarily non-negative its integral, the cumulative hazard, is necessarily non-decreasing. That ...


1

Final formula appears to be sum(i = 1 to length(pattern) : if (first i flips of pattern match the last i flips of pattern) then P(series of i flips matches first i flips of pattern)^-1 else 0) (For a fair coin, this reduces to Conway's algorithm; it can probably be proven by a similar method.)


31

At any given point in the game, you're $3$ or fewer "perfect flips" away from winning. For example, suppose you've flipped the following sequence so far: $$ HTTHHHTTTTTTH $$ You haven't won yet, but you could win in two more flips if those two flips are $TH$. In other words, your last flip was $H$ so you have made "one flip" worth of ...


0

ARMA model coefficient Interpretation Regarding the AR part, in my view, them have a purely correlational interpretation only. Therefore them are transformation of total or partial linear correlation coefficients and maintain them interpretation too. In the case of $MA$ part a non observable series is involved (errors) and I'm not sure if the same ...


17

First, you can refactor your R code to be (IMHO) a little more legible, also using pbapply::pbreplicate() to get a nice progress bar: n_sims <- 1e5 library(pbapply) results <- pbreplicate(n_sims,{ flips <- NULL while(length(flips)<3 || !identical(tail(flips,3),c("H","T","H"))){ flips <- c(flips,...


7

Here is a somewhat clumsy brute-force method to obtain the probabilities and order statistics. Getting the mean will take more work. So first just generate the possible sequences and associated probabilities where "HTH" are the last 3 flips (with that sequence not occurring previously). Then look for patterns. For integer patterns the go to ...


1

I think your logic here is correct. If you model the event times parametrically you can produce a smooth hazard and cumulative hazard function. When the event times are modeled non-parametrically the hazard function is a blip function. You could try plotting this with vertical bars rather than a smoothing function.


1

spatstat.core::CDF() can be used to to create a cumulative density function from a given output from density(). set.seed(123) x <- rnorm(10000000) x_density <- density(x, n = 10000) x_cdf <- spatstat.core::CDF(x_density) sds <- c(-2, -1, 0, 1, 2) names(sds) <- sds # check cdf at different values setNames( x_cdf(sds), sds) #> ...


Top 50 recent answers are included