Hot answers tagged r
7
From the look of your link below, you obtained this plot:
It looks like you're just plotting y vs its index in Libre Office Calc.
Let's look at the difference - first in R:
R gives (using plot(x,y)):
But if you plot y vs its index and plot lines (plot(y, type="l", col=4)), you get:
If you don't tell Libre Office Calc to do a scatter plot, by ...
7
You got the same results because you didn't specify a special variance or correlation structure in the gls function. Without such options, a GLS behaves like a OLS. The advantage of a GLS model over a normal regression is the ability to specify a correlation structure (option correlation) or allowing the residual variance to differ (option weights). Let me ...
6
The results don't necessarily show that a gamma distribution is a good fit. You might as well as assume that a normal distribution must fit well because you have estimated a mean and a standard deviation.
Your graphical assessment is already showing that a good (if not necessarily the best) gamma shape for these data is systematically different from your ...
5
The mode is a bit of a red herring. Here is a very simple solution to this problem that circumvent the need to define the mode precisely. I'm surprised it has not been proposed earlier. The constraint on the mode can be easily satisfied by drawing samples from a symmetric distribution and scaling them suitably:
...
5
"Dichotomous Predictor Variables", there are two ways to code dichotomous predictors: using the contrast 0,1 or the contrast 1,-1.
This is factually wrong. There is no limit to the number of ways they can be coded. Those two are merely the most common (indeed between them, almost ubiquitous), and probably the easiest to deal with.
I kind of ...
4
I think you're specifying a great deal of random effects which are better modeled as fixed effects since you'll directly estimate the odds ratio for their association with outcome. For instance, the question of whether partisanship modifies the relationship of discount and product is simply done as a test of interaction:
m1 <- glmer(chose ~ partisanship ...
3
Why would a major drug company even want to convert to R from SAS? SAS costs millions but it is nothing to a drug company. However, converting all the stable reporting systems from SAS to R would cost 50-100 times more.
SAS has phenomenal support system: every time I needed help they were able to provide it within few hours.
And what exactly does R have ...
3
i) First, a recommendation:
Use pchisq( -2*sum(log(p-values)), df, lower.tail=FALSE) instead of 1- ... - you're likely to end up with more accuracy for small p-values. To see that they're sometimes going to give different results, try this:
x=70;c(1-pchisq(x,1),pchisq(x,1,lower.tail=FALSE))
ii) Yes, it's one-sided. Small values of the chi-square ...
3
You need to use the matrix \$loadings, not \$rotmat;
x <- matrix(rnorm(600),60,10)
prc <- prcomp(x,center=TRUE,scale=TRUE)
varimax7 <- varimax(prc$rotation[,1:7])
newData <- x %*% varimax7$loadings
The matrix $rotmat is the orthogonal matrix that produces the new loadings from the unrotated ones.
3
Here is a trick I've used before, although I don't think I've ever published it. If x is your monthly time series, then you can construct annual totals as follows.
y <- filter(x,rep(1,12), sides=1) # Total of previous 12 months
To get the forecasts of the annual totals:
library(forecast)
fit <- auto.arima(y)
forecast(fit,h=12)
The last forecast ...
2
Your error is in the following line:
x<-(2*ret*pnorm(a)) + (2 * exp(b) * (ret * (sd^2)) * pnorm(c))
Change the multiplication (ret * (sd^2)) to an addition (ret + (sd^2)) and you're doing fine: the function returns the value $0.713$ which matches the value of the function in the graph on page 4.
2
It sounds like a two component gaussian mixture model would be exactly what you're looking for.
I don't tend to use R, but GMMs are so commonly used I can't imagine it they're not implemented somewhere. If you really can't find them anywhere then they're not too hard to implement yourself. Wikipedia has the exact expectation maximisation update steps ...
2
My first observation is that you did not lag the inputs relative to the closing price and that is why you observed such good fit. The SMA (simple moving average) uses the closing price in its calculation and the high low range encompasses the closing price, so using them to predict the closing price imparts a look ahead bias. My opinion is that if you are ...
2
If you don't have them already, start by setting up some predictors, $x_1$, $x_2$, ...
Choose the population ('true') coefficients of your predictors, the $\beta_i$'s, including $\beta_0$, the intercept.
Choose the error variance, $\sigma^2$
generate the error term, $\varepsilon$, as an independent random normal vector, with mean 0 and variance $\sigma^2$
...
2
Make a glm fit to the marginals with Poisson errors (yielding a log-linear model) and then use predict on expand.grid data.frame from the the row and column values based of the second sample. (There's no particular advantage that I can see in using IPF to estimate a log-linear model of this sort.)
require(reshape2)
Loading required package: reshape2
> ...
2
You forgot the intercept and did not calculate the residual variance correctly.
Here is how I would do it:
#get the design matrix
Designmat <- model.matrix(eval(eval(m$call)[-2]), as.data.frame(X))
#this is the same as cbind(1,X) for your model
#calculate the standard errors
predvar <- diag(Designmat %*% vcov(m) %*% t(Designmat))
SE <- ...
1
It depends on what you want to use this for. I can easily imagine situations in capacity planning where you would be most interested in extreme occurrences, as these peak events are what strains capacity most. If that is the case, then your tail behaviour would be a problem. I can also imagine other situations where the system is somewhat flexible so that ...
1
I think your reasoning is sound. The interaction basically tests if your lines are parallel (i.e. if the slopes differ for plants that did have a LWR on the site vs. plants without LWR). In your graph, the lines are almost perfectly parallel which supports the conclusion from your regression containing the interaction. What I might add: The adjusted ...
1
I have to agree with Glen.
It is axiomatic in control system's engineering that there is no such thing as "best" without a measure of goodness.
Some (weak) examples of candidate bests include:
Best = robust indicator of central tendency
Best = robust indicator of variation around central tendency
Best = fastest to compute
Personally, when trying to ...
1
auto.arima can select the order of differencing automatically, and make appropriate forecasts taking account of these differences. It can also include regressors via the xreg argument, and select the appropriate model order taking account of the regressors. If you include regressors, these will be differenced along with the response variable as part of the ...
1
The Standard error is an estimate of the variance of the strength of the effect, or the strength of the relationship between each causal variable and the predicted variable. If it's high, then the effect size will have to be stronger for us to be able to be sure that it's a real effect, and not just an artefact of randomness.
The t-statistic is an estimate ...
1
It sounds like you need a decent basic statistics text that covers at least basic location tests, simple regression and multiple regression.
Std. Error,t value and Pr.
i) Std. Error is the standard deviation of the sampling distribution of the estimate of the coefficient under the standard regression assumptions. Such standard deviations are called ...
1
You can use a likelihood ratio test (LRT) to test whether a random effect is significant. First, fit the random effects model. Then fit the model without the random effect. Extract and store the log-likelihood for each model using logLik and calculate the twice difference between the log-likelhood of the mixed effects model and the normal model. Use a ...
1
Here is an approximate numerical answer. It can easily be made more precise.
Let $\{U,V,W\} = {X,Y,Z}/(X+Y+Z)$, where $X,Y,Z$ are i.i.d. with a trapezoidal density on $[0,1]$:
$f(x)=1+a-2ax.$ $U,V,W$ will have identical marginals.
Given a numeric 'a', I used Mathematica to get the cdf of $U$:
F[u_] = Assuming[0 < u < 1, Simplify@Integrate[
Boole[x ...
1
If you fit an arima model with external regressions, you MUST provide newxreg to the predictions function. This goes for arima, auto.arima, and Arima. You cannot provide external regressors to ets.
xreg should contain the exogenous variables for the training set, and newxreg should contain those variables for the prediction set. If the nexreg contains ...
1
There are two main ways to do binomial glms; you mention one (we'll call that the first way).
There's also a third way I'll cover at the end, which by the sound of it is perhaps what you're after.
The second way
This is clearly described in the help on GLMs (?glm), which is to supply the outcome for each individual bernoulli (0-1) trial. The help ...
1
I think that your approach is correct. Model m1 specifies a separate intercept for each subject. Model m2 adds a separate slope for each subject. Your slope is across days as subjects only participate in one treatment group. If you write model m2 as follows it's more obvious that you model a separate intercept and slope for each subject
m2 <- lmer(Obs ~ ...
1
Another possible solution is as follows:
# Create an "acf" object called z
z <- acf(dummy)
# Check class of the object
class(z)
# View attributes of the "acf" object
attributes(z)
# Use "acf" attribute to view the first 13 elements (1 = lag at 0)
z$acf[1:13]
# Get rid of the first element (i.e. lag 0)
z$acf[2:13]
# Plot the autocorrelation function ...
Only top voted, non community-wiki answers of a minimum length are eligible
