New answers tagged

1

I agree with @whuber that results are clear without a formal test. However, it you are interested in checking the validity of your testing procedure in R, some results from R may be helpful. Data: liv = c(654, 57); die = c(1438, 52) tot = liv + die prop = liv/tot prop [1] 0.3126195 0.5229358 Test of binomial proportions. You can do a test to see if the two ...


1

Decision trees (generally) choose splits to minimize impurity within the nodes. If a split can change a 75% "no" node into two children with 99% no and 51% no (respectively), then it is reasonably likely to do so, even if the final predicted (majority) class in both children is "no".


0

No. What you did above is simulated $(X,Y,Z)$ where $X$, $Y$ and $Z$ are independent. The sample covariances will simply converge to 0 as $n$ converges to infinity. To simulate any random variables jointly, you have to have an idea how they are related


0

If you want the intercept to be the mean of the levels, then effects coding will give you an unweighted mean as your intercept. It will have a standard errors. Moreover, you can now use the standard errors produced by your beta estimates to calculate confidence intervals to see if they cross the unweighted grand mean. It can be done in R using the contrasts ...


0

my answer will not be specific to R but I think you can find the way of doing this by looking elsewhere. UPDATE: you say you just want to see if they learn at all. Then you don't need to use $k$, but check if $r_0 \neq r_\inf$. For this actually, you don't need to fit a logistic. I think a simpler approach should suffice. In the end, the best way would be to ...


0

Have a look this : http://gradientdescending.com/state-space-models-for-time-series-analysis-and-the-dlm-package/ it says "The main output here is s. Again there are 13 columns, the first 2 from the linear trend model and the last columns from the seasonal component."


0

I don't know if this is the best method, but one method I've used in the past is to train an Autoencoder for this task. How this works Train an Autoencoder (AE) to model your data $X$. What you essentially want from this step is for the AE to learn the distribution of $X$. We will use this to examine if every new sample belongs to the same distribution or ...


1

It is difficult to analyze this implementation with pomp as it is not a standard package (and I can not get it installed for the version of R on my computer, alternatively you could use deSolve to compute the ODE, as this example is a deterministic case). However, it is not neccesary to run your code as it already looks peculiar based on a quick scan. It ...


1

Good question. Let $(y,x)$ be your test data with $y = x^T \beta^* + \varepsilon$ and assume that $x$ is random as well (independent of $\varepsilon$). To simplify, let $$\delta = \beta^* - \widehat \beta.$$ Then, your mean-squared prediction error is \begin{align*} \mathbb E [(y - x^T \widehat \beta)^2 \mid \mathcal T] &= \mathbb E [(x^T \beta^* - x^T \...


0

lme fits linear mixed models, not non-linear ones. For certain types of non-linear mixed models, you can use nlme. If area is nested in the levels of id, you can use the code model = lme(dv ~ intervention, random=~1|id/area, data=Data, method="REML") to fit the nested mixed effects model. Edits: If you just specify random ...


0

In gamlss BEINF0(mu, sigma,nu) Y has a point probability at zero p_0 = P(Y=0) = nu/(1+nu) and has a BE(mu, sigma) distribution on (0,1) with probability (1-p0). Hence E(Y) = mu/(1+nu) On the predictor scale eta = logit(mu) you will get parallel lines against x. But on the mu scale you will get S shaped curves shifted horizontally, (but only shown for x from ...


1

I just ran across this and thought I should try to clarify a couple of points. It it true that the the QR solution to this "regression through the origin" problem is one of the y_i/x_i values, but it is a weighted quantile with weights |x_i|. Floyd and Rivest's SELECT is somewhat quicker than Hoare's quickselect for ordinary sample quantiles, an ...


0

Outliers are observations that are drawn from a different distribution (process) than the majority of the data. If you are removing points that are from a distinct process than the one that you are trying to model, then it is generally acceptable to remove them or model them somehow if that is feasible. If you can justify your decision to remove those points ...


1

With respect to the last issue, you treat them as comparisons. R has made the "none" the base category, so we interpret it as follows: For animals with no treatment, as compared to those with low treatment, there is a -239.50 change in value, on average. Repeat the comparison with None vs other condition. You can change what becomes the base "...


4

The way to represent models can vary a lot with the application, and I'm not sure if there is any name rather than "vector" form or "observation level" form. We need to know if we are numbering the individuals within groups or not. This only makes sense if it was a designed study with blocks (one categorical variable) and a treatment ...


0

In your example, everyone gained weight and you tested the hypothesis that everyone lost weight (alternative = "less"), so it VERY VERY UNLIKELY that the participants lost weight, so $p=1$. Try instead "greater" or "two-sided". > wilcox.test(weight ~ group, data = my_data, paired = TRUE, + alternative = "...


2

You are indeed simulating the data backward. You should simulate the potential outcomes first and then the outcomes from them. The "consistency equation" $Y = AY_1 + (1-A)Y_0$ describes the data-generating process for $Y$, the observed outcomes. The potential outcomes "occur" before the treatment is assigned in the formulation of causal ...


0

Model Diagnostics With a parametric survival regression model you do have residuals. You just might have to work a bit harder to get them than you do with ordinary least squares. Chapter 18 of Frank Harrell's course notes and of his book on regression modeling strategies goes into detail about diagnostics for parametric survival models. A parametric model ...


0

As the comment from @whuber notes, you can readily test whether the order of variable entry matters in terms of the model itself. (It won't.) There might, however, be an issue if you try to use some forms of anova() analysis on the model. I'd recommend using an anova() function that doesn't depend on the order of variable entry, like the that in the R rms ...


0

First of all, while the advice you received in your course reflects the general sentiment of the stats community, there is nothing wrong with a transformation, and due to simplicity, I would prefer this solution to a GLM in your case. So, before you move to complicated GLM / GLS solutions, try e.g. if https://www.rdocumentation.org/packages/forecast/versions/...


1

As mentioned in the comment, the first n days in normalized days are disturbing to me, so I'd like to stay with discrete days. The data then becomes: data <- data.frame(day = 1:7, y = c(-1.77840188804283, -0.766238080609954, -0.174160208674574, 0.245925726822919, 0.571769693063722, 0.838003598758299, 1.06310115868241)) ...


2

You can, of course, use any numerical clustering (k-means, hierarchical clustering, spectral clustering...) on the projected data (the points). But why should it be better than clustering in the original, high dimensional feature space?


0

I do not know why day is negative in your data and not in sequence from 1, 2, ... I guess those values are obtained after transformation. For simplicity, assume that we are day in whole number sequence. Suppose that you want to compute $sum(y)$ up to day 10, i.e., $sum_{y10}=y_1 + y_2 + ... + y_{10}$. You need to know $y_8$, $y_9$, and $y_10$. By inserting $...


0

A few suggestions with respect to working with crime data. First, your plot is investigating the relationship between unemployment rates and total crime counts across many geographic regions. It is unclear what a "region" represents, but it appears they are rather large aerial units (e.g., counties, states, countries, etc.). Some regions experience ...


2

There are now simpler ways to do this. For example, the preprocess function of the caret package library(caret) preproc <- preProcess(data, method = c("center", "scale") scaled.new <- predict(preproc, newdata = new) or scale_by in the standardize package or using the receipes package library(recipes); library(dplyr) rec <- ...


0

The facebook prophet package supports multiple seasonalities. Yearly, weekly and daily seasonalities are built-in but custom seasonalities can be specified. Here is a custom monthly seasonality: df <- ... # data to build model on or decompose future <- ... # data to make forecasts on m <- prophet(weekly.seasonality=FALSE) m <- ...


1

The term is estimated marginal means. In R there is a package emmeans that does exactly this. Once you have the estimated marginal means then you can further use those to calculate different contrasts than the contrasts that were originally specified in the model.


0

You need to put the bs = 'cr' bit inside each of the s( ) terms; at the moment it is being ignored as it isn't an argument to gam() but is getting passed along via ... and silently ignored. If you run summary(modeltest) you'll see the parametric terms listed for the factors. To see them in the plot, you need to add all.terms = TRUE to also have them plotted ...


0

How to interpret the random effect?? Generally we don't usually interpret the random effects. We specify random effects to account for non-independence of observations with repeated measures. Random slopes are used to allow for the "effect" of a variable to differ between levels of a grouping variable (in this case ID). According to the question ...


1

I'm confused as to why you'd fit 3 day as a cyclic smoother when you only have three months of data? The point to cyclic smooths is that there is a natural point that is both small and large that are essentially the same thing because the covariate space wraps around on itself; think wind direction, aspect, and day of year where you have observations from ...


1

The heart of this question seems to be how to treat the trial variable. It is clear that random intercepts are needed for subject to account for the repeated measures. By specifying trial as a random slope for subjects we are allowing the "effect" of trial to vary for each subject. Since trial is absent as a fixed effect in your models, this means ...


0

I have a very similar related question. Not an expert, but based on what I've researched so far and what is highlighted on this page, here are some points that may be helpful in answer to your last question Is it possible that the analysis could be valid despite the non-proportional hazards? As others have mentioned, large sample size may be a factor ...


0

I figured it out, so for those who are interested: the model.matrix command does not include the offset, so it needed to be added manually. Once the offset is taken into account, the predictions of predict, type = "response" and the manual approach provided by Zuur coincide.


0

Your diagnosis is probably correct, and the big standard errors is because the loglikelihood function in this case is far from quadratic. Confidence intervals based on likelihood profiling might be better, as discussed at Why is there a difference between manually calculating a logistic regression 95% confidence interval, and using the confint() function in ...


1

Gene is a factor with three levels. You can confirm this by checking levels(factor(data$gene)). The reference level is rolled into the intercept term, and the coefficients corresponding to the other two levels are displayed next to gene1 and gene2. The model fit looks unreliable due to the exceeding high standard error on gene2. I suggest taking a look at ...


0

Although this question already has an accepted answer, I would like to address the yet unadressed problem of a possible non-convergence of k-means. This can indeed occur due to ties and cyclical assignemnt of ties to different clusters. To avoid this issue, you must always use the same deterministic tie breaking, e.g., always assign ties to the cluster with ...


2

$n\sigma^2$ is the variance of the final value, i.e. $y_n$, however you calculate the variance based on the historical values of $y_t$, that's one reason you get lower values. If you change your code to the following, you'll obtain values closer to $250$, estimated over $180$ realisations of $y_t$. Of course, to build an histogram of variances, you'll need ...


0

I think what you see is mostly caused by the forecast $\hat y_t$ being a function of $y_{t-1}$ only which (sometimes but not always, see below) makes $\hat y_t$ more similar to $y_{t-1}$ than $y_t$. If we for simplicity ignore the uncertainty in the parameters, the one step ahead predicted values based on all data up to time $t-1$ for an AR(1) model $y_t=\...


1

If anyone looks at this post in the future, I have collected answers from experts, so I'm all set. The proportion mediated is negative if the direct and indirect effects have opposite directions. It makes more sense to discuss the proportion mediated if the direct is the same for both. If the outcome is binary, the ACME and ADE are absolute percent changes, ...


0

You can try Y <- as.matrix(Fish[,2:7],nrow=nrow(Fish),ncol=6)%*%as.matrix(FishPC,nrow=6,ncol=2) However, be aware of how the PCs were computed. PCA is often performed on data that has been centered and sometimes scaled. If that is the case for your PCs, then you should also center and/or scale Fish as appropriate.


0

geom_line by default produces smoothed lines within plotted points. I'd like to know if you tried the following, which can probably solve this for you: plot + stat_smooth(method = 'lm', se = F) If you don't specify that your method derives from a linear model (method = 'lm'), stat_smooth function will produce the same connect the dots line as geom_line .


0

I actually found the problem. The higher moments IV approach requires the endogenous regressors to be non-normal. In my example above, the logs of $Z_1$ and $Z_2$ are almost normal. If I change rlnorm(N,0,1) + 10 in the data generating process to e.g. rlnorm(N,0,1.5) + 10 also the logs of $Z_1$ and $Z_2$ remain non-normal and the IVs for the log-log model ...


1

For example, If I have these coefficients, How can I manually predict (by hand, not using R) the outcome (1, 2 or 3) for x1 = 26 and x2 = 55 ? You cannot do it by hand. The conditional probabilities are given by 2D integrals with the integrand being a multivariate normal distribution density function. There is no closed form solution. See the vignette ...


0

How many levels are there for each IV? It is hard to interpret this without fully understanding the structure of your data. I will still try even though I don't know the structure of your data. If your IVs are categorical, (Intercept) will be the first level of the factor. If your IVs are continuous/numerical, (Intercept) will be the expected (mean) value of ...


1

From documentation for lm Non-NULL weights can be used to indicate that different observations have different variances (with the values in weights being inversely proportional to the variances); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations (including the case ...


2

Assume the actions of these individuals are independent. (If you don't, you need to specify--mathematically--how they are interdependent.) Assume that the chance any individual leaves during any time period does not change over time. (This can be modified, but then you will have to specify how these chances vary over time.) These assumptions imply that the ...


0

In 2017, Dahlin & Schön published a really good article with R scripts for both Kalman and Particle filters, as well as Metropolis-Hasting algorithm: https://arxiv.org/pdf/1511.01707.pdf It is not the most efficient code (its in R only), but very well made.


0

The same problem occurs in image processing with mean filtering or other convolution operations, and it is solved in different ways by different "border treatments", which can be, e.g., repeat, reflect, zero, or wrap around. As the igraph library function running_mean does not support a "border treatment" option, you must implement it ...


0

I suggest you read Discovering Statistics Using R by Andy Field. He explains gives information about how to perform various robust ANCOVA in R step by step and it is really easy to follow.


0

A simple approach would be to draw $X\sim U[0,1]$ and $Y\sim N(x^2,x)$, i.e., the variance of $Y$ is $x$. Then the $0.25$ quantile of $Y$ would be $z_{0.25}\sqrt{x}$. Here are simulated data, with the expectation as a black and the quantile as a red line: R code: set.seed(1) xx <- runif(1e3) yy <- rnorm(length(xx),xx^2,sqrt(xx)) plot(xx,yy,pch=19,las=...


Top 50 recent answers are included