title: ANOVA on percentages

text: **Please read edit 3 first** I am trying to find out the significant factors in a dataset of percentages, a sample of which are below. The difficulty is that the data violates the assumptions of ANOVA, and most of the data are fairly close to 100%. Please note that glm + binary would not work: the samples used to calculate each proportion are not independent. I do have access to the denominator, if that helps. Any direction where to start? I've read quite a few things here and elsewhere (notably trying to use some transformation such as arcsin, etc...) and also some other methods I never heard of ("contingency table approaches"). As in a "textbook ANOVA" I would like to know which factors are significant, and how much of the variability they explain.               data = c(0.79,0.98,0.95,0.95,1,0.98,0.99,0.97,0.99,0.99,0.98,0.99,0.99,0.94,0.94,       0.86,0.84,0.86,0.97,0.96,0.53,0.87,0.97,0.81,0.99,1,0.99,0.87,0.98,0.97,0.93,0.8,       0.7,0.94,0.89,0.98,0.89,0.98,0.96,0.98)      **Edit:** sorry for the lack of clarity. Here's how my percentages are calculated: I basically throw a number of particles (known only a posteriori) in a funnel and count how many make it through/how many remain stuck in the funnel. The percentage is the ratio of the particles having made it through divided by the total number of particles.If a particle which comes at the beginning of the trial remains stuck, the odds of a subsequent particle to remain stuck are higher. As such, I don't think I can apply a generalized linear model, specifying binomial as familly (in R I mean). But again, my statistical insights may be utterly wrong. **Edit2:** regarding independance, I guess my comment was misleading. Each sample in the vector above is independant from the others. However, as I explain in the edit above, the samples used to calculate each percentage are not. **Edit3:** Well, I reckon I probably made a mess of my question, think it may help if I rephrase the problem and show some data (a fraction of my whole dataset). Besides, I have progressed a bit, hopefully in the right direction. I do not know if I should do a full edit of my original question, or even abandon and start anew, I'll add this as an edit for the interest of history (let me know if I should do differently). My response variable is the percentage of particles having made it through a funnel (the total number of particles is different for each percentage). If a particle at the beginning of the trial remains stuck, the odds of a subsequent particle to remain stuck are higher. The (potential) explanatory variables are 1) the type of particles, 2) the funnel type, 3) the funnel position and 4) the total number of particles. In a first stage, I want to find which of these actually impact the response variable and by how much. In a second stage, I would like to use the current dataset as a reference for the analysis of future samples. Precisely, I would like to know if the percentage of particles having made through the funnel is significantly different from the reference dataset and by how much. Plotting the data indicates me that each population may have a different mean and a different variance.               # Libraries ####     library(ggplot2)     library(betareg)     library(lmtest)          #Create data####     df5 = structure(list(type = c("Type1", "Type1", "Type1", "Type2", "Type2",                                    "Type2", "Type2", "Type2", "Type2", "Type3", "Type3", "Type3",                                    "Type1", "Type1", "Type1", "Type3", "Type3", "Type3", "Type2",                                    "Type2", "Type2", "Type1", "Type1", "Type1", "Type1", "Type1",                                    "Type1", "Type3", "Type3", "Type3", "Type2", "Type2", "Type2",                                    "Type1", "Type1", "Type1", "Type3", "Type3", "Type3", "Type1",                                    "Type1", "Type1", "Type3", "Type3", "Type3", "Type2", "Type2",                                    "Type2", "Type3", "Type3", "Type3", "Type2", "Type2", "Type2"),                           funnelType = c("fType1", "fType1", "fType1", "fType1", "fType1",                                          "fType1", "fType2", "fType2", "fType2", "fType1", "fType1", "fType1",                                          "fType2", "fType2", "fType2", "fType1", "fType1", "fType1", "fType1",                                          "fType1", "fType1", "fType1", "fType1", "fType1", "fType1", "fType1",                                          "fType1", "fType2", "fType2", "fType2", "fType1", "fType1", "fType1",                                          "fType2", "fType2", "fType2", "fType1", "fType1", "fType1", "fType2",                                          "fType2", "fType2", "fType2", "fType2", "fType2", "fType2", "fType2",                                          "fType2", "fType2", "fType2", "fType2", "fType2", "fType2", "fType2"),                           position = c("a", "b", "c", "a", "b", "c", "a", "b", "c",                                        "a", "b", "c", "a", "b", "c", "a", "b", "c", "a", "b", "c", "a",                                        "b", "c", "a", "b", "c", "a", "b", "c", "a", "b", "c", "a", "b",                                        "c", "a", "b", "c", "a", "b", "c", "a", "b", "c", "a", "b", "c",                                        "a", "b", "c", "a", "b", "c"),                           total = c(420L, 293L, 324L, 549L,                                     527L, 603L, 533L, 571L, 438L, 496L, 534L, 604L, 489L, 360L, 383L,                                     524L, 560L, 606L, 493L, 513L, 572L, 530L, 527L, 543L, 616L, 471L,                                     554L, 392L, 530L, 443L, 561L, 529L, 599L, 529L, 481L, 521L, 621L,                                     567L, 609L, 447L, 398L, 462L, 528L, 574L, 522L, 654L, 531L, 556L,                                     642L, 569L, 684L, 372L, 540L, 345L),                           percentage = c(0.98808, 0.95569, 0.9784, 0.57741, 0.81017, 0.82919, 0.93809, 0.94041,                                          0.94744, 0.93352, 0.98129, 0.99006, 0.97339, 0.94728, 0.98695,                                          0.9542, 0.84989, 0.92574, 0.79116, 0.92782, 0.98077, 0.96605,                                          0.97155, 0.96503, 0.97076, 0.95756, 0.99097, 0.96429, 0.9434,                                          0.99097, 0.78253, 0.94518, 0.97664, 0.99056, 0.96261, 0.98273,                                          0.88402, 0.96824, 0.92118, 0.95524, 0.97991, 0.9762, 0.83144,                                          0.95643, 0.98085, 0.95107, 0.94162, 0.98741, 0.83635, 0.94372,                                          0.98099, 0.88447, 0.94817, 0.95365)),                      .Names = c("type", "funnelType","position", "total", "percentage"), class = "data.frame",                      row.names = c(2L,                                    3L, 4L, 18L, 19L, 20L, 50L, 51L, 52L, 66L, 67L, 68L, 98L, 99L,                                    100L, 114L, 115L, 116L, 130L, 131L, 132L, 146L, 147L, 148L, 162L,                                    163L, 164L, 194L, 195L, 196L, 210L, 211L, 212L, 226L, 227L, 228L,                                    258L, 259L, 260L, 306L, 307L, 308L, 322L, 323L, 324L, 354L, 355L,                                    356L, 370L, 371L, 372L, 386L, 387L, 388L))          ggplot(df5,aes(x=total,y=percentage))+geom_boxplot(outlier.shape = NA)+geom_jitter(aes(col=position))+facet_grid(type~funnelType)      It does seem that I am in a situation comparable the dyslexic example in http://psycnet.apa.org/journals/met/11/1/54/. I consequently think that I should be able to apply the same method to analyse my data. To do so I use the package betareg in R (no interactions for the sake of the example).               full <- betareg(percentage ~ type+funnelType+position+total|type+funnelType+position+total,data = df5)     summary(full)      This yields the results that all explanatory variables but "total" have a significant effect on the mean. From what I know explanatory variables are likely to be significant in large data sets (I have a relatively large number of observastions in my full dataset). How do I judge the magnitude of the effect of each factor? What should I do if two of the explanatory variables are correlated? In my second stage I want to compare a new sample to the subpopulation in my dataset that has the same predictor values. How should I perform this? Unequal variance t-test?


title: Best approach to classifying 3-point trajectories?

text: I have a sample of about 300 subjects who have been measured at 3 different times (morning, afternoon, evening). The variable of interest can be assumed to be approximately normal. It appears that most subjects have an increase between between morning and afternoon, followed by a decrease from afternoon to evening. Some however show an opposite pattern (devrease->increase), while yet others remain approximately the same. What I am interested in is clustering, or classifying the subjects according to their trajectories. After a bit of googling, I have discovered GMM. I don't really understand what goes on behind the scenes, but it looks like the classification is done based on straight line fits when we only have 3 points. This seems highly inappropriate to me since straight line does not capture the increase followed by a decrease type of behaviour. Another thing is, people seem to be using Mplus package, which I am not familiar with and would rather avoid buying and learning (I am very comfortable with R and Matlab).


title: SVD on a 65 million by 3.4 million sparse matrix

text: What is the best way compute SVD on a very large positive matrix (65M x 3.4M) where data is extremely sparse?   **sparse** : less than 0.1% of the matrix is non zero.   **best way** :   * will fit into memory (I know that online methods exists)   * will be computed in a reasonable time: 3,4 days   * will be accurate enough however accuracy is not my main concern and I would like to be able to control how much resources I put into it. It would be great to have a Haskell, Python, C# etc. library which implements it. I am not using mathlab or R but if necessary I can go with R.


title: Trend Analysis of feature importance over time in R

text: I'm running an experiment on a Streaming Classification Model (an Online Random Forest) that I've created. If that is a completely foreign concept to you here is a presentation I did on it recently: http://www.slideshare.net/CvilleDataScience/online-random-forest- in-10-minutes. Essentially to test my model's responsiveness to learning new signals, I've simulated a data-set that contains two signals. The model is trained on one signal for a period of time and then I start training on the other signal (in an immediate step change fashion). Incrementally throughout this training process I ask the Online Random Forest for reports on feature importance. The result is a time series of feature importance for the different features/predictors over time. Here is a visualization of this time series with the black line representing when I changed the signal: ![Feature Importance Over Time](http://i.stack.imgur.com/nQ8Iz.png) The experiment I'm running changes a number of the different model parameters and generates a bunch of these time series matrices. What I'd like to measures is the time to stabilization for each of the matrices. I'm interested in any approaches for doing this. What I'm tried so far is reducing the entire matrix time series down to a single time series based on:   1. The subsequent Euclidean Distance between feature importance vectors   2. The angle between feature importance vectors After getting this single time series I've been running it through quality control processes in R using the qcc library. Using that library I measure how long the process stays in the "out of bounds" state. This works pretty well, but isn't as robust as I'd like, meaning sometimes the out of bounds violations don't happen sequentially even thought "I think" they should. The detailed part of my question is this. Given the time series for one of the features that I know I changed from "unimportant" to "important", how could I dynamically fit a line from where the feature importance goes from the stable "unimportant" state to the stable "important" state. Essentially how could I break apart the signal into 3 components? How could I do this dynamically given I don't know "when" the process reaches a new stabilization point? Here is an example of this line (the lower plot is the line smoothed with the smooth function in R): ![enter image description here](http://i.stack.imgur.com/UJpy8.png) Any references to literature or R packages is much appreciated!!


title: What's a good measure of spread near a threshold?

text: I am choosing summary statistics for percentage data that is clustered near 100%. Each subject takes the same test twice. Each test consists of at least 40 questions. I want to be able to predict which subjects will improve more. I have a hypothesis that subjects with higher variability of performance in the first test will improve more. To test this, I computed the variance of each subject's per-question scores, and indeed those with higher variance improved more. But the problem is that the presence of the threshold introduces an artefact: subjects with near-100% performance to start with can't improve much and also may have lower variance simply because the upper tail of question scores is limited by the maximum score threshold. Inter-quartile range would suffer from the same problem. I will try "one-sided variance" (using only question scores lower than the mean), and "lower quartile to median range". Any other ideas, and does anyone know of a use in the literature of something like this?


title: How to produce a vector based on two probability distributions?

text: I am working on distribution on plants and I need to enter the chances that a plant starts reproducing at a certain age in a vector. Only this has to be done for more generations. Now, Imagine I have a plant that starts reproducing at least at 80 years and at last at 81 years. I assume that the chance for a plant to reach its reproducible age is equally distributed, so the chance a plant starts reproducing at 80 is 0.5 and the chance a plant starts reproducing at 81 is also 0.5. However, I want to know these chances for 2 generations of a plant. So I make a tree diagram, where the chance that 2 generations last 160 years (2*80) is 0.25 (0.5*0.5). The chance that 2 generations last 161 years (80+81 or 81+80) is 0.5 (2*(0.5*0.5)). The chance that 2 generations last 162 years is also 0.25 (0.5*0.5). The vector that I need here is a cumulative vector of chances: c(0.25, 0.75, 1) However I have a plant that starts reproducing at least at 80 years and at last at 149 years. I assume again that the chances for a plant to reach its reproducible age is equally distributed, so that is 1/70 a plant starts reproducing at 80, 1/70 at 81 etc. I want to know these chances for 20 generations. (You can of course make a tree diagram of this as well, but that will be one enormous tree diagram). This will become a vector of 1400 values, where the first value represents the chance that all generations start reproducing young and that is in (80*20)= 1600 years and the last value represents the chance that all genereations start reproducing at ist oldest age (150*20)= 3000 years. In this vector each value stands for the chance in one year (3000-1600=1400). **So here's my question: Is there a way in R to calculate such a cumulative vector of chances?**


title: Inverse correlation of time series data

text: I have 3 hourly power generation data for around 600 locations for a year. (i.e. 8 data per day for 365 days for each location.) I want to find out a way where out of this 600 locations, I can say suppose choose suppose 10 locations which are producing power at different time periods. suppose:               Location:   Hour 0, Hour 3, Hour 6, Hour 9, Hour 12,  Hour 15             A          30       00      50      70       00      20            B          50       20      70     100       00      40        C          00      100      20      00       40      30        D          20       15      10      00       40      30      Here, A and B are highly correlated where A and C are inversely correlated. So, is there a way where I can distinguish the locations which are producing powers at different time periods so that, by accommodating these few location, the overall power generation can be sort of stable? By sort of stable, i mean something like this which I mentioned in the comment later on to clarify my question. Since A and C are inversely co-related, so if i choose power from A and C out of four locations, I would get some power all the time (always above 30 in this scenario). But if i take A and B, I would get very high amount of power when I have power (like 70+100=170 in hour 9), and I wont get any power at Hour 12 (since A and B produces 00 power in Hour 12). I want to avoid choosing A and B, and I want to choose A and C. my objective is to choose 10 locations out of 600 where by adding powers of each hour of 10 locations, total power for each hour should be above a certain threshold, like 30.


title: Autocorrelation of daily time-series using acf and pacf plots

text: I have a time series dataset with daily entries over the past 20 years. I would like to see whether the series is autocorrelated. To do so I have done the following>               dataset=read.csv(file="https://dl.dropbox.com/u/22681355/dataset.csv", head=TRUE)     par(mfrow=c(2,1))     dataset=diff(dataset)     acf(dataset)     pacf(dataset)      From the plots it seems that there is no autocorrelation, however, I'm not sure whether the acf and pacf plots are suitable for my dataset. What is a suitable test for stock data ?


title: Visualising plausible values in multiple imputation

text: I have been checking density plots to get a feel for the plausibility of values that have been imputed using the `mice` package in R. I would be grateful for some advice/guidance/comment on the following problem. The imputations are created by a call to `mice()` :               require(mice)     imp <- mice(...)       The details omitted as it a lot of code and a large dataset. Hopefully this won't detract from the question. This is the plot which is generated by the built-in `densityplot` function in the `mice` package: `densityplot( x=imp , data= ~ age)` ![enter image description here](http://i.stack.imgur.com/qOw43.png) The blue line is the observed data, the red lines are the imputed data. This caused me some alarm. Particularly:   1. Since the observed data are fully contained in each imputed dataset, how can some values which appear in the observed data have a zero density in the imputed data ?    2. There are 7019 observations of which only 11 are missing, so I would expect the imputed densities to be nearly identical to the observed.   3. In a general sense, how can the plots look so different ? So I compared it to a plot using `ggplot`:               require(ggplot2)     require(reshape)          fortify.mids <- function(x){      imps <- do.call(rbind, lapply(seq_len(x$m), function(i){            data.frame(complete(x, i), Imputation = i, Imputed = "Imputed")          }))          orig <- cbind(x$data, Imputation = NA, Imputed = "Observed")      rbind(imps, orig)     }          x11()     ggplot(fortify.mids(imp), aes(x = age, colour = Imputed,         group = Imputation)) +     geom_density() +      scale_colour_manual(values = c(Imputed = "#000000", Observed = "#D55E00"))      ![enter image description here](http://i.stack.imgur.com/uqes3.png) And as you see, all the densities overlap (you can't even really notice that there is more than 1 line). Can anyone explain what is going on with the first plot ? Note that the data used to generate these 2 plots is the same.


title: Show that a comparison of the two way drugs can be made by means of the statistic $(b-c)^2/{b+c}$

text: > suppose two drugs are administered to each of _n_ patients suffering from > headaches. The reactions are summarized in the following Table:                                            Drug 1                       Good      Not Good   Total               Good     a          b         a+b     Drug 2  Not Good   c          d         c+d               Total    a+c        b+d         n      > Starting from the marginals (a+b) and (a+c), show that a comparison of the > two way drugs can be made by means of the statistic $\frac{(b-c)^2}{b+c}$, > which is approximately a $\chi_1^2$. Please give me some hints so that I can solve this homework. The problem is a large sample problem.


Do you include a reference category for a series of dummy variables in a probit regression? If so, how would you interpret the reference category? The question underlying my confusion is that I don't understand what purpose the constant term serves in a probit regression. I think the constant is simply the Y assumed as 1. Given that Y is assumed as 1 serving as the constant, what are the marginal effects of the other variables on whether this will occur? If there is a reference category then it seems that the constant is equal to Y assumed as 1 times the probability of the reference category occurring when Y is equal to 1. Please explain this to me and cite your sources so I can read more about it. Thanks!


title: Maximum Entropy and Multinomial Logistic Function

text: I have a newbie question. I tried to find the answer from Google but couldn't get a clear answer. Is the MaxEnt model exactly same as multinomial logistic regression (i.e., softmax regression)? It looks like both try to estimate the parameters of the softmax function. Just wondering, then, what are the difference between them? Do they use a different learning method?


title: Percentile calculation for the combined group

text: I was struggling to find some solution of my problem which is as below: Let say, I have 20 groups and for each group I have say 10the percentile value. Then is it anyway possible to calculate the 10th percentile for the combined group using those 20 individual 10th percentile? If there is not any exact solution then any approximation will also be helpful. Best regards,


title: When does l1 regularisation give a sparse solution?

text: I was maximising a likelihood function, which is convex. I know that the system has a K-sparse solution. I wanted to know the conditions (or some sufficient conditions) on the likelihood function under which l1-regularisation is guaranteed to give me the sparse solution. In other words, I am trying to solve the following problem $$\max_{x \in \mathbb{R}^n} f(x),$$ where $f$ is a concave function. However I know that the true solution is K-sparse. I wanted to know the conditions under which the following optimisation returns the correct solution $$\max_{x \in \mathbb{R}^n} f(x) - \lambda ||x||_1.$$ I found online that if $f=||y-Ax||_2^2$, then it is sufficient for $A$ to satisfy a property called the restricted isometry property. I wanted to know if anything is known in the general case or even special cases where such guarantees were known.


title: Binary Models (Probit and Logit) with a Logarithmic Offset

text: Does anyone have a derivation of how an offset works in binary models like probit and logit? In my problem, the follow-up window can vary in length. Suppose patients get a prophylactic shot as treatment. The shot happens at different times, so if the outcome is a binary indicator of whether _any_ flare-ups happened you need to adjust for the fact that some people have more time to exhibit symptoms. It seems that the probability of a flare-up is proportional to the length of the follow-up period. It's not clear to me _mathematically_ how a binary model with an offset captures this intuition (unlike with the Poisson). The offset is a standard option in both Stata (p.1666) and R, and I can easily see it for a Poisson, but the binary case is a bit opaque. For example, if we have \begin{equation} \frac{E[y \vert x]}{Z}=\exp\\{x'\beta\\}, \end{equation} this is algebraically equivalent to a model where \begin{equation}E[y \vert x]=\exp\\{x'\beta+\log{Z}\\}, \end{equation} which is the standard model with the coefficient on $\log Z$ constrained to $1$. This is called a _logarithmic offset_. I am having trouble figuring out how this works if we replace $\exp\\{\\}$ with $\Phi()$ or $\Lambda()$. Update #1: The logit case was explained below. Update #2: Here's an explanation of what seems to be the main use of offsets for the non- poisson models like probit. The offset can be used to conduct likelihood ratio tests on index functions coefficients. First you estimate the unconstrained model and store the estimates. Say you want to test the hypothesis that $\beta_x=2$. Then you create the variable $z=2 \cdot x$, fit the model dropping $x$ and using $z$ as an non-logarithmic offset. This is the constrained model. The LR tests compares the two, and is an alternative to the usual Wald test.


title: Drawing numbered balls from an urn

text: PROBLEM There is an urn with a set of balls where each ball is labeled with a different integer. The numbers on the balls are known and are not a range of integers. For example the set of balls could be $B_k$ = {1,4,67,3,12}. There are two robots, $A$ and $B$, which draw a set of $n$ balls with replacement. Let $a_i$ represent the set of number values on the balls drawn for robot $A$ where $i$ indexes the $n_a$ draws. $b_i$ and $n_b$ are analogously defined for $B$. The robots have some information about why they draw each ball and how many they will draw which we know nothing about except that the method for decision is similar between $A$ and $B$. Robots $A$ and $B$ draw sets of balls for $N_a$ and $N_b$ rounds where $N_a \approx N_b \approx 10^6$. Let $j$ index the rounds, so $n_{aj}$ balls are drawn each round and the ball numbers are recorded in $a_{ij}$ for $A$. $n_{bj}$ and $b_{ij}$ are defined analogously for $B$. The objective is to choose which robot is expected to be the optimal choice for a new sampling $x_{ij}$ where optimal maximizes the sum over all balls, $\Sigma_j \Sigma_i x_{ij}$. Simply put, which robot would you expect to draw the highest sum of ball numbers in the next million rounds? ADDITIONAL INFORMATION Since we have $a_{ij}$ and $b_{ij}$ we can say a bit about their distributions. The samples drawn on each round $a_i$ and $b_i$ have 0 balls drawn about 95% of the time. The second most common number of balls drawn per round is 1 and then 2 and 3 and so on so that a histogram balls per round, $n_j$, is monotonically decreasing. Some rounds have multiple draws of a highly numbered ball which leads to outliers in the distribution of sum per round, $\Sigma_i a_{ij}$. This means that the sum of ball values per round is very hard to model analytically. There could be more information I am not thinking of using since this is a vast simplification of the real problem so please ask. STANDARD SOLUTION The standard solution given to me is to make the decision based on the mean and the standard error on mean for the sum of the balls in each round. Define: $\mu_a = \Sigma_j \Sigma_i a_{ij} / N_a = mean_j(\Sigma_i a_{ij} )$ and $ \sigma_{\mu_a} = \sqrt{var_j(\Sigma_i a_{ij} )/ N_a}$. The decision and confidence of that decision is calculated by comparing gausians $\mathcal{N}(\mu_a,\,\sigma^2_{\mu_a} ) $ and $\mathcal{N}(\mu_a,\,\sigma^2_{\mu_a} )$. The motivation for this choice over $\mu^*_a = \Sigma_j \Sigma_i a_ij / \Sigma_j n_{aj}$ is that the probability for a second ball to be drawn in a round given that one already has ($\approx 0.2$) is much higher than the probability the first ball is drawn ($\approx 0.05$). This smells of Bayesian motivation but I am unable to do the derivation showing $\mu$ to be the empirical Bayes estimator. QUESTION Is the mean $\mu_a = \Sigma_j \Sigma_i a_ij / N_a$ really the best estimator and why? What about the sensitivity to outliers? Why not the mean over all the balls drawn, $\mu^*_a = \Sigma_j \Sigma_i a_ij / \Sigma_j n_{aj}$, which is less sensitive to outliers? It seems a lot of information is lost by not using the knowledge of the set $B_k$ so is something more Baysian advisable?


title: Graphing versus further testing interactions

text: Following the question here, I have come across different (textbook) recommendations. Some authors suggest that the effect of the moderating variable is best understood in a graphic form. Others suggest that after graphing the effect, further significance testing should be undertaken to ascertain if the pattern (in the graph) is real or not (unsure if 'real' is the right word here). My confusion is that if I have found that the interaction of the moderating variable on the the X-Y relationship is significant, **why should this further testing be required.** In other words, the graph should be a sufficient explanation of the moderating effect. After all, it does show the pattern and the initial test showed that the pattern is significant. Are there any best practice in this regard. I understand, of course, that this depends on the context of a study but it always good to know what the common thinking is on this matter. As always, many, many thanks.....


title: How to deal with identical p-values with the Benjamini-Hochberg method for correcting for multiple testing

text: I'm testing for mutations in DNA and I'm using the Benjamini-Hochberg method to modify the threshold I'm testing my p-values against. The method is basically to rank the p-values and compare them to a new threshold defined by $q(k) = k/n*\alpha$, where $\alpha$ is your original threshold, commonly 0.05, $k$ is the rank of the p-value you compare with, and $n$ is the total amount of p-values (== total amount of tests). Due to the nature of my data, many p-values are identical. The following ordered set of p-values could be an example: $p_1$ = 0.01 $p_2$ = 0.03 $p_3$ = 0.03 $p_4$ = 0.03 $p_5$ = 0.09 Assuming a threshold of $\alpha$ = 0.05, what I have to compare my values to are: $q_1$ = 0.01 $q_2$ = 0.02 $q_3$ = 0.03 $q_4$ = 0.04 $q_5$ = 0.05 And thus I accept the second test and reject the third and fourth even if they have the same original p-value. I tried to search the literature for a solution to this, but nothing came up. Is there an established method for this? And, if not, what is your preferred ad hoc solution?


title: t-stats in OLS regression

text: Very basic question: suppose I run an OLS regression $Y$ ~ $X_1 + X_2$ and obtain t-stats of $t_1$, $t_2$ respectively. If $t_1$ >> $t_2$ (such as an order of magnitude larger), is it safe to conclude that $X_1$ has significantly more predictive power on $Y$, than $X_2$? Thanks.


title: Finding a linear relationship without linear model

text: (If this doesn't make sense it is because I am inexperienced at statistics) Assume I have a broad universe of actions $A$ I could pursue, and for each $a_i \in A$ there is an associated characteristic $x_i \in X$, where $X$ is drawn from a normal distribution centered at 0 with standard deviation 1. Let's also say there is an associated value $v_i$ for each $a_i$ that in this case simply follows the equation: $$v_i = .13x_i+\epsilon,$$ where $\epsilon$ is just a noise term, drawn from a normal distribution with mean $0$ and standard deviation $\sigma$. Further assume that I am not aware of this relationship, but I want to figure it out. One more constraint: I am not allowed to use linear regression. What I've been doing is selecting $n$ $a_i$'s with characteristic $x_1$ and $n$ $a_j$'s with characteristic $x_2$, taking the means of their values $v_i$ and $v_j$: $u_1$ and $u_2$, and inferring the line from the points $(x_1,u_1)$ and $(x_2, u_2)$. The problem is, the line is usually off by a lot, especially when the standard deviation of the noise, $\sigma$, gets larger. To see why this is I constructed a worst case scenario. Basically we are just sampling from the normal distributions for the $v$'s corresponding to $x_1$ and $x_2$. Therefore we can construct confidence intervals around each point. I figured the worst case scenario would be inferring a line from the bottom of one confidence interval to the top of the other. I want to be 95% confident that the actual line will be within these bounds, so the confidence interval around each point needs to have confidence $1-p$ where $p^2/2=.05$, meaning $p=.316$ and we need a 68.4% confidence around each mean. Now using Chebyshev's inequality, we can solve for the number of points we have to draw to get within an acceptable error $\epsilon_a$: $$P(|u_i-\bar{v_i}|\geq\epsilon_a) \leq \frac{\sigma^2}{n\epsilon_a^2} = .316.$$ Solving for n we get: $$n = \frac{\sigma^2}{.316\epsilon^2_a}.$$ Now putting actual parameters in, for the real standard deviation we have $\sigma = .1$, for acceptable error, how about within 1% of the actual $\Delta v$, which is a function of the difference in $x$'s, which we'll say is 1 (regression can easily do this), so we get: $$n = \frac{.1^2}{.361(.01*.13*(x_1-x_2))^2} \approx 16000.$$ The problem is, it takes me anywhere from 5 seconds to 5 minutes to generate each triple of $(a_i, x_i, v_i)$ (there are a lot of things going on in the background), so doing up to 32000 iterations is very tedious and unmanageable. Is there a more efficient way of doing this, besides using a linear model? Also, is there a name for the algorithm that I am using so I can better research it?


title: "Brute force" expected deviance for logistic regression?

text: A commonly used goodness of fit statistic for logistic regression is the deviance. This is also known as the likelihood ratio chi-square statistic. It is defined as: $$D=\sum_{i=1}^{N}d_i^2$$ $$d_i^2=2y_i\log\left(\frac{y_i}{n_ip_i}\right)+2(n_i-y_i)\log\left(\frac{n_i- y_i}{n_i-n_ip_i}\right)$$ Where $N$ is the number of binomial cells, $y_i$ is the response, $n_i$ is the number of units in the ith cell, and $p_i$ is the predicted probability, given by $\log\left(\frac{p_i}{1-p_i}\right)=x_i^T\beta$. I have been taught that we have, approximately: $$E(D)\approx N-\dim(\beta)$$ This is based on using the chi-square approximation to the likelihood ratio. This relationship is then used to test for "over-dispersion" in the model. If the ratio $\frac{D}{N-\dim(\beta)}$ is significantly different from $1$ (using the chi-square distribution) we conclude the model exhibits over-dispersion. Is there any reason why we don't simply use a brute force approach and directly average out $y_i$ from the deviance residual? That is, for each deviance residual calculate the condition mean exactly via: $$E(d_i^2|p_i)=\sum_{y_i=0}^{n_i}{n_i\choose y_i}p_i^{y_i}(1-p_i)^{n_i- y_i}\left[2y_i\log\left(\frac{y_i}{n_ip_i}\right)+2(n_i- y_i)\log\left(\frac{n_i-y_i}{n_i-n_ip_i}\right)\right]$$ This is not a intensive calculation to do (or to code), and can easily be used to calculate the conditional variance of $d_i^2$ as well, given that most binomial cells are small. Further we can replicate the model dispersion scenario by calculating $\frac{d_i^2}{E(d_i^2|p_i)}$, or calculate T-like quantities $\frac{d_i^2-E(d_i^2|p_i)}{\sqrt{V(d_i^2|p_i)}}$ and search for outliers. I feel like I must be missing something. Any pointers?? has anyone seen this done elsewhere?


title: How best to represent a web app's usage?

text: I've been gathering data for the usage of a web app. For each day at the same hour and over several months I have:   * (A) total number of members    * number of online members and its % to (A)   * number of members who connected during the day and its % to (A)   * number of new members    * number of forum members    * number of new forum members    * number of chat participants  Plus the periods where it was on sale (there's one time fee to access it). What can I show with those data and how best to represent it graphically ? (I'll be using Apple Numbers.)


title: Comparings BIC of lme models with and without a correlation function

text: I have created an `lme` model using the same predictors both with and without a specified correlation structure based upon distance between the points (lat/long). I know that I am able to compare the AIC/BIC between different correlation structures (corSpher vs. corExp vs. corGaus), but am unsure if I can directly compare the AIC/BIC values of a model with no correlation (`lme(y~x1+x2, random=~x1|A, data = K)`), and one with a correlation (`lme(y~x1+x2, random=~x1|A), corr = corGaus(form=~long+lat), data = K)`) structure? If I am not able to compare the AIC/BIC, are there any recommend methods in comparing the fit/best model? Any advice would be greatly appreciated. Thanks


title: Good books on text mining?

text: Hi I wanted to know if there are some good books on text mining and classification with some case studies?. If not some papers/journals accessible to public would do. If they illustrate their examples with R even better. I am not looking for step by step manual but something which illustrates the pros and cons of various text mining approaches to various classes of problems.


title: Normalizing constant & rejection sampling

text: I think I understand what a normalizing constant is. Say for example you have a pdf $f(x)$ with support $0 \le x \le 5$. If you wanted to truncate the pdf and only look at $0 \le x \le 3$ you would need to use a normalizing constant $c$ such that $\int_0^3 cf(x)\,dx = 1$. The problem I'm dealing with states: > In class we gave a rejection method for the case when the target > distribution π(x) is only known up to a normalizing constant. Show that the > accepted samples from that algorithm indeed follow the target distribution > π(x) What does it mean when it says we only know the target distribution up to a normalizing constant?


title: Why isn't this code for prop.test giving me the right answer for test statistics?

text: I have used the following line of code but it doesn't give me the right answer for (`chi-squared`) `test statistic` while it gives the right answer for `p-value` regarding the following question: "Consider testing H_0: p = 0.5 against H_1: p ≠ 0.5, given that there were 38 successes in a sample of size 56. Find the value of the (chi-squared) test statistic. (Do not use the continuity correction.)"               > prop.test(38,n=56,p=1/2)         1-sample proportions test with continuity correction          data:  38 out of 56, null probability 1/2     X-squared = 6.4464, df = 1, p-value = 0.01112     alternative hypothesis: true p is not equal to 0.5     95 percent confidence interval:      0.5390880 0.7935168     sample estimates:             p      0.6785714       Can you give some explanation in your answers? I am not coming from a statistics background.


title: Missing stems in stem-and-leaf plot

text: ![Stem-leaf plot](http://i.stack.imgur.com/ObRA8.png) Why is there no 17| xxx in the stem leaf plot? Thanks.


title: Computing Cronbach’s Alpha for Test with Random Items

text: I need to find out Cronbach's Alpha for a test that pulls 10 questions randomly from a pool of 40 questions for every subject that takes the test. As a result, almost every test has a different set of questions. I thought that I could compute CA for the 40 questions and leave unanswered items as missing values, but I'm afraid this might affect the variance of the items. The sample is not big enough to find a subset of subjects that has answered exactly the same questions (n=50).


title: Zero-inflated negative binomial

text: I am trying to understand zero-inflated negative binomial regression. My impression is that if a zero-inflated negative binomial model does not contain any logit part, the model is identical to the one can obtain with just ordinary negative binomial regression. Is this correct? Thank you PS: the logit part I was talking about - well - zero-inflated model assumes that the 0s within the dataset are generated based on two different process: one is negative binomial and the other is, if I remember it correctly, poisson. By "no logit part" I meant what if we take out the effect of the poisson distribution from the zero-inflated model? would it be same as ordinary negative binomial regression?


title: Why don't dummy variables have the continuous adjacent category problem in cluster analysis?

text: I know that if we use categorical variables in cluster analysis we would assume that the scale is continuous and we don't have this concept of distance between two adjacent categories. But what is the difference when you use dummy variables? The zeros and ones will be used anyway to calculate the distances in the cluster analysis. In a nutshell, why don't 0's and 1's have this same issue? Any references about it? Thanks


title: Tips for rejected Ljung-Box stats? Can Minitab model ARIMAX?

text: Two questions:   1. What are some basic tips on going back to the drawing-board after a seasonal ARIMA fit shows significant coefficients (1,0,0)x(0,1,0) but fails in its Ljung-Box stats?    2. Can minitab be used to model an ARIMA with regression?


title: Interpreting regression results for aggregate data: does it make sense to refer to the "average of average"?

text: I don't have access to individual data, so I have to rely on aggregate data for a large number of cities in a country. If my outcome variable is the log of city average earnings and the predictors are dummy variables for n-1 national regions, does it make sense to interpret coefficients in this way: "Average earnings are on average 10.3% higher for Region_B compared to Region_A" ? In other words, does it makes sense to speak of "average advantage in average earnings" , or is that like saying "average of an average" ?


title: Clustering with 3 attributes

text: Please bear with me because I am very new to data mining. I have a database of 3 attributes: latitude, longitude and temperature. I want to find clusters for the temperature data and I also want to include the effect of latitude and longitude in that so that temperature is not the only determining factor for the clusters. If I, let's say, build K-means clusters using all these attributes (in WEKA for example), what do the resulting clusters tell me? Can I get any interpretation of how latitude and longitude information is related to clusters of temperature? What is the correct way to go here?


title: Learning hidden Markov model where transition/emission/initial probabilities aren't independent

text: I'm working on a problem that I've cast as an HMM, except that unlike the "traditional" case where the transition probabilities $a(i,j) = p(s_i = j \,|\, s_{i-1}=i)$, emission probabilities $b(j,o) = p(o_i | s_i = j)$ and initial probabilities $\pi(j) = p(s_0 = j)$ are all just independent numbers, in my case they are (relatively complicated) functions of a set of common parameters. I'm reasonably sure that this means I can't rely on Baum-Welch to find the optimal parameters (because Baum-Welch assumes that the three sets of probabilities can be optimized independently), and I've tried my hand at directly optimizing the likelihood by just computing the likelihood using the forward algorithm. However, this has proved slow (my sequence is extremely long) and somewhat unreliable (the likelihood surface is a bit weird looking). A friend suggested that both of those problems could be solved if I used an EM algorithm (using forward-backward probabilities in the E-step and numerically optimizing in the M-step). However, I'm having a hell of a time figuring out exactly what the Q function should be. I believe it is $$ Q(\theta\,|\, \theta^{t}) = \sum_j\log(\pi(j))p(s_0 = j \,|\, x, \theta^{t}) \\\\+ \sum_t\left[\sum_{i,j} \log(a(i,j))p(s_t = j, s_{t-1} = i \,|\, x, \theta^{t}) \right]\\\\+\sum_t\left[\sum_j\log(b(j,o))p(s_t = j \,|\, x, \theta^{t}) \right]. $$ However, it's not clear to me that this will make life any easier, as it still looks like it's $O(nm^2)$, where $n$ is the length of the sequence and $m$ is the number of hidden states, time to evaluate. And that's in addition to forward-backward run to calculate all the posterior probabilities necessary during the E step. So I feel like I must be doing something wrong or else missing the point of using the EM algorithm rather than directly optimizing the likelihood in this problem.


title: How to understand a convolutional deep belief network for audio classification?

text: In "Convolutional deep belief networks for scalable unsupervised learning of hierarchical representations" by Lee et. al.(PDF) Convolutional DBN's are proposed. Also the method is evaluated for image classification. This sounds logical, as there are natural local image features, like small corners and edges etc. In "Unsupervised feature learning for audio classification using convolutional deep belief networks" by Lee et. al. this method is applied for audio in different types of classifications. Speaker identification, gender indentification, phone classification and also some music genre / artist classification. How can the convolutional part of this network be interpreted for audio, like it can be explained for images as edges?


title: Computing confidence intervals for population variance from a sample in R

text: Is there a package available for R (on CRAN, github, r-forge, etc.) that computes CIs for the population variance, given a sample of data, 95% CI parameter, etc.? The `Rmisc` package computes the CIs around the mean, but I'd also like the CIs around the variance. I found a few papers and an online source that contains the R code necessary to do this, and I've already used that code to roll my own implementation. However, given how large and rich the R package ecosystem is, I'm heavily biased towards not rolling my own solution, particularly when it comes to computing statistics measures.


title: Grouping subject data for within-subject analysis in SPSS

text: I have over 30 subjects who each did the same task about 50 times. Each of these data point can be in condition a1 or a2 of variable A or condition b1, b2, b3 or b4. So for each subject, I have datapoints for a1b1, a2b2, etc. I have three dependent variables (v1, v2, v3) for each data point. I understand that this is a within-subject analysis, done through a repeated measure test. But the way I understand it, if I want to compare a1 to a2, I have to make a column for a1 with each subject's average for a1, the same for a2, and compare those two columns. And I have to do it for each dependent variable, and then again for b1-b4. And if I want to study the interaction, I have to do the same for a1b1, a2b1, a1b2, and so on, again for each dependent variable. This seems quite impractical. As of now, my data is set with a column for subject's ID, a column for A condition, a column for B condition, a column for each dependent variable v1 v2 and v3, and a row for each data point. Is there a way to do the repeated measure analysis as follows?               "(average of v1 when A=a1 for each subject) vs (average of v2 when A=a2 for each subjet)"


title: What graphical techniques are used in Structural Equation Modeling?

text: I'm curious if there are graphical techniques particular, or more applicable, to structural equation modeling. I guess this could fall into categories for exploratory tools for covariance analysis or graphical diagnostics for SEM model evaluation. (I'm not really thinking of path/graph diagrams here.)


title: When do you use AIC vs. BIC

text: How do you know when to use AIC or BIC for determining model fit? Is it just a judgment call? Is there an intuitive explanation as to which heuristic is better than the other?


title: Interpretation of AIC, BIC and KIC

text: > **Possible Duplicate:**   >  Is there any reason to prefer the AIC or BIC over the other? Can anyone please interpret each term in AIC, BIC and KIC. And the difference between the three. Thanks in advance,


title: Post-hocs for within subjects tests?

text: What is the preferred method for for conducting post-hocs for within subjects tests? I've seen published work where Tukey's HSD is employed but a review of Keppel and Maxwell & Delaney suggests that the likely violation of sphericity in these designs makes the error term incorrect and this approach problematic. Maxwell & Delaney provide an approach to the problem in their book, but I've never seen it done that way in any stats package. Is the approach they offer appropriate? Would a Bonferroni or Sidak correction on multiple paired sample t-tests be reasonable? An acceptable answer will provide general R code which can conduct post-hocs on simple, multiple-way, and mixed designs as produced by the `ezANOVA` function in the `ez` package, and appropriate citations that are likely to pass muster with reviewers.


title: What does "equal a priori class probabilities" mean?

text: I am trying to solve a problem about my homework. The problem says that > Assume a two-class problem with equal a priori class probabilities Does it mean, mean vectors and covariance matrices should be equal? The question of the homework is: Assume a two-class problem with equal a priori class probabilities and Gaussian class-conditional densities as follows: $$p(x\mid w_1) = {\cal N}\left(\begin{bmatrix} 0 \\\ 0 \end{bmatrix},\begin{bmatrix} a & c \\\ c & b \end{bmatrix}\right)\quad\text{and}\quad p(x\mid w_2) = {\cal N}\left(\begin{bmatrix} d \\\ e \end{bmatrix},\begin{bmatrix} 1 & 0 \\\ 0 & 1 \end{bmatrix}\right)$$ where $ab-c^2=1$.


title: How to analyse data on metric space?

text: I just ran an experiment and I'm not sure how best to analyse the data. My data are distance values between objects in a metric space bounded by [0,1]. I have drawn up a probability density estimate as follows:![enter image description here](http://i.stack.imgur.com/pqIhZ.png) This looks to me like a kind of hybrid between an exponential and normal distribution, what tests can I run to make better sense of this? * * * _More background:_ This metric space is a finite corpus of Documents, the distance between which indicates their similarity: if the distance is 0 they are identical, if 1 they have no commonality. This sample represents all the distances between 1,000 randomly selected documents.


title: Can I use log-likelihood distance on data of only continuous variables?

text: I have to run a SPSS two-step cluster analysis. All my 4 variables are continuous scalar standardized parameters (with normal distribution). The dataset includes 10,000 cases. SPSS suggest to use euclidean distance with such a dataset, but the resuls are not significant (2 clusters: 99% and 1%), while using the log-likelihood distance option the clusters seem much more meaningful (both if I specify a fixed number of clusters and if I do not). **Question:** Which may be the reason of such a meaningless results with euclidean distance? maybe noise handling? And is it incorrect to use the log-likelihood distance even if my variables are all continuous?


title: Controlled experiment with binary outcome variable

text: I have data from an experiment where students were tested on multiple words for the correct pronunciation of a specific linguistic phenomena. The experiment was done with a control and experimental group, and all students were tested prior to a specific instruction and at the end of the experiments (roughly 8 weeks later).               > head(act1)       studentid groupid itemid test type result     1        1B       B      1    0    1      0     2        5B       B      1    0    1      1     3        6B       B      1    0    1      1     4        8B       B      1    0    1      0     5       11B       B      1    0    1      1     6       15B       B      1    0    1      1          > levels(act1$groupid)     [1] "B" "D"      D is the experimental group while B is the control group. First, I want to test the hypothesis that the experimental group has had a significant improvement in the production of the studied linguistic phenomena. I ran McNemar's test on each group individually:               > act1wide <- reshape(act1, idvar=c("studentid","groupid","itemid","type"), timevar=c("test"), v.names=c("result"), direction="wide")     > act1wideb <- subset(act1wide, groupid=='B')     > act1wided <- subset(act1wide, groupid=='D')     > mcnemar.test(act1wideb$result.0, act1wideb$result.1)              McNemar's Chi-squared test with continuity correction          data:  act1wideb$result.0 and act1wideb$result.1      McNemar's chi-squared = 0.0556, df = 1, p-value = 0.8137          > mcnemar.test(act1wided$result.0, act1wided$result.1)              McNemar's Chi-squared test with continuity correction          data:  act1wided$result.0 and act1wided$result.1      McNemar's chi-squared = 9.0312, df = 1, p-value = 0.002654      This seems to show that group 'B' did not improve while group 'D' did. **Questions:**   * Is this a valid test?   * How can I have a more meaningful test that would compare group data at the same time?   * The factor _type_ has three levels. How can I break down the analysis to get specific information by type? I am a bit of a stats noob so I greatly appreciate detailed explanations on how to make sense of my data.


title: How to test for increases in categorical variable values measured over time?

text: attached is a (German...) bar chart which shows the five possible values of the one variable "Culture": Nix, BV, Candida, Mix, Unspecific.   The variable is measured at three different times, 1 (blue), 2 (green), 3 (something). I am looking for some kind of test to determine that the number of measured individual values increase, or decrease. For example, I'd like to show that over time (i.e. from time=1 to time=3), the value "Nix" is measured more often, whereas "Candida" and possibly "Mix" go down. I am a bit confused because we have three time points and five possible values, this makes it a bit unclear for me. ![Bar chart](http://i.stack.imgur.com/pjeBq.jpg)


title: Metrics of correlation for 3+ variables

text: I have a basic question, which hopefully you all can resolve for me. What is the best way to determine correlations between 3 or more variables? I have a dataset in which 5 continuous variables each correlate (Pearson, R ~= 0.7) with another continuous variable (which I'll call `Z`). I know I can use partial correlation between groups of 2 variables and `Z`, but I don't think this is exactly what I'm looking for, if I understand it correctly.


title: Sample size calculation for t-test

text: I need help calculateing the sample size for a study. I am using a group of students in a particular gradelevel. My study will be pretest-posttest type. I only have about 45 students in that grade, but my advisor said I have to perfrom a sample size. I will like to know the effects of Polya's four-step problem solving on the mathematical understanding of word problems. Can someone please help me. Thank you in advance.


title: Markov Chain process - Missing values

text: I have a vector of sequences with presence (`1`) and absence (`0`), from were I have calculated the first order Markov Process. This is how the data looks like:               dataset=c(NA,NA,0,0,0,0,1,0,1,1,0,1,1,NA,NA,NA,NA,NA,0,               1,1,1,1,1,1,0,NA,1,1,1,1,1,1,0,0,NA,NA,0,1,1,               NA,NA,0,NA,0,0,0,1,1,1,0,1,1,1,1,0,1,NA,0,1)      For thsi vector I have calculated:   * Probability of Presence (`P1`) and Absence (`P0`)   * Probability of having Presence followed by Presence (`P11`), Presence followed by Absence (`P10`), Absence followed by Presence (`P01`) and Absence followed by Absence (`P00`) The transition probabilities were obtained by using a loop the checks for sequence of 2 values: to calculate `P_00` I am using                  P_00 : dataset[j]==0 & dataset[j+1]==0 , etc.      These are the results:               P_0=0.3913043     P_1=0.6086957     P_0+P_1=1          P_00=0.1538462     P_01=0.2307692     P_11=0.4615385     P_10=0.1538462     P_00+P_01+P_11+P_10=1      The idea is to populate the NA with presence/absence sequences according to the Probability values obtained for the sequence. The problem is that for now I am not being able to find the best and more adequate process for this problem and I am a newbie with this type of problems. Even if there is already a package that can do this for me, I would prefer to understand how can I do it.


title: Controlling unobserved determinants in time series regression

text: I'm trying to estimate what effect the tariff rate in a specific country has on evasion of tariffs (or rather the underreporting of the value of goods in customs). My dependent variable is a proxy for evasion and the explanatory one is tariffs. I have data spanning continuously over about 15 years. Now I want to control for other possible unobservable determinants of evasion (such as an increase in customs enforcement or some technological improvement) as robustness and am not sure what would be appropriate (pretty new to econometrics). Would year dummies do the trick (and is this equivalent with "time/year fixed effects"?)? And finally, what would be the difference with estimating a model in first difference? (Can you do this for multiple years as in my dataset?) Any input would be much appreciated! Thanks Oscar


title: Selection of priors for a BYM spatial regression model

text: I am using a BYM model in WinBugs to describe the distribution of a non- infectious disease. The model at present is a standard enough BYM model without much modification, (a Poisson-gamma hierarchical model with terms for structured and unstructured residuals) as follows:               model     {       for(i in 1:N)       {       observed[i] ~ dpois(mu[i])         log(mu[i]) <- log(expected[i]) + alpha + u[i] + v[i]            theta[i] <- exp(alpha + u[i] + v[i])       u[i] ~ dnorm(0, tau.v)            }            v[1:N] ~ car.normal(adj[], weights[], num[], tau.v)            alpha ~ dflat()       tau.u ~ dgamma(0.5,0.0005)       tau.v ~ dgamma(0.5,0.0005)            mean <- exp(alpha)          }      My question is about the selection of priors for `tau.u` and `tau.v`. They are both set at (0.5,0.0005) in all the examples I can find. How are these values determined? Thanks.


title: Acceptance Rate of Statistical Journals

text: I am interested to know if there is any list that indicates the **Acceptance Rate** (i.e. Total # of submitted manuscripts / total # of accepted manuscripts) of the Statistical Journals? By Statistical Journals I mean the journals like the ones listed here. Note: I am well aware of the IF (impact factor) of a given journal but I am interested in the Acceptance Rate.


title: Proof that regression residual error is an unbiased estimate of error variance

text: Consider the least squares problem $Y=X\beta +\epsilon$ while $\epsilon$ is zero mean Gaussian with $E(\epsilon) = 0$ and variance $\sigma^2$. I need to prove that $\frac{V(\hat{\beta})}{N-(n+m)}$ is an unbiased estimate of $\sigma^2$ with $V(\beta) = ||Y-X\beta||$ . I wasn't able to find the answer online. I just got confused by a thousand different ways to write things down. EDIT:   $Y = \begin{pmatrix} y(0)\\\ \vdots \\\ y(N-1)\end{pmatrix} \quad$ $X = \begin{pmatrix} x^T(0)\\\ \vdots \\\ x^T(N-1)\end{pmatrix}\quad $ $\beta = \begin{pmatrix} a_1\\\ \vdots \\\ a_n\\\ b_1 \\\\\vdots \\\ b_m \end{pmatrix}$


title: z-test vs t-test under assumptions

text: A traditional t-test to compare two means assumes independent random variables, $X_1,...,X_n \sim n(\mu_1,\sigma_1 ^2)$ and $Y_1, ...,Y_m \sim n(\mu_2,\sigma_2 ^2)$. One usually assumes normal distributions with $\sigma_1$ = $\sigma_2$. Then the statistics $$\frac{\bar{X}-\bar{Y}}{S\sqrt{1/n+1/m}}$$ is used and it's observed value is compared to $t_{n+m-2,1-\alpha/2}$ where $S^2$ is the pooled estimate of variance, valid under both the null and alternative, given the assumption of equal variances. If the variances are known one can use $$\frac{\bar{X}-\bar{Y}}{\sqrt{\sigma_1^2/n+\sigma_2^2/m}}$$ and compare with a $z_{1-\alpha/2}$ test. If we assume $n=m$ and for example $\sigma_1^2=1$ and $\sigma_2^2=2$ and investigate using simulation which test is better. How would we do that. It seems wise to check level and powers but how would one do that?


title: Why does R use standard deviation instead of variance in normal distribution?

text: The title says it all. When working with the normal distribution, we use variance, e.g. $X \sim N(\mu, \sigma^2)$ but working in R, we have to specify the standard deviation:               qnorm(0.5, mean = mu, sd = sqrt(sigma^2))      Why?


title: Quick and easy way to remove outliers

text: I have a set of data of locations and associated rent prices. Now there seem to be several outliers which I would like to get rid of so that a plot of my original data gains more meaning. In the "world" of statistics, would it be acceptable if I did this by eliminating any prices that who deviate from the mean by more than twice the standard deviation? The aim of what I'm working on is to test out several machine learning techniques. I don't have to be extremely accurate but I wouldn't like to do something that is totally unacceptable.


title: Hierarchical decomposition of an imbalanced multiclass classification problem

text: I have a heavily imbalanced multiclass text classification problem: one class is very probable a priori (`P`), while the remaining four ones are about equally improbable (`I1` to `I4`). I have performed many experiments using a single multiclass Random Forest, trying to model all classes at once (with the goal of optimizing the negative loglikelihood). The heavily imbalanced structure of that particular problem made me think however that a better strategy would be to decompose it in a hiercharchy, by conflating the improbable classes into one (`I`): a first-level RF could be used to model the `P` vs `I` binary problem, while a second-level RF could focus on the specialized four-class problem. It's then easy to combine the outputs of the two models into one. My intuition was that this decomposition should work well, because each model has a simpler task than the overall one, with more training information (i.e. being less diluted in a wider and unevener set of classes). But even though I've been very careful to make my two implementations very similar and comparable, I found that the single multiclass model clearly outperforms the hierachical one, contrary to my intuition. I'd like to know if there's a way to understand this result.


title: What is the difference between function createFolds and createDataPartition in caret

text: I think the two functions are all to split the data, but I really can't get the difference between both of them, even I has read the help manual of them.


title: Finding underlying distribution of proportion data

text: I have a large dataset of count data --- count of positive cases and total number of cases distributed by distance (in meters). The empirical distribution of proportions and beginning of data are attached (there are 444260 rows, red line is the proportion of cumulative sum) I want to find an underlying, theoretical distribution that fits best the data. I would like to have a function that estimate probability of the positive case given distance in meters. However, due to sparsity of data (small total number of cases on a specific distance), some proportions on far distance are high because of randomness. What would be the best way to transform dataset, to smooth the distribution? I’ve tried various binning methods, Fayyad Multi-Interval Discretization, log scale, but all those methods seems to “rough” to capture the steep decline on beginning. What would be the best candidates for theoretical distribution --- Poisson, negative binomial, Pareto, …? Thanks, David * * *                   D   K   CNT      1  0   167 1622      2  1   0   12      3  2   0   1      4  3   0   1      5  4   0   1      6  5   0   3      7  6   1   7      8  7   0   8      9  8   2   14      10 9   1   16      11 10  0   10      12 11  0   17      13 12  1   13      14 13  1   31      15 14  0   36      16 15  0   57      17 16  1   77      18 17  3   60      19 18  0   31      20 19  0   56      ![Distance-Proportion plot](http://i.stack.imgur.com/KWevq.png)


title: Determining the reliability of weather forecast

text: basically, I would like to answer the following question: > How well do weather forecasts predict the future weather? For example, I would like to know how likely it is to be raining when forecasts predict 70% raining risk 3 days ahead of time. I know there are plenty of weather data sets, but are there any data sets which contain information about _forecasts_ specifically? Has anything like this been done before? Thank you


title: How can I improve my cross-correlation methodology?

text: I will break my question appropriately to make it easier to explain. This is actually a short question but looks longer due to my plots so I apologize for that. Any input would be greatly appreciated. ### Setup I have a few locations that I am monitoring with sensors. The dataset that I have is a set of events that the sensors measured at different locations. The measurements themselves look like the following: ![enter image description here](http://i.stack.imgur.com/igor8.png) ### First take I wanted to figure out if any of the events recorded by these sensors were correlated so I ended up computing the cross-correlation of these series (plotted below). The plot shows me that there seems to be some correlation between the sequences `V3` and `V10` but nothing to show a strong correlation. So I ended up revising my approach: I used a window to observe the events. ![enter image description here](http://i.stack.imgur.com/QQA6Q.png) I took a time window of 40 seconds and plotted the time series again along with their cross-correlation with each other. This time, it gave me something encouraging. I can see a perfect correlation between 3 and 10 and 3 and 6. ![enter image description here](http://i.stack.imgur.com/7XlMk.png) and its corresponding cross-correlation plots. ![enter image description here](http://i.stack.imgur.com/BzvZ3.png) ### Question It looks like I should not be computing the cross-correlating the entire timeseries but rather only across a small window. This makes sense because I am not really looking to see if the entire time series is the same anyways. What I am really looking for is if there is a correlated event. If we assume that events beyond a certain time period need not necessarily be correlated then we can perhaps compute the cross-correlation for a sliding window that would be more accurate and notify only if it is greater than 0.8. Well, at this point, I feel like I am reinventing the wheel and am hoping there are standard techniques to do this. Most of my data points would be zero with a sudden burst now and then. In that case, how can I make I conclusions (that there is indeed a correlation) stronger? I head about FFTs and DFTs. Would they be of any use here? Any comments?


title: Is there a way to represent datasets using Predictive Model Markup Language?

text: This seems like a dumb question, but does PMML have a way to represent a data set? And if not, why not? There's detailed support for defining the name, type, and legal values of each feature that will appear in the data, as well as which of those features to use in the trained model, and what transformations to apply to the values of those features before using them. However, I haven't been able to find an example, or any discussion, of representing raw data (either training sets or test data) in PMML. Or for that matter representing data in XML at all. _PMML in Action_ (pp. 7-8) explicitly says "...raw input data is usually formatted as a flat file in which columns represent data fields and rows records or transactions." The sample datasets at the DMG website are all in .csv format. This omission seems particularly strange, since PMML has a model format for support vector machines, which are (almost) just weighted sets of examples. So it would have been only the tiniest bit of effort to define format for training sets and test sets.


title: Text mining software (beyond R)

text: I am looking for software for teaching text mining to business students. Unfortunately, I cannot consider programming languages like R. Does anyone have experience with other text mining software? Can you recommend it? (or not!) Ideally, software that is easy to learn for non-technical people, not overly expensive (not necessarily free), which can be used to convey the concepts of text mining clearly (rather than getting bogged down in software details). Thanks!


title: Is ANOVA good for non-normally distributed series?

text: I just read Finding Seasonal Spreads, by Paul Teetor. As far as I understand how ANOVA works, it should be accurate for NORMAL distributed series. Financial series are _not_ normally distributed, or maybe am I wrong? Could someone explain why he used ANOVA in this case?


title: Rejection Sampling

text: Suppose we are using rejection sampling and we want to sample from a distribution, say $p$. In order to calculate the acceptance probability we use the ratio: $$P(u < \frac{p(x)}{Mq(x)})$$ Therefore we still use the distribution of $p$ for the randomly generated values $x$. So why do we use rejection sampling exactly? Why don't we sample directly from $p$? In other words, what makes a distribution difficult to sample from? Thank you.


title: Predictive power (or $R^2$) adjusted for certain variables

text: I will frame this question for Ordinary Least Square (OLS) regression, but my question is for both OLS and Logistic. Let's say we data over 10000 different individuals. For each person we have three variables - $x_1$, $x_2$, $x_3$ and a target variable $y$ which is how much they have spent on the product. Let's say that the data spans over two products, $p=1$ and $2$. **Objective** Create a prediction model for $y$ (spend) using $x_1$, $x_2$, $x_3$, which should be product agnostic. I will use the model only to determine top spenders. Hence I will be happy as long as the model rank orders (i.e. it gives higher score to person A than B if person A indeed is likely to spend more), and am not concerned about the point estimate's accuracy much. **Approach 1** Since there are two products and one product may have intrinsic higher spend (eg. offers more reward points), I will introduce a dummy variable, $d=1_{p=1}$. I will build a model $y=a_0+a_1x_1+a_2x_2+a_3x_3+a_4d+\epsilon$ I will then ignore $a_4$, i.e. set it to zero while scoring since I want product agnostic rank ordering. **Approach 2** Don't use the product dummy. I will build a model $y=a_0+a_1x_1+a_2x_2+a_3x_3+\epsilon$ Which approach is better? I feel Approach 1 should do a better job, since it adjusts for intrinsic difference in spend across products, which I want to discount in the analysis. Metrics like $R^2$ will definitely be more in raw form for Approach 1, since there is an additional variable. However since I cannot use the dummy while scoring, I would like to discount its effect on measuring predictive power. How can we characterize the power of the model discounting the effect of dummy in Approach 1?


title: Shrinkage estimation for regression with ARMA errors

text: I was wondering if someone knows a R-package or function library for the topic of shrinkage for regression with ARMA errors. Please let me know if you came across something related. Thank you! Regards,   Patricia Tencaliec


title: Automating model selection criteria production

text: I want to perform model comparison according to several criteria using R. My dataframe's name is `df`               head(df)             Y  X1  X2   X3   X4     1 18 307 130 3504 12.0     2 15 350 165 3693 11.5     3 18 318 150 3436 11.0     4 16 304 150 3433 12.0     5 17 302 140 3449 10.5     6 15 429 198 4341 10.0      I want to compare all possible combinations with Y as the dependent variable. I'll give you what I've done so far to get the R.squared for all possible models. I would like your comments for better coding, especially on getting all those model formulas               # First get all possible combinations of the 4 inependent variables     comb_list <- lapply(1:4,function(i) combn(4,i))          # now create the formulas     comb_list_forms <- lapply(unlist(       sapply(1:length(comb_list), function(i)       sapply(1:dim(comb_list[[i]])[2], function(x) {nam <- names(df[1+comb_list[[i]][,x]])       formul <- "Y~"      sapply(1:length(nam), function(y) formul <<- paste(formul,nam[y],sep="+"))      formul <- sub("~+", "~", formul,fixed = T)      }))),as.formula)          # finally get r.squared     attach(df)     models.r.sq <- sapply(comb_list_forms, function(i) summary(lm(i))$r.squared)      Thank you


title: Proper way to analyse purchase data

text: I have a dataset containing information about various factors influencing the amount purchased of different products. factors include the number of people who bought it, the degree to which the product was on sale (90 or 80%), the number of times the product was on sale, etc. I would like to see how much such factors influence the popularity of each product (number of people buying it) I was wondering if a simple regression or anova would do the job or whether there are some other specific statistics more suitable for such problems. further problems: I have data from different companies offering different products. I want to make inferences about the relationship between the type of product and the amount purchased by pooling all the data available from the different companies. However, my problem is that the different companies differ in size and therefore the amount they have sold. Can this be solved somehow?


title: Is there other types of mixture distribution besides the normal mixture

text: There are quite a lot of study on the normal mixture distributions, say, $X=Y*Z$,where $Z$ is a normal r.v. and Y is a r.v. follows other distributions and $Y$ and $Z$ are independent. Some well-known distributions are Normal inverse Gaussian, generalized Hyperbolic etc. I would like to know that in probability, is there some other mixture distributions in such form (product of two independent r.v.s) that are well- studied? For example, let $Z$ is possion or exponential and then is there some well-studied probablity distribution can be generated from this form. If possible, could anyone so kind give me some textbook or monograph on this topic? Thanks so much!!


title: How to describe or visualize a multiple linear regression model

text: I'm trying to fit a multiple linear regression model to my data with couple of input parameters, say 3. \begin{align} F(x) &= Ax_1 + Bx_2 + Cx_3 + d \tag{i} \\\ &\text{or} \\\ F(x) &= (A\ B\ C)^T (x_1\ x_2\ x_3) + d \tag{ii} \end{align} How do I explain and visualize this model? I could think of the following options:   1. Mention the regression equation as described in $(i)$ (coefficients, constant) along with standard deviation and then a residual error plot to show the accuracy of this model.    2. Pairwise plots of independent and dependent variables, like this: ![enter image description here](http://i.stack.imgur.com/hHbN7.png)   3. Once the coefficients are known, can the data points used to obtain equation $(i)$ be condensed to their real values. That is, the training data have new values, in the form $x$ instead of $x_1$, $x_2$, $x_3$, $\ldots$ where each of independent variable is multiplied by its respective coefficient. Then this simplified version can be visually shown as a simple regression as this:  ![enter image description here](http://i.stack.imgur.com/r58G1.png) I'm confused on this in spite of going through appropriate material on this topic. Can someone please explain to me how to "explain" a multiple linear regression model and how to visually show it.


title: Fitting Mixture Distribution of Poisson and Negative Binomial (Delaporte Distribution)

text: I am trying to fit a model to some simulated data. The idea is to use ML- Estimation. However, I am totally lost. I have a dependent variable which is a sum of two (unknown) variables. The first part is Poisson-Distributed and the second part is Negative Binomial Distributed. In the sum this leads to a Delaporte model. I already implemented parts of the model based on a log likelihood-function. However, I am missing the weighting matrix.               set.seed(1234)     x1 <- runif(5000, min = 0, max = 5)     x2 <- runif(5000, min = 0, max = 5)          # Model 1     ## beta0=0.5, beta1=0.3, beta2=0.4, eta=0.5, delta=2          alpha1  <- rgamma(5000, shape = 0.5, rate = 0.5)     lambda1 <- alpha1 * exp(0.5 + 0.3 * x1 + 0.4 * x2 )     Q1      <- matrix(1:5000, ncol=1)     for (i in 1:5000){       Q1[i] <- rpois(1, lambda=lambda1)     }     Q1 <- apply(Q1, 1,as.numeric) # converts matrix object to numeric vector     u1 <- rpois(5000, lambda=0.5)     y1 <- Q1 + u1          # ML Estimation     n     <- length(x1)     nlogL <- function(par) {       n     <- length(x1)       beta  <- par[1:3]       eta   <- par[4]       delta <- par[5]       -(-n*eta*log(eta)*n*mean(y1)+sum(delta*log(delta/(delta+lambda)))+sum(log(???)))     }     par0 <- as.vector(c(0.2, 0.1, 0.2, 0.2, 1.3))     out<-nlm(nlogL,p=par0, hessian = TRUE)     out      Where $???$ is needed! $$ ??? = \sum_{q_i = 1}^{y_i} w(q_i) $$ \begin{align} w &= w(q_i) \\\ &= w(q_i| \beta, \eta, \delta) \\\ &= \frac{\Gamma(q_i + \delta)}{\Gamma(\delta)\Gamma(Q_i+1} \bigg[\frac{\lambda_i}{\eta(\lambda_i+\delta)}\bigg]^q_i \frac{1}{\Gamma(y_1-q_i+1)} \end{align}


title: How to divide ordinal set into bins?

text: What is the best way to divide Cross Validated/Stack Overflow/Server fault reputations into separate bins, considering reputation systems tend to be highly skewed? For example: Suppose I want a set of 64 bins - 0-250,250-500,500-100...10,000-11,000. Is there a method that takes into account the number of users that have a specific reputation wen specifying the bin limits or can I just assign arbitrary bins? slotishtype


title: Any transformation or technique to better represent data clumped in the bottom left corner of scatter plot?

text: I have a scatter plot where the data are distributed very unevenly. If you can imagine the ares of the scatter plot divided into 4 quadrants, about half of the data (i.e., half of the data points) like in the bottom-left or southwest quadrant. Apart from deleting the other half of the data as "outliers", is there a transformation or technique to represent this distribution? Edit: I am looking for a general rule-of-thumb -- or heuristic if the $10 word makes you feel better -- similar in concept to the type of transformation heuristic that textbooks provide and practitioners use when faced with non- normally distributed data. The point -- like distribution transformations -- is to get the data to act more linearly, because many parametric statistical techniques are linear "under the hood". Thanks in advance!


title: Correcting for sampling in arbitrary distributions

text: Let's say I've got a series of Gaussian random variables $^{1}x_{i} \sim N(x_i,\sigma_i)$, with each $\sigma_i$ different, but known _a priori_. I've got another identical series of Gaussian random variables $^{1}y_{i}$, which may (or may not) be correlated to the set $^{1}x_{i}$. Both series of random variables are drawn from a population, which here I'll call '1'. I've got a finite number of each $x_i$ and $y_i$, and the cardinality of $x_i$ and $y_i$ is the same. I then go away and draw a similar (but independent) series of random variables from population 2; $^{2}x_{i}$ and $^{2}y_{i}$. The cardinality of the set $\\{^{2}x_i\\}$ is not necessarily the same as $\\{^{1}x_i\\}$. I then wish to calculate an arbitrary -- but known -- function of the $x_i$, say, for example, $f(x,y)=\arccos\left(e^{-x_i/y_i}\right)$. It's clear that, in general, $f(x,y)$ is not normal, but has a distribution that can easily be numerically determined (e.g. by Monte-Carlo methods). My question is this: how do I determine if the expectation value of $f(x,y)$ is different for population one as opposed to population two? In other words, what's the best way to do hypothesis testing with the most power when one has a known, but not-normal sampling distribution?


title: Detecting poaching by age class proportion - suggested test for detecting decrease in given age class?

text: My dataset consists of observations of individuals, recording the date and the individual's age class. Age classes vary in bin size:               0-4     5-9     10-19     20-40      Example data would be:               Date       Age class     01/01/10   0-4     01/01/10   10-19     02/01/10   20-40     04/01/10   0-4     04/01/10   5-9     04/01/10   5-9      Not every day in a given month is sampled, and not every month in a given year is sampled (due to lack of observers), though I have data from each month of the year for the years 2005, 2006, 2007, 2008, 2012, and 2013. I have roughly 3000 data points from each of these years. I would like to detect a decrease in the observations of a given age class, e.g. has there been a decline in the number of individuals from the 20-40 age class observed. To avoid problems of sampling effort, it seems sensible to use proportion of total individuals observed from a given class, e.g. has the number of individuals from the 20-40 age class decreased as a proportion of the total number observed each year from 2005 to 2013. This could be daily, monthly, or yearly proportions, though sometimes very low numbers are observed in a given day, which might do funny things to proportions. For a bit of added complexity, there is seasonal (i.e. intra-year) variation in proportion of a given age class, natural and unrelated to poaching or a real decline in a given age class. Ideally a test that uses as much data as possible would be used, though if necessary, the years 2005-2008 and the years 2012-2013 could be grouped and tested against each other. I would appreciate any advice on suitable tests I could use on this data. I'm currently manipulating it in Excel, but could use R if necessary.


title: What differences and relations are between common random numbers, antithetic variates, and control variates

text: common random numbers, antithetic variates, and control variates are all variance deduction methods. From their Wikipedia articles, it seems to me that they are all used when estimating some quantity of two distributions, and they reduce variance of estimators by controlling the covariance between the random variables of the two distributions. So I wonder what the differences and relations are between common random numbers, antithetic variates, and control variates? Thanks! * * *


title: How to determine if short strings of text are closely related to a larger text?

text: I have 1 short string of text (let's say it's a tweet, max 140 characters):   `"A review of my beloved Roku 3 media player"` I also have a larger body of text (like a blog article, hundreds of words) which I know is related to the tweet:   `"The Roku 3 media player is a great way to watch your favorite ...."` The tweet and blog article are both about the **Roku 3 media player** specifically. Same author, and they share many of the same phrases, words, collocations, etc. It's likely that the string "Roku 3" appears in the text, along with variations like ""Roku 3 streaming media player", "Roku 3 player" etc I then have 10 other tweets, some which are related to the "Roku 3 media player", some of which are not (but very similar): RELATED "A good review of the **Roku 3 media player** "   UNRELATED "The Roku 2 media player review"   RELATED " **Roku 3** is amazing"   RELATED "The **Roku 3** is better than the Roku 2 by far"   RELATED "The **Roku** version **3** streaming media player, fully reviewed"   UNRELATED "A comparison review of the top **3** media player boxes. **Roku** , Android, Toshiba"   RELATED " **Roku 3 streaming media player** reviewed" Those are some examples, and I would have 10 in total. All of the tweets contain "Roku", some are about "Roku 2" and are unrelated, one uses "Roku version 3" and is related etc. Obviously, this is a very small data set. What is the best method to classify each of the 10 tweets as relevant or not, in relation to the first tweet and blog article? What sort of features would be useful?


title: Pairwise deletion in multiple regression

text: Approximately 50% of cases are missing data on one of my predictor variables. With the default option selected (listwise treatment of missing data), the models produced are weak. This is probably because the listwise option reduces _n_ substantially. The alternative (pairwise exclusion), when selected, produces a strong model (the total variance explained is about 50%) with a number of significant predictors (the variable with 50% missing data is a significant predictor in this model). However, this sounds a bit too optimistic. I've read that when pairwise exclusion is selected, SPSS will base degrees of freedom for significance testing on the number of cases with complete data (in this case, 32) rather than on the total number of cases. From what I understand, this means that the significant effects may be exaggerations. Am I right to be concerned about the potential for exaggerated effects when pairwise exclusion is selected? Or are the parameter estimates (and the model as a whole) still trustworthy?


title: Plotting data from several files on one plot

text: I have more than one file of data from different experiments. For the sake of the argument here, let's say I have three files `E1.txt`, `E2.txt`, `E3.txt`. I anticipate that I will have more files in future, say `E4.txt`, `E5.txt`, etc. Each of these files has two columns, Time(h) and Growth. Plotting data from one file, for instance, E1.txt, is simple. I just read the data from `E1.txt` using the `read.table` function in R. I then assign the data from `E1.txt` to variable `E1`. Then I just use the plot function this way, `plot(E1)`. Now, what do I need to write in a script that will plot all the data from the `E1.txt`, `E2.txt`, and `E3.txt` files on the same plot? I thought of writing `plot(c(E1,E2,E3))`, but it doesn't work. Thanks for the help.


title: Seasonal differencing in Arima function in forecast package in R

text: I just want to ask about the `Arima` function in forecast package. The usage of it is,               Arima(x, order=c(0,0,0), seasonal=list(order=c(0,0,0), period=NA),         xreg=NULL, include.mean=TRUE, include.drift=FALSE,          include.constant, lambda=model$lambda, transform.pars=TRUE,          fixed=NULL, init=NULL, method=c("CSS-ML","ML","CSS"), n.cond,          optim.control=list(), kappa=1e6, model=NULL)      My data has trend and seasonality, so I applied seasonal differencing using the codes below,               Diff1LogCP <- diff(LogChickenProd, lag = 4, differences = 1)      Now, I want to have a model SARIMA(1,1,1). Using the `Arima` function, I'm not sure with the `order` and `seasonal` argument. Since my model is SARIMA(1,1,1), so what I did is,               SARIMA111 <- Arima(ChickenProd, seasonal = list(order = c(1,1,1), period = 4))      Notice that I ignore the `order` argument before the `seasonal` argument. This is because I'm using a seasonal ARIMA model, and not the ordinary ARIMA. So I use the `seasonal` argument only. Is that correct? Or do I need to include the `order` argument, that makes my new codes be               SARIMA111 <- Arima(ChickenProd, order = c(1,1,1), seasonal = list(order = c(1,1,1), period = 4))      Or should it be like                SARIMA111 <- Arima(ChickenProd, order = c(1,0,1), seasonal = list(order = c(0,1,0), period = 4))      Thank you in advance!


title: Generalizing Add-one/Laplacian Smoothing

text: Let us assume we are estimating a proportion or rate of "hits". If we have $h$ hits and $m$ misses, the obvious estimator is $\dfrac{h}{h + m}$ In order to avoid unreasonable estimations of $0$ or $1$ when our sample size is small, we can do some (Add-1/Laplacian) smoothing: $\dfrac{h+1}{h+m+2}$ I have read that this has a Bayesian interpretation of having a 50/50 prior over the hit-rate. A couple ideas about generalizing this spring to mind, but I'm uncertain as to the theory.   1. ### Level of confidence If I'm very confident that the hit rates are 50/50, I could add $2$ instead of $1$ to the hits and misses. Or if I'm less confident, I could add $1/2$. What doesn't immediately make sense to me though is what the Bayesian interpretation (if any) is. Isn't the prior just $p = 0.5$, and that's that? Or is there a natural way to represent concentration? If so, what level of concentration does the Laplacian add-one smoothing correspond to? If not, why doesn't this variable-confidence scheme make sense?   2. ### Different prior probabilities Instead of a uniform prior, we could have some other prior over the hit rate. To accomplish this, we could add some number to the hits and some number to the misses such that the proportion worked out. However, I don't immediately know how to parameterize it. For instance, if I have a prior of 1/4, should I add $0.5$ and $1.5$ to the hits and misses, or should I add $1$ and $3$? This ties into the previous question about level of confidence. I'd like to parameterize this so I can change the prior probability without altering the confidence (if such a concept makes sense).


title: Sample size for counterbalanced design

text: I would like to run an experiment with 30 participants in a treatment group. Each participant will work in each of three different conditions: two weeks in 1st condition, two weeks in the second condition and two weeks in the third one. Is my sample too small? Let's say the expected effect size might be about 0.80. Will it also have to be taken in consideration.


title: Only allow positive MCMC-samples: Implications for credible interval

text: We calculate a Bayesian model and only expect positive values for our parameters. Our prior is however a uniform prior—we get negative samples from MCMC. For Bayes-factor calculations we use the Savage-Dickey density ratio. In order to account for the restriction to positive values, we normalize the positive area of prior and posterior distribution so they each have an area of one and calculate the density ratio. _Question:_ what does this mean for our credible interval? At the moment we use the interval calculated by JAGS. But JAGS doesn't know about our restriction. I suppose the credible interval changes when discarding all negative samples? How should we calculate the correct credible interval? Thanks a lot!


title: How to compute expectations from a probability density function?

text: How to find a tax/subsidy in an income probability density function situation? I am asked the following question: > Suppose all families with $Y \lt 20$ are given transfer payments equal to > $(20-Y)/2$; and suppose the cost of these transfers is raised by taxing all > families with $Y \gt 20$ an amount $t(Y-20).$ > > Solve for $t.$ I am trying to understand how to develop an equation to solve for $t,$ but none of my answers work. The probability density function is triangular: ![income probability density function](http://i.stack.imgur.com/aaXda.jpg) * * * The values I have solved for are:   * **Height** : 1/30   * **Mean** : 26.67   * **Median** : 25.36   * **slope of $Y\gt 20$** : 1/600   * **slope of $20\lt Y\lt 60$** : 1/1200 Any help would be greatly appreciated.


title: Standard deviation of revenue where numbers are only positive

text: Is there a way for standard deviation to ignore negative numbers? In other words, how can you make it so that 2 or 3 deviations do not extend into negative territory? When looking at a set of only positive numbers, within one or two standard deviations, I end up with negative numbers. It's essentially impossible for a data point to be negative. Am I doing something wrong? Also, when using deviations from the mean, do you add or subtract the full deviation from the mean or divide by two and plus and minus that to get one deviation from mean I'm sure this embarrassingly simple for this group but this is not my strong point! Please let me know if I can clarify. My recent experience with statistics is relatively limited. Thanks!


title: Visually random model residuals, yet heteroskedastic? ( very small Breusch-Pagan Test P-Value)

text: Can anyone explain why the BP, Breusch-Pagan, test rejects homoscedasticity with such an apparently randomized plot of residuals? ![enter image description here](http://i.stack.imgur.com/VaaQi.png) ![enter image description here](http://i.stack.imgur.com/hJTcw.png)


title: Why do political polls have such large sample sizes?

text: When I watch the news I've noticed that the Gallup polls for things like presidential elections have [I assume random] sample sizes of well over 1,000. From what I remember from college statistics was that a sample size of 30 was a "significantly large" sample. It was made to seem that a sample size over 30 is pointless due to diminishing returns.


title: Computing time issues when using ordinary kriging

text: I am having time issues when using ordinary kriging. I have a spatial grid of 200x300 where I want to calculate the interpolated values. For that I am using ordinary kriging. Now for each point in the grid, I calculate the variogram based upon a certain window size. So, instead of creating a variogram from the whole available dataset I use a window. I do not intend to use the whole dataset since it is of size 10,000 and even for calculating distance pairs it will consume around 10,000*9999/2 location. So, I intended to use the windowing method. After I create a variogram, I again use another window smaller in size to define the neighbors of the point where the data is to be interpolated. Then I calculate the weights for these neighboring points using the variogram of the larger window. I have managed to get the interpolated value. However for each cell it takes around 0.04 second so for a grid of 200x300 it takes around 200x300x.04 i.e .67.hrs. I have around 400 sets of grids. So it will take me around 268 which is not acceptable. Is there any other efficient way to accomplish this?


title: How are the constraints of the parameters of a GARCH model calculated?

text: For a simple GARCH(1,1) model: h(t) = K + G*h(t-1) + A*e^2(t-1), The parameters need to be constrained such that G + A < 1\. Can anyone please tell me how this is derived since I've seen this in a lot of books but there's not much explanation given about this. Thanks!


title: Find most likely fit for subgroup size and average density from product and supergroup size

text: I need to find the most likely sets of _group subset percent_ and _frequency_ for a _given result_ having the _super group size_. Let me illustrate with an example. Assume I have the following:                Entity | Product | Price | Total Sales | Entity Customers     -----------------------------------------------------------      BusA   | Yo-Yos  |    $6 |    $100,000 |         106,000      BusB   | Yo-Yos  |    $4 |    $173,000 |          72,000      BusC   | Yo-Yos  |    $7 |    $110,000 |          96,000      BusD   | Yo-Yos  |    $9 |    $ 80,000 |          55,000      Having only this, my charge is to back into the most likely sets of:   1. Percent cohort buying yo-yo's out of total entity customers.   2. Average yo-yo's purchased at each entity's price. When I model a function, I end up with:               f(percentBuyer, purchaseDensity) = {       totalSales / (         totalCust * percentBuyer * purchaseDensity       )     }      `purchaseDensity` representing an average per buyer. I'm seeking for the best fits against a test. A sample set:               [ 5%, 1],[ 5%, 2],[ 5%, 3],     [10%, 1],[10%, 2],[10%, 3],     [15%, 1],[15%, 2],[15%, 3]     ... et cetera ...          And so on. My sense is that I can graph these against the `Price` and `Total Sales` in the different rows and which values per set fit best the `Price` at that entity, I'll have my best fit. I know I could be getting something fundamentally wrong up to this point; I hope not. The issue I run into is I'm unsure how to graph this where the best fits can be determined and a range for best fit be developed that present a reasonable prediction of characteristic fluctuation. Where do I go from here? Go back? Am I overlooking a direct way to handle this type of concern? * * * Note, this is **not** homework; it is related to a project I am assigned at work, and involves a risk management cost-benefit depiction. If there are more appropriate tags for this, please feel free to edit them in.


title: If one group has 23 participants and the second group has 117 participants in it, can we rely on the result of t-test?

text: I was evaluating a research paper in which the authors have 23 male participants and 117 female participants. They have applied t-test to calculate Gender differences and have even concluded on significance of the difference. I am not sure when the groups are so markedly difference in terms of their size, we can use it for drawing such conclusions. Please help me!!!


title: Power analysis for observational study of secondary data

text: I am involved in a project that is attempting to measure the impacts of person-level, county-level and medical service area-level predictors on individual probability of returning to a hospital within a year window. I see this as a 3 level random effect logistic regression model. There are lots of predictors, probably 10 at the person level, 10 at the county and probably 5 at the hospital level. There are nearly 300,000 individuals nested withing 3000 counties and 154 hospitals. I would think that I should have good statistical power to detect differences in my outcome. My question, is: since the sample is already determined and I'm not collecting data (this is all secondary data), is a power analysis even appropriate? If it is appropriate, but my PI can't specify the effect size they are looking for, then what? I would appreciate any input the list may have on this problem. Best, Corey