New answers tagged

0

Splitting up the long data entry command should get part of the way there, in case that command causes anyone errors as it did for me: data1 = data.frame(B = c("m","m","m","m","m", "m", "f","f","f","f","f"), G = c("s","s","s","u","u", "u", "k","k","k","r","r"), ZN =c(78,82,34,67,98,56,37,45,27,18,34), GFR=c(120,100,90,60,100,110,100,90,95,87,96), g1 = c(35, ...


1

You should test for Granger (non-)casuality in the underlying VAR-model in (log) levels, rather than the VECM representation of it. This code will reproduce (part of) the example from Dave Giles' blog post using lm() and the lmtest package. I found the data here (where another implementation in R is also suggested) library(lmtest) library(xts) ### READ ...


0

To subset your dataset, assuming it is a dataframe, you can use rain1 <- rain[rain$country == 1, ] Programming questions are better suited to stackoverflow by the way.


1

Your response score does not follow any generalized linear model family distribution. If the 30 categories that make up the score were statistically independent and of equal difficulty, then the score would follow a binomial distribution but, as it is, neither of these things is true. In particular it appears that the categories vary considerably in ...


1

The motivation of Calinski and Harabasz in the original paper is based on Analysis of Variance theory which involves sums of squared distances from the (cluster) means. This corresponds to k-means (and Ward's method) but not to any clustering method (such as k-medoids) that isn't based on the sum of squared distances to the cluster mean. As far as I know (...


1

This is the classical difference-in-differences (DiD) model. It is a standard interaction model. You interpret the estimate of the intercept as the mean of your outcome for the control group (i.e., men) in the year(s) before the law was enacted (women = 0 and law = 0). The coefficient on women is the expected mean change in $y$ between treatment and ...


0

When you include a term like (1|country/region) in your model, that is equivalent with including the following two terms: (1|country) and (1|region:country). See the Cross Validated post Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4? for details. This means that you can replace a model formulation like: ...


0

This is an old post, but there is a flaw both in the code and the solution that one should be aware of. When you define design = svydesign(id = ~ psu, strata = ~ stratum, weights = ~ sample.weight, nest = TRUE, data = {df %>% filter(!smoker, age.months >= 50*12, age.months < 85*12)}) you are getting rid of all the ...


0

@horseoftheyear answered this one for me: So the types variable is a count for the number of different types of ice cream available at a location (from your description). Hence the estimated coefficient represents the effect of a one unit increase in types on the outcome variable, revenue. Ergo an additional type of ice cream is associated with an increase ...


0

You make y with c(before, after) and group with c(rep(0, length(before)), rep(1, length(after))), so this is pretty easy to structure. However, I see a few reasons not to bother with this test. 1) What will you do if the test comes back and says that the variances are unequal? What about if the test does not report a significant difference? Would you ...


1

This would be true if the model were estimated by conditional least squares, but here it is likely estimated by maximum likelihood. The initial values are estimated together with the model coefficients. For immediate details of the R implementation of ARIMA, see the help file for the R function arima from the stats package (scroll down to the section ...


2

The issue here is that the last coefficient in your model is being dropped because of collinearity. Essentially, the model is recognizing that two or more of your predictors are identical, or perfectly predicted by the combination of the other two predictors. This means that you cannot include all three terms in your model. This question will be useful ...


1

What is your aim of clustering? It all depends on what you want to know. If your aim is to reproduce the given classification, the ARI will tell you how well you do that. However I wonder why then would you want to cluster these data if you know the true clusters already? What you do here gives some sensible validation information for your clustering if, ...


1

I think you are doing it right... "I am not sure how to express the covariance matrix": but you specify variance to be 1, right? This is what you want? You just need to draw each time conditionally on drawing $X's$ in the right order: first draw first 10 $X's$, then the other five conditional on first five $X's$, and lastly draw $Y's$ condtional on 1, 5, 7,...


0

From my understanding, transfomring categorical to dummy variable may not help you in the end. https://online.stat.psu.edu/stat504/node/169/ This article shows specific steps about poisson regression modelling process in considering overdispersion and engineering variables. Using offset may help to address over-dispersion, the article is worth a good read,...


0

recently something occurred to us when we also perform test. The fatality rate does not quite describe number of death by certain disease to start with. We performed test against Covid-19 patient who have/not have chronic disease. it turns out patient who have chronic disease would have more chance to develop Pneumonia and acute respiratory. it might not ...


0

The warning should go away if you use family=beta_family(). As for the significant KS test, perhaps the large amount of observations you have makes it very sensitive to the slightest deviation from uniformity? (Also, specifying quantreg=T in the simulateResiduals() will (eventually) give you a more readable residual vs predicted plot.)


2

You do not compute the point estimates wrong, but your manual procedure of computing IV estimates in procedure 3 produces the wrong standard errors, so that you should trust the standard errors (and hence significances) of the results of ivreg. The issue is mentioned (but not fleshed out) here. It is discussed in many introductory econometrics textbooks like ...


0

If you're working with lots of time series, it is much easier to use tsibble objects than ts objects. Here is an example. library(tidyverse) library(tsibble) library(feasts) library(lubridate) library(stringr) set.seed(999) # Create data frame with 47 columns df <- matrix(exp(rnorm(47 * 64)), ncol = 47) colnames(df) <- c(letters, LETTERS)[1:47] df &...


0

Depending on the restriction you set. You can set any of them as 0 or set the sum as zero.


0

In principle Monte Carlo works as follows for integration: $$I=\int_a^bf(x)dx=\int_a^bf(x)\frac{b-a}{b-a}dx=(b-a)\int_a^bf(x)\pi(x)dx= (b-a)E[f(x)]$$ where $\pi(x)$ is PDF of uniform distribution on $[a,b]$. Therefore, we can approximate the integral by sampling from uniform: $$I\approx \frac {b-a} n \sum_{i=1}^n f(\xi_i)$$ where $\xi\sim U(a,b)$ Now you ...


0

A suggestion, employ interpolation plus a randomization component. This is likely required to effectively address one (or more) of what I would frame as a missing data connection link. Here is a source discussing various methods in use: Sixty-eight papers (83%) described how they dealt with missing data in the analysis. Most of the papers excluded ...


0

There may be confusion between inverse (cdf) sampling which means simulating from a distribution with cdf G by taking the quantile transform of a uniform $$X=G^{-1}(U)\qquad U\sim\mathcal U(0,1)$$ and importance sampling which means replacing the expression of an expectation under a given distribution with pdf $f$ by an identical expectation under a ...


0

Thank you EdM for your answer about the choose of the time point. For the code itself, I've finally found something like a slide presentation of the book of F. Harrell, made by himself and i've found the solution, written below. We have to declare the number of months in the cph function with time.inc =5*12 for 5 years for example (5 years of 12 months). ...


2

The interpretation of the coefficients is as follows: time = -0.16462 denotes how much the average RNA viral load changes (here decreases) with unit of time for subjects with GENDER = 1. GENDER2 = -0.32416 denotes the difference in average RNA viral loads between subjects with GENDER = 2 and GENDER = 1. time:GENDER2 = 0.07643 denotes the difference between ...


4

SPSS LOGISTIC does not handle sampling weights correctly for computing standard errors. If you have weights $w_i$ for each observation, SPSS will work out the loglikelihood contribution $\ell_i(\beta)$ for each observation, and maximise the weighted sum $\hat\ell(\beta) = \sum_i w_i\ell_i(\beta)$. So will R. The point estimates will agree exactly. SPSS, ...


3

Importance sampling approximates the integral $\mathbb E_{\pi(x)}[h(x)]$, so $\pi(x)$ is the original distribution. It makes sense if the original problem is to calculate an expected value of a function wrt distribution, $\pi(x)$. Instead of sampling from $\pi(x)$, we sample from a proposal distribution, $g(x)$ and calculate the expected value of the ratio $\...


0

This looks like a homework problem, so I'll be a bit coy with my answer. If you're allowed to use adjusted $R^2$, you can do it. $$R^2_{adj} = 1 - (1-R^2)\dfrac{n-1}{n-p}$$ Here, $p$ is the number of parameters estimated in the $\beta$ vector including the intercept. The printout gives you $R^2_{adj}$. The printout also gives you the number of parameters ...


1

Yes. Both those options work with mgcv if you don't have factors in the model. If you do have factors in the model and want to get rid of the constant, gam() has an argument drop.intercept, which you can set to TRUE to force no constant in the parametric part of the model.


2

Why are my estimated coefficients different? They should be different. You're no longer modeling count data. You're modeling rates. The offset is just like any other predictor in a linear model, the coefficients of the other terms shouldn't change when it is uncorrelated. No. The offset is not your typical covariate. The offset is a predictor whose ...


2

I can see two main options: random splines for the combinations given by the interaction of the factors, or by factor smooths, with the interaction of the factors used as the by variable I would create f1f2 in your data frame from the result of interaction(f1, f2, drop = TRUE), then you could do: Option 1 y ~ t2(x, y, f1f2, bs = c('tp', 're'), d = c(2,1),...


2

The error term $\epsilon_i$ is something that you don't observe. You assume a linear model $y_i = \alpha + \beta x_i + \epsilon_i$ for some $y$ and $x$, then you can estimate $\alpha$ and $\beta$ using lm(y ~ x), you can then compute the predictions $\hat{y}_i=\hat{\alpha}+\hat{\beta}x_i$ and the residuals $e_i = y_i-\hat{y}_i$.


0

Following my comment above, here is a reproducible answer: library(CARBayesST) ################################################# #### Run the model on simulated data on a lattice ################################################# #### set up the regular lattice x.easting <- 1:10 x.northing <- 1:10 Grid <- ...


0

First, should it be model_lok<-lokerns(y~x, data=data_test) instead of model_np<-lokerns(y~x, data=data_test)? I ran the example but get a different result I think the difference here might be (i) the decimal different between your input and (ii) the seed for random number used in R.


1

Another potential way of solving this problem is with a Poisson Distribution. In this case $\lambda$ is 1/1000. So the probability of of observing 0 events with a single trial is: ppois(0, 1/1000) #[1] 0.9990005 Thus it is 99.9% chance that no events will occur. Now if we are interested in testing a 1000 independent trials then the overall probability ...


3

Ah. Well, this seems to be a very badly formulated way of saying: We perform 1000 independent trials, each with a success probability of 0.001. Then we have a chance of 37% to observe no successes at all. This is a true statement, pbinom(0,1000,1/1000). Note how the statement in the question does not mention the number of trials, only the proportion of ...


1

With this type of model, the choice of time point prediction is up to you. There are no strata with potentially different baseline hazards over time, there are no time-dependent covariates, and I assume that the proportional hazards (PH) assumption is met, so the hazard ratios among different combinations of covariate values are the same at all times. That ...


1

A GARCH model assumes the standardized errors (shocks, innovations) are i.i.d. with zero mean and unit variance. After having fit a GARCH model, it makes sense to test whether this is the case. Some common checks are to examine presence of autocorrelation and/or autoregressive conditional heteroskedasticity in the standardized errors; under the i.i.d. ...


1

residual deviance equals to -2 times log-likelihood, and it also equals to the sum of squared residuals of the regression model I fit. The second statement is only true in general for ordinary least squares/linear regression/MLE with Gaussian residuals ... however, there are a variety of different ways of computing residuals. The residuals that are stored ...


3

The difference between these two plots is whether or not you constrain the model to satisfy proportional hazards -- the difference is in the Cox model, not in the adjusted curve The first model has four curves for the four groups. The hazard for individuals in different groups is proportional; you can see this because the curves all have steps at the same ...


0

You need to square it before summing, otherwise it doesn't make sense.. And you should not sum up the cross-validation errors, they should be keep separate. Something like below should be ok, I simplified some of the code: errorfun<-function(testData, prediction){ return(sum((testData[,6] - prediction)^2)) } cv<-function(data, model, d){ idx = ...


2

I have no problem with @GregSnow's answer (+1), However, from discussions in comments, it's not clear to me that you understand what one can and can't claim after doing a statistical analysis---especially when the result is not to reject the null hypothesis that all means are equal. In your experiment, you have effects of soil type and amount of rain in ...


1

One way to understand the solution to the problem - the answers by carlo, whuber and comments already say a lot of this - is to re-express the logit expression as $\exp(\beta_1 (\gamma+X))\over 1+\exp(\beta_1(\gamma+X))$, where $\gamma={\beta_0\over \beta_1}$. Doing so, you can maximize the likelihood $$ \max_{\beta_1,\gamma} E\left [\mathbf{1}(X>c)\...


0

No, unequal sample sizes are not a problem. I don't know if there is a canonical citation for the unbalanced data claim, but this paper mentions it, and it comes from one of the authors of the lme4 package. Pinheiro, J. C. (2014). Linear mixed effects models for longitudinal data. Wiley StatsRef: Statistics Reference Online. NA Just check the regular ...


2

Pearson Residuals As tosonb1 points out, "The Pearson residual is the difference between the observed and estimated probabilities divided by the binomial standard deviation of the estimated probability". I just wanted to mention that Pearson residual is mostly useful with grouped data i.e, say, there are $n_i$ trials at setting i of the explanatory ...


1

This link may be helpful for you: https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-how-do-i-interpret-odds-ratios-in-logistic-regression/ Typically, odds ratios are interpreted as scalars that represent the increased/decreased likelihood of association to the response variable. For example, if this GLM was from a logistic family, we could say that,...


1

It's saying that the odds of females graduating is 4.31 higher than the odds of males, holding constant disorder. That doesn't mean females are 4.31 times more likely to graduate, but it does mean females are more likely to graduate. You can think of it as that among females, the ratio of the number who graduate to the number who don't graduate is 4.31 times ...


0

I think you want to first use R's predict function to first generate Prediction intervals or Confidence intervals around the independent values. Read this link for more info. R prediction info The next step is to use those intervals to generate possible values with a random number generator scaled appropriately for those values.


1

One way to look at the data: x <- c(2, 4, 5, 6, 7, 8, 9, 9, 4, 5) y <- c(0, 1, 0, 0, 0, 1, 1, 0, 1, 0) stripchart(x ~ y, meth="stack", ylim=c(.5,2.5), pch=19) abline(h=1.5, col="green2") From the stripcharts, we can see that the data values for the two groups are not much different. If integer values given in x are rounded continuous variables, ...


0

A simple method is bootstrapping. Estimate regression, get the residuals. Get an empirical cdf of x. Then use inverse of it to generate new x set. Plug new x into the regression model and add bootstrapped residuals That was if you think x causes y. If there’s no causality then it’s easier. Get eCDFs of x and y. Then estimate correlation of eCDF outputs ...


Top 50 recent answers are included