title: Different p-values for anova() and summary() in glm.nb

text: using this code               library(MASS)     n = c(300, 200, 100)      group = rep(1:3, n)     x=c(rnbinom(300, size=2, mu=2.47), rnbinom(200, size=2, mu=2.27),     rnbinom(100, size=2, mu=2.27))     glm1 = glm.nb(x ~ factor(group))      I have created three neg. binomial distributions and analysed them using a glm. The independent variable (group) is meant to be categorial. When I use `anova(glm1)`, I get an p-value (Pr (>Chi)) of 0.07. However, when I use `summary(glm1)`, I get a p-value (Pr (>|z|) of 0.2 for factor 2 and 0.02 for factor 3. Thus, depending on if I use `summary()` or `anova()`, the factor is significant or it's not. I have three questions and would be very glad if someone could help me with one of these.   1. Can I conclude if group is a significant factor or not?   2. Residuals are not distributed normally. Can I use `anova(glm1)` nonetheless?   3. Apparently `summary()` gives me an intercept and two factors. Do the two factors relate to the difference between group 1 and 2 respectively between group 1 and 3? I would really appreciate your help, even only a small hint on one of these question would be great. Please tell me, if any information is missing or if I didn't express myself comprehensibly.


title: longitudinal correlation in SPSS

text: My clinical study aim is looking at treatment effect on two outcomes. Each subject underwent 3 measurements at 3 different time points at baseline, 6 and 12 months after initiated treatment. We measure two dependent variables at each time point. I would like to perform correlation analysis to see if there is any correlation between the two dependent variables in each subject using baseline variable for each subject as reference. I am wondering if Spearman's rank correlation is appropriate for this.


title: Best K in K-fold cross validation

text: I'm using `K-fold cross validation` technique for generating `train`,`test` and `validation`indexes for neural network. My sample size is `230~700`. What is best `K` for cross validation here. Now I'm using `10-fold cross validation` but i think it is too high. what is your idea? Thanks.


title: Some questions on Principal Component Analysis (PCA)

text: I am trying to understand some descriptions about PCA (the first 2 are from Wikipedia): 1) Principal components are guaranteed to be **independent** only if the data set is **jointly normally distributed**.   * Is the independence of principal components very important? How can I understand this description? 2) PCA is sensitive to the relative **scaling** of the original variables.   * What does 'scaling' mean there? Normalization of different dimensions? 3) The transformation is defined in such a way that the first principal component has the largest possible variance and each succeeding component in turn has the highest variance under **the constraint that it be orthogonal to the preceding components**.   * Can you explain this constraint? Thank you for your help in advance...


title: modeling error when estimating the parameters of the beta distribution

text: After reading _Basketball Beta and Bayes_, I started to think about parameter estimation with observational error. If you are counting the number of successful and failed free-throws how can you account for score keeping errors? I have a very similar problem where in addition to a set of outcomes (success or failure) I have a likelihood that the observation was correct. It would analogous to the free-throw example if instead of seeing the free-throws you received two likelihoods (score or miss). How can I utilize this type of data to estimate the parameters of the beta distribution? Ideally the beta distribution would be almost uniform if the likelihoods indicate a high probability of error.


title: what is difference between ordinary least squares and residuals

text: I think the ordinary least squares are the sum of the vertical distance between the observed data to the model line(regression). And residual is calculated by add up all vertical distance between each observed data and the model line(regression), not by sum of squares? I am still not clear about these concept. Thanks!


title: What is the best chart to display an entity's magnitude?

text: I am looking to display article ranking through a bubble like chart, which is focused on magnitude, so the x and the y axis are irrelevant. I have looked at d3 but didn't find my chart type. Essentially, I want the bubbles to float around lightly, just some animation. Any data-viz experts aware of such a chart type/library that I could find useful in this regard. Thank you


title: Summarizing a lognormal distribution with geometric mean and standard deviation

text: I have some data that I strongly suspect are lognormally distributed, and I'd like to summarize the distribution using the mean and standard deviation. I've read that with lognormal distributions the goemetric mean and standard deviation should be used, but using them produces slightly strange results. Specifically, when using the _sample_ mean and standard deviation I calculate that 86% of my data lie within +/- 1 standard deviation about the mean. However, when I use the geometric mean and standard deviation I calculate that only ~5% of my data lie within +/- 1 standard deviation about the mean. I have over 30,000 pieces of data, and the distribution looks strongly lognormal, but I don't really understand these results. Is it valid to use the geometric mean and standard deviation here?


title: Is cross-validation the most important measure of a predictive model's effectiveness?

text: Why bother with p-values, R squared, etc. ... Model size is not a factor with the computing power available now so why not just run multiple iterations of all possible sets of input variables and see which one has lowest cross- validation error?.


title: Determine which group is more important

text: I have a question how to analyze mean differences where I am unsure whether I am lacking the statistical vocabulary or simply ignoring something important. I have a 5 point satisfaction variable (very satisfied to totally unsatisfied) as dependent variable and three group independent variable to predict the satisfaction. The independent variable has three groups: "Top 2 users", "Neutral", "Bottom 2" - these groups are defined according to some other detailed satisfaction questions which are aggregated before (e. g. "Do you like your telephones color? Do you like the camera? etc.). I have used a Oneway test to determine whether or not there are differences in the first place - and the differences are highly significant and very obvious (mean in the Top 2 group at 1.9 and in the Bottom 2 group 4, ignoring the "neutral" group which is not of real interest). But know my supervisor asked me which of the two groups is more important - e.g. is it more important to gain Top 2 users or to avoid having Bottom 2 users. I am unsure if I can answer the question with the data I have at hand. Do you have an idea if this question is even answerable? I have calculated the effect size, which is quite strong with -0.76, but is this really a good answer to the question at hand? My method to get the effect size: df is a dataframe, var is the independent variable, testvar the dependent variable.               getEffectSize <-function(df, var, testvar)     {       out <- by(df[,testvar],df[,var], stat.desc)            effectSize <- mes(out[[1]][["mean"]], out[[3]][["mean"]], out[[1]][["std.dev"]], out[[3]][["std.dev"]], out[[1]][["nbr.val"]], out[[3]][["nbr.val"]])           print(effectSize)     }


title: Categorical inputs in a linear regression

text: if I use values (1,0) for a categorical column in X... does it lead to any ill conditioned matrices during OLS/ robust regression?.


title: Approximating a Complementary Cumulative Distribution Function via a piece-wise function

text: I hope this is not too much to read, but I tried to give you a specific overview over my problem. I am currently trying to model the German electricity market, with a special focus on balancing power. Considering this setting I have a data set consisting of hourly observations of balancing power in Megawatts for the year of 2013. As the application of this data set would not represent any uncertainty in the data I wanted to use a little trick to "assume" uncertainty. While this is really only a basic approximation my model has to have some sort of "uncertainty", else there will be a big mistake. In order to do this, I divided my data set into five blocks (too many more would be to much for my calculation) of varying Megawatt-steps calculating the "probability" of a specific value to be in one of these blocks using the Complementary cumulative distribution function. What I mean by that is that I took the aforementioned function and tried to best approximate it in a (5-)piece-wise fashion. Unfortunately I had to do this in my example by a "sense of proportion". I hope the picture makes it a little more clear (On the x-axis the values are in Megawatts). **EDIT:** To be clear, I am not trying to best-possibly approximate the whole- function in a sense of a non parametric boundary estimation as user603 kindly described in his answer. (Or maybe this can be done using non parametric boundary estimation...?) I just want to fit the five (or whatever number) blocks best-possible to the function in a fashion that, e.g. the Sum of squared Errors is at a minimum. Also: i would love the knot location to be estimated. ![\(5-\)Piece-wise approximation of the CCDF](http://i.stack.imgur.com/xh9Q2.png) The following code example should further explain what I do:               # This function divides the input data into 5 different blocks.     # v, w, x and y are all Megawatt values dividing the data set     getCall_block <- function(v, w, x, y, data){            n <- nrow(data)       #' parameters for function       data_v <- subset(data, data[,2] <= v, select = c(colnames(data)[1],                                                        colnames(data)[2]))       data_w <- subset(data, data[,2] > v & data[,2] <= w, select = c(colnames(data)[1],                                                                       colnames(data)[2]))       data_x <- subset(data, data[,2] > w & data[,2] <= x, select = c(colnames(data)[1],                                                                       colnames(data)[2]))       data_y <- subset(data, data[,2] > x & data[,2] <= y, select = c(colnames(data)[1],                                                                       colnames(data)[2]))       data_r <- subset(data, data[,2] > y, select = c(colnames(data)[1],                                                       colnames(data)[2]))       #'        df_v <- getAverage(data, data_v, 0, v)       df_w <- getAverage(data, data_w, v, w)       df_x <- getAverage(data, data_x, w, x)       df_y <- getAverage(data, data_y, x, y)       df_r <- getAverage(data, data_r, y, max(data[,2]))            return(rbind(df_v, df_w, df_x, df_y, df_r))     }      And the getAverage function...               #This function is called in the above mentioned function.      #It calculates the mean average of the points of the CCDF lying in one block.      #v and w are the input parameters in Megawatts.      #They can be understood as starting point and end point of a given block.          getAverage <- function(data_c, data, v, w){       #' Imposing greater equal on v if v > 0       if(w > v){         if (v > 0){           v <- v + 1         } else {           v <- 0         }       #' Initializing div, df, mean, n, sum         div <- 0         mean <- 0         n  <- nrow(data)         sum <- 0         df <- data.frame(ID = data[,1], call_freq = numeric(n))       #' Creating inverse ECDF based on dataset         F1 <- ecdf(data_c[,2])       #' Building sum over every data point in the interval [v;w]/(v;w]         for (i in seq(v, w, 1)){           sum <- sum + (1 - F1(i))           div <- div + 1         }         mean <- sum / div         df$call_freq = mapply(function(x) return (x),  mean)         return (merge(data, df, by = "ID"))         } else {           cat("Check input parameters. v is greater than w.")       }     }      My questions are the following: I) Does this way of approximating a CCDF sound okay for you guys? What did I overlook or completely mess up? II) Is there a way to automatically approximate the block sizes, in order to not rely on a sense of proportion. E.g. minimizing the sum of squared errors between the piece-wise function and the original function. Maybe this could be done by using a clustering algorithm like k-means? But then how to cluster this CCDF...? III) Is there a far more easier way to approximate any function in a piece- wise fashion? Thanks in advance!


title: Difference between disease conditions by non-parametric ANOVA

text: I've the following simplified data as shown.               PATIENT B1  B2  B3  B4     A1  0.19    0.18    0.12    0.85     A2  2.8 0.24    0.06    0.11     A3  0.03    0.24    0.24    0.07     A4  0.65    1.6 0.08    0.31     F1  0.17    0.07    3.86    2.41     F2  0.11    1.74    0.51    0.34     F3  0.11    2.28    0.57    4.06     F4  0.23    0.68    2.51    0.31     S1  0.5 0.19    2.13    0.09     S2  1.21    0.25    2.02    0.2     S3  2.06    0.05    0.16    0.4     S4  0   1.02    0.01    0.37     J1  2.64    0.68    0.1 0.3     J2  2.7 3.89    0.15    0.34     J3  0.09    0.22    0.17    2.74      We have measurements for several biomarkers B1, B2, B3, B4 and for different disease subtypes A, F, S etc. for 4 (3 for J) different patients in each disease group A1-A4, F1-F4, etc. The data is non-Gaussian so I wish to do a non-parametric Kruskal Wallis test to show if the different disease subtypes are different or not. I've consolidated the data in the following format (Pastebin Link) throwing away the individual biomarker information and did a Kruskal Wallis test. Is there any way I can use the Biomarker value in the Kruskal Wallis given the biomarker info as shown in this Pastebin link. Can a Friedman Test be done for this case?


title: Predicting the next value

text: I've recently read a blogpost where someone tries to predict the outcome of eurovision. Quoting from the summary of the site: > Essentially, we can look at people’s voting preferences in the Eurovision > Song Contest as composed of two components: **song quality, and a > “friendship” score** , which takes into account how much the voting country > likes or dislikes1 the country being voted on. If we want to know whether a > voter V will rank country A or country B higher, we add up the song quality > and the friendship score in each case, and subtract the two. Then we can > take this difference, feed it through a logistic curve and use the result as > a probability. > > I’ve taken voting results from both the Eurovision finals (going back to the > introduction of televoting in 1998) and the semi-finals (going back to their > introduction in 2004). I’ve then used a * _Markov Chain Monte Carlo sampler_ > *2 to calculate the song qualities and friendship scores, assuming that > they’re both normally distributed. > > Once I’ve got the parameters, it’s relatively straightforward to **run a > simulation** of this year’s contest, including the semifinals and all the > voting procedures. Last year I ran 10,000 simulations like this, and looked > at both who the likely qualifiers from the semifinals were, and the overall > winner. The model managed fifteen out of twenty of the qualifiers, as well > as the eventual winner. Although this level of success, particularly > predicting the winner, was probably mostly luck, this year I’d like to do > the same. I was wondering if the same kind of approach can be used in order to solve my problem. I want to be able to predict the next state based on the previous state(s). There are 7 states. I do not know what the transition matrix is. I know that the state space consists of all permutations from 1..7 . I know I can probably do this with a HMM. Though I was wondering if this is possible with MCMC (as used in the blogpost). I have data, which are examples of these state transitions, though I'd like to see if it is possible to predict the next state based on my data with mcmc. My data for example looks like this:               6 ==> 1 ==> 3 ==> 5 ==> 2 ==> 6 ==> ...      And I'd like to predict the next value after the 6. ( Is this possible with the same approach like in the eurovision example ?) Additionally I can add extra data if necessary. This data would indicate which day of the week (1 to 7) and another parameter for the week of the year( 1 to 52). So then the data would look like this:               state: 6 | state: 1 | state: 3 | state: 5 | ...     day : 1  | day: 5   | day : 7  | day: 2   | ...     week : 3 | week: 3  | week : 3 | week: 4  | ...      Maybe this kind of data is better suited ?


title: How to construct a 95% confidence interval of the difference between medians?

text: My problem: parallel group randomized trial having a very right-skewed distribution of the primary outcome. I do not want to assume normality and use normal-based 95% CIs (i.e. using 1.96 X SE). I am comfortable expressing the measure of central tendency as the median, but my question is then how to construct a 95% CI of the difference in medians between the two groups. The first thing that comes to mind is bootstrapping (resample with replacement, determine median in each the two groups and subtract one from the other, repeat 1000 times, and use the Bias-corrected 95% CI). Is this the correct approach? Any other suggestions? Thanks.


title: When to aggregate for time series forecasting?

text: I have a some historical sales data for various product SKUs, including category information ("department" "category", "subcategory"). I want to use this to generate sales curve (a baseline forecast for future demand), no doubt using an appropriate exponential smoothing algorithm. The trouble is that the data for a particular SKU will often be too sparse to be able to infer a sales curve. I'm thinking that when that is the case, the sales data for the SKU could be aggregated with other SKUs for the same product and I could then use the sales curve for the product rather than the SKU (which would presumably be more reliable). If the product doesn't have enough sales history then I'd use the sales curve for the subcategory, and so on and so forth… potentially even having to resort to a sales curve for the top level category (a department in this case) if insufficient sales data is available to generate reliable sales curves for lower level categories. What I'm wondering is, how do I determine when I have enough data to be able to trust a sales curve? How many sales "data points" do I need for a SKU or a subcategory in order to be able to trust its sales curve? Perhaps I'm approaching this the wrong way though. An alternative that I'd considered was to use data from all levels of the hierarchy in the forecast with appropriate weightings for each. So a product's future sales would be predicted by `W1 * Forecast(Product) + W2 * Forecast(Subcategory) + W3 * Forecast(Category) + W4 * Forecast(Department)`.... Working out the weightings might be a bit complex but not impossible (and might be done using a regression analysis). I'm sure this kind of problem has been solved before but I'm a bit of a stats newbie. Are there "off the shelf" algorithms/techniques that could use for this problem?


title: How to output varimax-rotated data in R?

text: I have ran PCA on 25 variables and selected the top 7 PCs using 'prcomp'.               prc <- prcomp(pollutions, center=T, scale=T)      I have then done varimax rotation on those components.               varimax7 <- varimax(prc$rotation[,1:7])      And now I wish to varimax rotate the PCA-rotated data (as it is not part of the varimax object - only the loadings matrix and the rotation matrix). I read that to do this you multiply the transpose of the rotation matrix by the transpose of the data so I would have done this:               newData <- t(varimax7$rotmat)*t(prc$x[,1:7])      But that doesn't make sense as the dimensions of the matrix transposes above are 7*7 and 7*16933 respectively and so I will be left with a matrix of only 7 rows, rather than 16933 rows... does anyone know what I am doing wrong here or what my final line should be? Do I just need to transpose back afterwards?


title: How can I adapt ANOVA for binary data?

text: I have four competing models which I use to predict a binary outcome variable (say, employment status after graduating, 1 = employed, 0 = not-employed) for n subjects. A natural metric of model performance is hit rate which is the percentage of correct predictions for each one of the models. It seems to me that I cannot use ANOVA in this setting as the data violates the assumptions underlying ANOVA. Is there an equivalent procedure I could use instead of ANOVA in the above setting to test for the hypothesis that all four models are equally effective?


title: Probability of period without event

text: I have a data set of a list of invoices each of which have a date. I'm trying to detect when I might consider that a customer has stopped ordering a part. The way I'm approaching this might be completely wrong. I'm taking the invoices, ordering by date, and then counting the difference between dates of consecutive invoices. I then produce a nice graph and it lets me say something like "95% of all orders were within 15 days of each other". I don't know where to go next :) Is this concept at all related to the Pareto distribution? Reading up on the Poisson distribution I can take the same data and produce the probability of N invoices per day. I get something like 57% of days have 0 invoices, 23% have 1 invoice, 12% have 2 invoices, etc. At this point I get to the not knowing what I'm doing part :) The end goal is to be able to notice stopped orders, ideally with some kind of tweakable 'confidence' cutoff. I'm afraid of using the 'confidence' word because I assume to be using it wrong. Thanks!


title: In R, given an output from optim with a hessian matrix, how to calculate parameter confidence intervals using the hessian matrix?

text: Given an output from optim with a hessian matrix, how to calculate parameter confidence intervals using the hessian matrix?               fit<-optim(..., hessian=T)          hessian<-fit$hessian      I am mostly interested in the context of maximum likelihood analysis, but am curious to know if the method can be expanded beyond.


title: What is an appropriate method for providing bounds when performing maximum likelihood parameter estimation?

text: I have successfully implemented a maximum likelihood estimation of model parameters with bounds by creating a likelihood function that returns NA or Inf values when the function is out of bounds. I optimize the function using optim in `R`. detailed example available on github Quick example:               likelihood.fun<-function(par, ...){            likelihood<- -sum(dnorm(..., log=T))            if(any(c(par[1]<0,              par[1]>5,               par[2]>5,              par[2]<0)){likelihood<-NA}            return(likelihood)          }      Is this equivalent to box optimization or deprecated compare to box optimization? If this is not equivalent: How can I implement this using               optim(..., method="L-BFGS-B", lower=c(...), upper=c(...))               from example, this does not seem to work:     optim(..., method="L-BFGS-B", lower=c(0,0), upper=c(5,5))      or               constrOptim()      This is linked to this question on constrOptim.


title: How to set limits using constrOptim in R?

text: I am using constrOptim to minimize a log likelihood function for maximum likelihood estimation of parameters. I wish to set the bounds on my parameters, but to not understand the constrOptim definition of the feasibility region. The feasible region is defined by ui %*% theta - ci >= 0 I have a set of parameters with bounds [lower, upper]               a[0,5]  ie 0<a<5     b[0,Inf]     c[0,Inf]     e[0,1]          theta (starting values) = c(1, 1, 0.01,0.1)      What are the _ui_ (constraint matrix (k x p)) and _ci_ (constraint vector of length k) for these parameter bounds? Is there a straightforward way to get from a list of upper and lower bounds to a _ui_ and _ci_ value?


title: Can I use T-Statistic for constructing a Confidence Interval for sample from large population?

text: Can I use a T-Statistic to construct a Confidence Interval for a sample collected from a _large_ population and without knowledge about the population's standard deviation? I know the "model" answer is that a sample, even if the sample size $n$ is small, would have to use a Normal Z statistics if it was collected from a large population. But I have used a T-Statistic in a problem like this one because earlier, I felt that it would be more appropriate to use a T-Statistic since I'm having the standard mean error estimated by the sample standard deviation. And now, I am worried that I've done wrong. So, can I use a T-Statistic to construct a Confidence Interval for a sample collected from a _large_ population and without knowledge about the population's standard deviation? In my case, $n=8$. If no, why would it be wrong to use a T-Statistic in this case? If yes, also, what are the reasons for being correct to use?


title: Does the presence of an outlier increase the probability that another outlier will also be present on the same observation?

text: ** **Edit:** (10/26/13) More clear (hopefully) mini-rewrites added at the bottom** _**I'm asking this from a theoretical/general standpoint - not one that applies to a specific use case._** I was thinking about this today: **Assuming the data does not contain any measurement errors, if you're looking at a specific observation in your data and one of the measurements you recorded contains what could be considered an outlier, does this increase the probability (above that of the rest of the observations that do not contained measured outliers) that the same observation will contain another outlier in another measurement?** _**For my answer I'm looking for some sort of theorem, principle, etc. that states what I'm trying to communicate here much more elegantly. For clearer explanations see Gino and Behacad's answers._** Example: Let's say you're measuring the height and circumference of a certain type of plant. Each observation corresponds to 1 plant (you're only doing this once). For height, you measure:               Obs 1  |   10 cm     Obs 2  |   9 cm     Obs 3  |   11 cm     Obs 4  |   22 cm     Obs 5  |   10 cm     Obs 6  |   9 cm     Obs 7  |   11 cm     Obs 8  |   10 cm     Obs 9  |   11 cm     Obs 10  |   9 cm     Obs 11  |   11 cm     Obs 12  |   10 cm     Obs 13  |   9 cm     Obs 14  |   10 cm      Since observation 4 contains what could be considered an outlier from the rest of the data, would the probability increase that measured circumference also contains an outlier for observation #4? I understand my example may be too idealistic but I think it gets the point across...just change the measurements to anything. # _**Edited in attempts to make more clear:_** (10/26/13) **Version 1 Attempt of Abbreviated Question:** > In nature and in general, is there a tendency (even a weak one) that the > **_probability is greater_** that the "degree of variance from the mean in > any attribute(s) of an observation" will be _similar_ to the "degree of > variance from the mean in any other* specific attribute of that same > observation" **_in comparison to the probability_** that it will INSTEAD be > _more similar_ to the "degree of variance from the mean in that same* > specific attribute of _any other observation_." _* next to a word means I was pairing what they reference.   "Quotes" used above mean nothing and are used simply to help section parts together/off for clarity._ **Version 2 Attempt of Abbreviated Question:** > In nature and in general, is variance from the mean across observations for > one attribute¹ correlated (even with _extremely_ loose correlation) to the > variance from the mean across observations for all attributes¹? _¹Attribute meaning measurement, quality, presence-of-either of these, and/or nearly anything else that the word "attribute" could even slightly represent as a word. Include all synonyms of the word "attribute" as well._


title: Why is the Tukey confidence interval wider than the T-test confidence interval?

text: I have two confidence intervals for a comparison of two means. The first is a t-test CI and the other is a Tukey CI. The Tukey CI $(2.36, 26.5)$ is wider than the t-test CI $(6.46, 21.33)$. What is the reason for this?


title: Is Cronbach's alpha interchangeable with $\rho$?

text: I am trying to calculate the reliability of my change scores using the following formula. Is $r$ ($\rho$) in this instance interchangeable with Cronbach's alpha? $$ \frac{.5(r_{xx}+ r_{yy}) - r_{xy}}{1 - r_{xy}} $$ Where reliability of difference scores is a function of the score reliability of each measure ($r_{xx}, r_{yy}$) and the correlation between the measures ($r_{xy}$).


title: Power of McNemar Test

text: I have been playing around with the GPower tool to calculate the power of various statistical tests, and what I noticed is that if I can change the proportion of the discordant pair in the McNemar test that I can increase the power of the test substantially. What I wonder is, why can I not delete the discordant pair and carry out the test since the test is only carried out with those anyway? The only change it would have is to change the sample size, but the OR would stay the same, whilst the power of the test could increase. I have not found anything online about this topic yet.


title: $R^2$ in multivariate regression

text: I'm trying to determine how the population $R^2$ value is defined in the multivariate regression model where we have $Y_i = \mu_y + B^\prime(X_i - \mu_x) + err$ Where $Y_i \in \mathbb{R}^q$ and $X_i \in \mathbb{R}^p$ and $q, p>1$


title: Negative values in the data set

text: I have a set of negative values and would like to put a fitted model in R. How would I approach it? Say, I want to use Gamma distribution as a fitted curve using method of moments but the data is all negative and thus, it leads to the negative rate.


title: Correlation of parts of two variables that miminize correlation with another variable

text: I have three data vectors $A$, $B$ and $C$ that are all more or less correlated. What I want is a meaningful definition for a measure of the degree of correlation between those parts of $A$ and $B$ that are not correlated with $C$. I do apologize if this is not fully understandable. I struggle to come to an easier formulation of the problem. However I can give an example: Example: Say $A$, $B$, $C$ could all be binary, i.e. only contain the values 0 and 1. I then would want to know the number of positions where $A$ and $B$ are 1 and $C$ is 0. As a side note: How does one "uncorrelate" variables, i.e. given two variables $A$ and $B$ how to split $A$ into a part that is maximally correlated with B and another that is maximally uncorrelated?


title: Classifier for uncertain class labels

text: Let's say I have a set of instances with class labels associated. It does not matter _how_ these instances were labelled, but _how certain_ their class membership is. Each instancs belongs to _exactly one_ class. Let's say I can quantify the certainty of each class membership with a nominal attribute that goes from 1 to 3 (very certain to uncertain, respectively). Is there some sort of classifier that takes into consideration such a certainty measure and if yes, is it available in the WEKA toolkit? I imagine this situation occurs quite often, for example when instances are classified by human beings which are not always completely sure. In my case, I have to classify images, and sometimes an image could belong to more than one class. If this happens, I give the class a high uncertainty, but still classify it with only one class. Or are there any other approaches to this problem, without a specialized classifier? E.g. only taking "certain" classifications for training? I fear that in this case, there will be more misclassifications because "border" cases are not covered.


title: Plotting variables in transformed space

text: Suppose $A = X_1/X_2$ and $B = X_3/X_4$. Why would one plot the data in $(\log A, \log B)$ space as opposed to $(A,B)$ space?


title: How do residuals affect the t-statistic?

text: Could a t-statistic for a simple regression relationship with a smaller sum of squared residuals but a smaller or less positive slope be larger than that of a plot with a larger sum of squared residuals and a large slope?


title: Why don't asymptotically consistent estimators have zero variance at infinity?

text: I know that the statement in question is wrong because estimators cannot have asymptotic variances that are lower than the Cramer-Rao bound. However, if asymptotic consistence means that an estimator converges in probability to a value, then doesn't this also mean that its variance becomes 0? Where in this train of thought am I wrong?


title: What is my experimental design?

text: I am very new to mixed/multilevel models. I have an experiment where we measured 2 scale variables (varA and varB) in 2 different groups of subjects (control and treatment) at 4 different time points. I suspect that the relationship between varA and varB is different for the 2 groups (control vs treatment), however it should be roughly consistent across the 4 time points. Researching on the internet I've come to conclusion that I need to model this using a multilevel (mixed effects) model. So, I am interested in the difference between the regression line between varA and varB for the 2 groups, but I want to account for the repeated measures. What is my experimental design in this case? It's not very clear to me what are the fixed and random effects. Additionally, is it also possible to test if the regression lines are significantly different across the time point? Would that call for a different model? Thanks very much, any help on how I should go about doing this in R/Matlab/SPSS would be greatly appreciated as well!


title: factor analysis: can irrelevant factors be identified by tweaking FA options?

text: So I ran a factor analysis (principal components method) on a dataset, first using the correlation matrix and then using the covariance matrix, both with varimax rotations. The results of both factor analyses indicate nearly identical first and second factors, and somewhat different third and fourth factors. I am curious if this could be interpreted as factors 3 & 4 not being useful? Or that only the first two factors are appropriate for analysis? Or is this just showing a pattern in the correlation/covariance structures of my data (i.e. correlated correlations/covariances?)


title: Gaussian neighborhood function and non linear learning rate for SOM in R

text: I've been working on SOMs and how to get the best clustering results. One approach could be to try many runs and choose the clustering with the lowest within sum of squared errors. However, I do not only want to initialize random values and have several tries, but also want to choose good parameters. I have read in "Influence of Learning Rates and Neighboring Functions on Self-Organizing Maps" (Stefanovic 2011) that if you do not know which parameters for the neighborhood function and learning rate to choose, it is probably the best option to choose a gaussian function and a nonlinear learning rate. My data is a time series lets say:               matrix(c(sample(seq(from = 10, to = 20, by = runif(1,min=1,max=3)), size = 5000, replace = TRUE),(sample(seq(from = 15, to = 22, by = runif(1,min=1,max=4)), size = 5000, replace = TRUE)),(sample(seq(from = 18, to = 24, by = runif(1,min=1,max=3)), size = 5000, replace = TRUE))),nrow=300,ncol=50,byrow = TRUE) -> data      which has 300 observations with 50 values each. 100 observations each tend to be more similar. I'm working with the kohonen package. The code:               grid<-somgrid(4,3,"hexagonal")          kohonen<-som(data,grid)          matplot(t(kohonen$codes),col=kohonen$unit.classif,type="l")      gives me clusters with values **between 10 to 22** , which is similar to the obersations ![enter image description here](http://i.stack.imgur.com/eNd7N.jpg) I tried the "som" package, too, which offers a gaussian neighborhood function and a inverse-time learning rate.               som<-som(data,4,3,init="random",alphaType="inverse",neigh="gaussian")          som$visual[,4]<-with(som$visual,interaction(som$visual[,1],som$visual[,2]))          som$visual[,4]<-as.numeric(as.factor(som$visual[,4]))          matplot(t(som$code),col=som$visual[,4],type="l")      ![enter image description here](http://i.stack.imgur.com/vaSJV.jpg) Here I get clusters with values **between 15 and 18** , so all clusters "shrink" and get more similar. With different input series I get the same phenomena My two questions:   1. Why do clusters from the self-organizing map with the som package get so extraordinarily similar and shrink to a much smaller range, even though it is said that you get good clusters with gaussian neigborhood function and non linear learning rate?   2. How can I avoid this range shrinking with a gaussian neigborhood function and a non linear learning rate in order to get appropriate clusters?


title: How to combine different feature values to build a score for ranking

text: I am doing some information extraction study on some data. We extract a "word" list from the data. And for each "word", it has a n dimensional feature vector in R ^n. Next, We want to rank the list according to a score. The question is, how can I combine all the feature values to build a score function for my ranking. The difficulty is we can only work in an unsupervised way. I think many people are facing to this question in bioinformatic or other area when doing IE, right? I had thought about max entropy, but it usually works in supervised way. Is there any good method to solve this problem in unsupervised way ?


title: Baum-Welch training example

text: I'm the author of a new Baum Welch trainer using MapReduce for the Apache Mahout project (https://issues.apache.org/jira/browse/MAHOUT-627) I'm looking for an example with a reasonably small data set which I can use to demonstrate the usage of the trainer I have developed. Does anyone here know of any online demonstrations for this training which go beyond the normal "toy" examples? I'm interested in typical HMM training application areas--bioinformatics, NLP, Signal Processing etc.


title: What is the difference between "coefficient of determination" and "mean squared error"?

text: For regression problem, I have seen people use "coefficient of determination" (a.k.a R squared) to perform model selection, e.g., finding the appropriate penalty coefficient for regularization. However, it is also common to use "mean squared error" or "root mean squared error" as a measure of regression accuracy. So what is the main difference between these two? Could they be used interchangeably for "regularization" and "regression" tasks? And what are the main usage of each in practice, such as in machine learning, data mining tasks?


title: Testing the significance of multiple variables in a data set

text: My data set relates to goods sold. I have a unit price, qty sold at that price and a number of other variable, such as sales person, geography, marketing support. the qty sold varies throughout any sampled period of time irrespective of price. It also 'appears' to vary more when the price is changed-as common sense would sugggest. What I would like to measure is the impact/significance of each of the other variables on the qty sold. I'm not a stats expert or coder butr have some experience/knowledge of using linear regression and correlation coefficients in my work


title: Statistical hypothesis testing - how to demonstrate that power $>$ size

text: Let $(X_1,X_2,\ldots,X_n)$ be a random sample from $N(\mu,\sigma^2)$, where $\sigma$ is known. Find the power function of the test of significance of $H_0:\mu=\mu_0$ against $H_1:\mu\neq\mu_0$ . Also, show that power${}>{}$size. My first step is to define the critical region as $W=\{(x_1,x_2,\ldots,x_n):|\sqrt{n}(\bar x-\mu_0)/\sigma|>\tau _{\alpha/2}\}$ where $P[Z>\tau_\alpha]=\alpha$, $Z\sim N(0,1)$, and $\alpha$ is the level of significance. The power function of the test of significance is then, $\beta(\mu)=P[(X_1,\ldots,X_n)\in W|\mu\in\Omega_1]=P[|\sqrt{n}(\bar x-\mu_0)/\sigma|>\tau _{\alpha/2}|\mu\neq\mu_0]$ and the size of the test of significance is, $\sup\limits_{\mu \in\Omega_0}P[ (X_1,\ldots,X_n)\in W|\mu\in\Omega_0]=P[|\sqrt{n}(\bar x-\mu_0)/\sigma|>\tau _{\alpha/2}|\mu=\mu_0]$. How do I proceed from here?


title: Is the number 20 magic?

text: I have references that advised to consider a sample size of minimum 20 for data fitting distribution. Is there any sense in this? thanks


title: Positive, Likely Positives and Unknowns in Random Forests

text: It is my first question here, so my excuses if my contribution is naïve. I'm facing a classification problem in which I have a minority class (~100 samples) labeled as "Positive", another class (~100-1000 samples) labeled as "Likely positive", and a third one (~3000 samples) labeled as "Unknown" or "Unlabeled" or, if you wish, "Likely negative". Therefore, my classes are imbalanced and categorical, but have some sort of incremental order. My variables are binary (presence/absence of a feature). I'm not an expert on machine learning, but so far I've found R randomForest package most suitable for many of the problems I've addressed. Random Forests seem quite robust and useful for different situations. Now, given a test set, in the end what I want to know, for each of the samples, is their probability of being positive. The intuition here is that the "Likely Positive" category should "help" somehow in assessing this probability. How would you guys address this problem? Thanks in advance! Miquel


title: Good resources (online or book) on the mathematical foundations of statistics

text: Before I ask my question, let me give you a bit of background about what I know about statistics so that you have a better sense of the types of resources that I'm looking for. I'm a graduate student in psychology, and as such, I use statistics almost every day. By now I'm familiar with a pretty broad array of techniques, mostly as they are implemented in the general structural equation modeling framework. However, my training has been in the use of these techniques and the interpretation of results -- I don't have much knowledge of the formal mathematical foundations of these techniques. However, increasingly, I've had to read papers from statistics proper. I've found that these papers often assume a working knowledge of mathematical concepts that I don't know much about, such as linear algebra. I have therefore become convinced that if I wish to do more than blindly use the tools that I have been taught, it would be useful for me to learn some of the mathematical basis of statistics. So, I have two related questions:   1. What mathematical techniques would be useful for me to know if I want to brush up on the mathematical foundation of statistics? I've encountered linear algebra pretty often, and I'm sure that learning about probability theory would be useful, but are there any other areas of math that would be useful for me to learn about?   2. What resources (online or in book form) can you recommend to me as someone who wants to know more about the mathematical foundations of statistics?


title: Providing CIs for comparing interaction point estimates across datasets?

text: I have more of a brainstorming type of question for statisticians and I would like to know what your thoughts are. I'll describe my data sets, so you can judge what type of statistical problem they involve. I have two different language experiments: one was done with Spanish speakers and one was done with English speakers. Both experiments have exactly the same design: a 2 by 2 within subjects design, crossing two factors, "grammaticality" (grammatical/ungrammatical word) and "number" (singular/plural word). I am interested in the interaction between these two factors. When I run a my statistical model (a factorial anova using lmer in R) I get a significant interaction in both cases: beta / the interaction coefficient is significantly different from 0, with p <0.05 My goal is to compare the value of beta between the two data sets (since I have two experiments, I have one beta value for the Spanish experiment and one beta value for the English experiment). In order to do that, I was thinking of providing the beta value for each experiment and getting a confidence interval for each beta, with the goal of showing that the two values (and their CIs) are not that different. I am trying to figure out which is the most appropriate way of calculating the 95% CIs. I thought of doing bootstrapping as a way of getting them. So I guess my questions are:   1. Is this way (providing the betas for English and Spanish and the 95% bootstrapped CIs?) a clear way of comparing across experiments? (that is, are there other more appropriate statistical ways of making the comparison across languages?)    2. Is bootstrapping an adequate way of getting the CIs? And if not, how would you get them? Let me know if this was unclear or if you need me to expand and thanks!


title: Quantify mediation effect of one variable on interaction of 3 others

text: I would like some quantitative measure of the degree to which an observed categorical variable predicts mediates the effect of an experimentally controlled factor on a continuous DV. The situation is made more complex by the fact that the effects of the experimental factor are interactions, not main effects. I am not sure that what I need is "mediation analysis" in the technical sense - I may be using that term incorrectly. Here are the details - sorry this is a long story: Participants are classified into two groups (High or Low previous experience with the target domain). They take a pretest consisting of several problems of different types. They then receive training in condition A or condition B. Finally, they take a posttest consisting of more problems of the same types found on pretest. My initial analysis of this data uses a mixed ANOVA with test score as the DV, test section (pretest/posttest) and problem type as within-subject factors, and training condition (A or B) and experience group (High/Low) as between- subject factors. I found an interaction of condition, group, and section, such that in the High group, those in condition A improved more from pretest to posttest than those in condition B, while in the Low group, this trend was reversed. The twist arises because I actually have two more factors. After training, before posttest, participants responded to an open-ended comprehension question. Responses were converted into two binary factors: Correctness (correct/not) and Abstractness (abstract/not). Both factors probably reflect some combination of pre-existing ability & knowledge and the effects of training. I added Correctness and Abstractness as factors in the model described above in order to answer several questions. (1) Do correctness and/or abstractness predict performance? To answer this, I just look at whether they have significant main effects. (2) Do they predict IMPROVEMENT? To answer this, I look at whether they interact significantly with test section. Right so far? Most important though is (3) Do they MEDIATE the effect of training on improvement from pretest to posttest? This is the one I'm not sure how to answer. I suspect I need to look at one or both of these things. Does adding Correctness and Abstractness to the model increase its explanatory power? Does adding them to the model DECREASE the estimated effect of training - remember, not a main effect but an interaction with group and section? I'm not sure how to actually answer these questions quantitatively though. Ideally I would like to have a number which I can compare to a standard and make a yes/no judgment for each of the questions I just posed. Ideas? For reference here are my actual models, in R, both before and after adding the additional factors: > >     ezANOVA( data=D, dv=.(accuracy), wid=.(subjid), within=.(type, section), > between=.(prevexp, training), type=3 ) >     ezANOVA( data=D, dv=.(accuracy), wid=.(subjid), within=.(type, section), > between=.(prevexp, training, soln_corr_bin, soln_abstr_bin), type=3 ) >


title: Calibrated boosted decision trees in R or MATLAB

text: In An Empirical Comparison of Supervised Learning Algorithms (ICML 2006) the authors (Rich Caruana and Alexandru Niculescu-Mizil) evaluated several classification algorithms (SVMs, ANN, KNN, Random Forests, Decision Trees, etc.), and reported that calibrated boosted trees ranked as the best learning algorithm overall across eight different metrics (F-score, ROC Area, average precision, cross-entropy, etc.). I would like to test **calibrated boosted decision trees** in one of my projects, and was wondering if anybody could suggest a good R package or MATLAB library for this. I am relatively new to R, although I have large experience with MATLAB and Python. I have read about R's **gbm** , **tree** , and **rpart** but I am not sure if these packages implement calibrated boosted decision trees or if there are others that implement them. Thanks


title: Validity of estimated survey weights

text: Suppose a sampling mechanism is used to draw a sample ($y_i$ for $i \in \{1,...,n\}$) of size $n$ from a population of size $N$ with probabilities $\pi_j$ where $j \in \{1,...,N\}$. Standard survey sampling tells us that we can estimate the population mean with either the horvitz-thompson estimator: $$ \hat{y}_{HT} = \frac{1}{N}\sum\frac{1}{\pi_i}y_i $$ or the generalized horvitz-thompson $$ \hat{y}_{GHT} = \frac{\sum\frac{1}{\pi_i}y_i}{\sum \frac{1}{\pi_i}} $$ My question is this: If the $\pi_i$ are unknown, but we have unbiased estimates of them ($E(p_i) = \pi_i$), can we use $p_i$ in place of $\pi_i$ in the above estimators and have some sort of consistency guarantee? If so, any references would be appreciated.


title: IV-RE xtoverid testing both IV validity and RE vs FE?

text: I am running an `xtivreg, re` and an `xtoverid` afterwards. My understanding of the help file and what I found online is that `xtoverid` basically compares IV-RE with FE (augmented RE). The test statistic I get is tiny and therefore p-value massive (0.8 to 0.9). My plan was to first test for endogeneity and then afterwards test either the RE-IV or RE versus the respective FE or FE-IV. However, the test seems to combine this into one. So does that imply that both the instruments are valid **as well as** that RE preferable to FE? I find that hard to believe because the FE results differ massively.


title: How to do a power analysis for an unbalanced mixed effects ANOVA?

text: I need suggestions for how to calculate the $n$ required for 80% and 90% power for >30% change from baseline ($T_0$) at $T_1$ or $T_2$, drug vs. placebo, 2:1 ratio; 20% CV in test; factors: subject, treatment, time, and two Gaussian covariates: one of the covariates is the baseline measurement of the endpoint. Please let me know if additional information is required for this calculation. It would also be really helpful if someone could provide a proper statistical description of this study design. I'm just a biologist and not a statistician, so where possible, please use lay terms. Added question: Is there a way to do this with G*Power?


title: Question with MLE

text: I'm having some problems with this question, and was hoping someone here could help. > Let $X_1,\ldots,X_2$ be $n$ determinations of a physical constant $\theta$. > Consider the model > > $X_i = \theta + e_i,$ $~~~~~~i = 1,...n$ > > and assume > > $e_i = \alpha e_{i-1} + \beta e_{i-2} + \epsilon_i$, $~~~~i = 1,\ldots,n$, > $~~~e_0=0, e_{-1}=0$ > > with $\epsilon_i$'s i.i.d standard normal, and $\alpha$ and $\beta$ are > known constants. What is the maximum likelihood estimate of $\theta$? > Carefully justify each step of your derivation.


title: Intuitive explanation of stationarity

text: I was wrestling with stationarity in my head for a while... Is this how you think about it? Any comments or further thoughts will be appreciated. > Stationary process is the one which generates time-series values such that > distribution mean and variance is kept constant. Strictly speaking, this is > known as weak form of stationarity or covariance/mean stationarity. > > Weak form of stationarity is when the time-series has constant mean and > variance throughout the time. > > Let's put it simple, practitioners say that the stationary time-series is > the one with no trend - fluctuates around the constant mean and has constant > variance. > > Covariance between different lags is constant, it doesn't depend on absolute > location in time-series. For example, the covariance between t and t-1 > (first order lag) should always be the same (for the period from 1960-1970 > same as for the period from 1965-1975 or any other period). > > In non-stationary processes there is no long-run mean to which the series > reverts; so we say that non-stationary time series do not mean revert. In > that case, the variance depends on absolute position in time-series and > variance goes to infinity as time goes on. Technically speaking, auto- > correlations to not decay with time, but in small samples they do disappear > - although slowly. > > In stationary processes, shocks are temporary and dissipate (lose energy) > over time. After a while, they do not contribute to the new time-series > values. For example, something which happened log time ago (long enough) > such as World War II, had an impact, but, it the time-series today is the > same as if World War II never happened, we would say that shock lost its > energy or dissipated. Stationarity is especially important as many classical > econometric theories are derived under the assumptions of stationarity. > > A strong form of stationarity is when the distribution of a time-series is > exactly the same trough time. In other words, the distribution of original > time-series is exactly same as lagged time-series (by any number of lags) or > even sub-segments of the time-series. For example, strong form also suggests > that the distribution should be the same even for a sub-segments 1950-1960, > 1960-1970 or even overlapping periods such as 1950-1960 and 1950-1980. This > form of stationarity is called strong because it doesn't assume any > distribution. It only says the probability distribution should be the same. > In the case of weak stationarity, we defined distribution by its mean and > variance. We could do this simplification because implicitly we assumed > normal distribution, and normal distribution is fully defined by its mean > and variance or standard deviation. This is nothing but saying that > probability measure of the sequence (within time-series) is the same as that > for lagged/shifted sequence of values within same time-series.


title: Why the introduction of a random slope effect enlarged the slope's SE?

text: I'm trying to analyze effect of Year on variable logInd for particular group of individuals (I have 3 groups). **The simplest model:**               > fix1 = lm(logInd ~ 0 + Group + Year:Group, data = mydata)     > summary(fix1)          Call:     lm(formula = logInd ~ 0 + Group + Year:Group, data = mydata)          Residuals:         Min      1Q  Median      3Q     Max      -5.5835 -0.3543 -0.0024  0.3944  4.7294           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         Group1       4.6395740  0.0466217  99.515  < 2e-16 ***     Group2       4.8094268  0.0534118  90.044  < 2e-16 ***     Group3       4.5607287  0.0561066  81.287  < 2e-16 ***     Group1:Year -0.0084165  0.0027144  -3.101  0.00195 **      Group2:Year  0.0032369  0.0031098   1.041  0.29802         Group3:Year  0.0006081  0.0032666   0.186  0.85235         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.7926 on 2981 degrees of freedom     Multiple R-squared: 0.9717,     Adjusted R-squared: 0.9716      F-statistic: 1.705e+04 on 6 and 2981 DF,  p-value: < 2.2e-16       We can see the Group1 is significantly declining, the Groups2 and 3 increasing but not significantly so. **Clearly the individual should be random effect, so I introduce random intercept effect for each individual:**               > mix1a = lmer(logInd ~ 0 + Group + Year:Group + (1|Individual), data = mydata)     > summary(mix1a)     Linear mixed model fit by REML      Formula: logInd ~ 0 + Group + Year:Group + (1 | Individual)         Data: mydata        AIC  BIC logLik deviance REMLdev      4727 4775  -2356     4671    4711     Random effects:      Groups     Name        Variance Std.Dev.      Individual (Intercept) 0.39357  0.62735       Residual               0.24532  0.49530      Number of obs: 2987, groups: Individual, 103          Fixed effects:                   Estimate Std. Error t value     Group1       4.6395740  0.1010868   45.90     Group2       4.8094268  0.1158095   41.53     Group3       4.5607287  0.1216522   37.49     Group1:Year -0.0084165  0.0016963   -4.96     Group2:Year  0.0032369  0.0019433    1.67     Group3:Year  0.0006081  0.0020414    0.30          Correlation of Fixed Effects:                 Group1 Group2 Group3 Grp1:Y Grp2:Y     Group2       0.000                                 Group3       0.000  0.000                          Group1:Year -0.252  0.000  0.000                   Group2:Year  0.000 -0.252  0.000  0.000            Group3:Year  0.000  0.000 -0.252  0.000  0.000      It had an expected effect - the SE of slopes (coefficients Group1-3:Year) are now lower and the residual SE is also lower. **The individuals are also different in slope so I also introduced the random slope effect:**               > mix1c = lmer(logInd ~ 0 + Group + Year:Group + (1 + Year|Individual), data = mydata)     > summary(mix1c)     Linear mixed model fit by REML      Formula: logInd ~ 0 + Group + Year:Group + (1 + Year | Individual)         Data: mydata        AIC  BIC logLik deviance REMLdev      2941 3001  -1461     2885    2921     Random effects:      Groups     Name        Variance  Std.Dev. Corr         Individual (Intercept) 0.1054790 0.324775                         Year        0.0017447 0.041769 -0.246       Residual               0.1223920 0.349846             Number of obs: 2987, groups: Individual, 103          Fixed effects:                   Estimate Std. Error t value     Group1       4.6395740  0.0541746   85.64     Group2       4.8094268  0.0620648   77.49     Group3       4.5607287  0.0651960   69.95     Group1:Year -0.0084165  0.0065557   -1.28     Group2:Year  0.0032369  0.0075105    0.43     Group3:Year  0.0006081  0.0078894    0.08          Correlation of Fixed Effects:                 Group1 Group2 Group3 Grp1:Y Grp2:Y     Group2       0.000                                 Group3       0.000  0.000                          Group1:Year -0.285  0.000  0.000                   Group2:Year  0.000 -0.285  0.000  0.000            Group3:Year  0.000  0.000 -0.285  0.000  0.000      ### **But now, contrary to the expectation, the SE of slopes (coefficients Group1-3:Year) are now much higher, even higher than with no random effect at all!** How is this possible? I would expect that the random effect will "eat" the unexplained variability and increase "sureness" of the estimate! However, the residual SE behaves as expected - it is lower than in the random intercept model. Here is the data if needed. ## Edit Now I realized astonishing fact. If I do the linear regression for each individual separately and then run ANOVA on the resultant slopes, **I get exactly the same result as the random slope model!** Would you know why?               indivSlope = c()     for (indiv in 1:103) {         mod1 = lm(logInd ~ Year, data = mydata[mydata$Individual == indiv,])         indivSlope[indiv] = coef(mod1)['Year']     }          indivGroup = unique(mydata[,c("Individual", "Group")])[,"Group"]               anova1 = lm(indivSlope ~ 0 + indivGroup)     summary(anova1)          Call:     lm(formula = indivSlope ~ 0 + indivGroup)          Residuals:           Min        1Q    Median        3Q       Max      -0.176288 -0.016502  0.004692  0.020316  0.153086           Coefficients:                   Estimate Std. Error t value Pr(>|t|)     indivGroup1 -0.0084165  0.0065555  -1.284    0.202     indivGroup2  0.0032369  0.0075103   0.431    0.667     indivGroup3  0.0006081  0.0078892   0.077    0.939          Residual standard error: 0.04248 on 100 degrees of freedom     Multiple R-squared: 0.01807,    Adjusted R-squared: -0.01139      F-statistic: 0.6133 on 3 and 100 DF,  p-value: 0.6079       Here is the data if needed.


title: Unidentified variables in multinomial model

text: I am building a simple multinomial logit (MNL) model. I've discovered that the particular variable that I am testing ("pctLS") causes the model to be "unidentified," but I am not sure how to correct the variable so that it no longer creates this problem. To be specific, the algorithm should converge to a point where the "Final Gradient Norm" is equal to zero. I am using Python Biogeme Here are the outputs from Biogeme: (1) MNL with 8 variables (EXCLUDING the variable "pctLS") (2) MNL with 9 variables (INCLUDING the variable "pctLS") (3) MNL with 1 variable ("pctLS") In the first model (without the problem variable), the algorithm converges appropriately. In the second and third, it doesn't ever find a clean solution. I understand that an unidentified variable can cause the log-likelihood function to have an indefinite solution. How is this? Second, how might I transform a problem variable to avoid such unidentification?


title: Multilabel logistic regression

text: Is there a way to use logistic regression to classify multi-labeled data? By multi-labeled, I mean data that can belong to multiple categories simultaneously. I would like to use this approach to classify some biological data.


title: Probability density function between -1 and 1?

text: I'm currently using Gaussian distribution as a mutation operator for my genetic algorithm. However, I only want to obtain values between -1 and 1. I also don't wish to truncate my Gaussian distribution, which leaves me with lots of 1's and -1's. What type of probability density function can I use to obtain values between -1 and 1, based on a mean value between -1 and 1? Here's an image for the distribution that I'm looking for with mean values of 0, -0.5, and 0.5: ![Distribution Function](http://i.stack.imgur.com/9hQLj.png)


title: How to identify a new pattern in a URL with a machine learning algorithm (Text mining)

text: I am trying to identify new patterns after analyzing a number of URLs. So let's say, I am investigating the hypothetical website Yoohle.com and their URLs have the following structure.   * domain = yoohle.com   * q= search phrase   * lan= language used   * pr= partner_id   * br= browser_id so a sample url will look like this               www.yoohle.com/test_folder/test_page?q=hello+world&lan=en&pr=stackoverflow&br=chrome      If I am investigating the web traffic of this website and seeing abnormal increase month over month, I would like to find out what's causing this. In this example I can just parse out the URL and look at the pr= value since it will tell me if there is a new partnership (maybe stackoverflow is going to be powered by yoohle.com and that drives the increase etc.) The question is, how can I build something robust that can compare 2 (or more) months and tell me exactly what's driving the increase. I want to get something like, "we are seeing an increase and it is driven by the following pattern"               www.yoohle.com/test_folder/test_page%pr=stackoverflow%      The tricky part is, you do not know anything about what the tokens mean unlike this example since I will not know what token stands for partner_id. Another issue is, if we look at token by token, this will be misleading because lan=en will also go up with a new partner assuming the users will still have English as the language. My idea is to analyze the tokens by looking at all the combinations but it is very costly, (4! in this example and probably 10+! for other websites). Also analyzing tokens itself is not going to solve the problem since I still need to analyze the values of the tokens. I tried k-means clustering, apriori algorithm did some research on URL/text mining but could not get what I want. Any ideas about how to approach building an algorithm will be beneficial. Imagine that you are seeing realtime data, so we are talking about analyzing around 100K URLs in a given month.


title: How can I perform an ANOVA from multiple CSV files with R?

text: I (or rather, my girlfriend whom I'm trying to convert to free software) am just getting started with R. She wants to do a two-way ANOVA on her data on mussels and the threads they release in self-defense. Her experiment has four conditions and in each condition two sets of data were collected (number of threads and average length). I am by no means a statistician; I only spent a few weeks learning the very basics of R as a programming language, but I know it is a powerful, industry- grade piece of software. I know I can import data through a command such as               > m10c = read.csv('.../mussels-only-10c.csv')     > m10c        Tank MusselNo ThreadsNo AvgLength     1     1        1         2      1.70     2     1        2         0      0.00     3     1        3         7      1.48     4     2        1        26      2.12     5     2        2         0      0.00     6     2        3        50      0.69     7     3        1         6      1.37     8     3        2        22      1.50     9     3        3         0      0.00     10    4        1         2      1.00     11    4        2         6      1.98     12    4        3        23      1.21      but I don't have a firm enough knowledge of statistics to know what to do with it. A few questions: What is her hypothesis? What are her dependent and independent variables? Is the data independent (did you control for pseudoreplication)? **Hypothesis** There will be significant differences in the thread count and length of mussels in the 10c condition than the 15c condition, and between those stored with snails and those without (but not as drastic). Expecting more threads in the warmer condition. **Independent Variables** Whether there are snails there or not; temperature **Dependent Variables** Number; length of the threads **"Is the data independent?"** She doesn't know enough about it to say whether or not she controlled it - this is an undergraduate research project. She doesn't _think_ it's relevant to an ANOVA.


title: Prerequisite for "Introduction to Statistical Thought"

text: I was browsing through http://www.math.umass.edu/~lavine/Book/book.html and found the book to be really exhaustive in terms of what I need it for but I was concerned about it's prerequisites. I have a background in Calculus (I/II/III) which were fairly applied (as opposed to proof based) and Linear Algebra. Not many other Math courses. I had taken a Probability and Stats course but dropped it midway because I didn't like the exposition. I remember things like average, mode, median etc. but not something to be relied upon. Can I read that book? As far as the book is concerned, the preface states _"This book is intended as an upper level undergraduate or introductory graduate textbook in statistical thinking with a likelihood emphasis for students with a good knowledge of calculus and the ability to think abstractly. "_ I'm not sure what upper level undergraduate means. The reason I'm confused is because most books start off with probability with simple heads and tails type questions. The distributions come in much later but in his book, they come up in the first chapter itself. (Succinct is better, so if this is OK, I'd love to continue reading the book)


title: Scaling in SVM (why and how to , plus references)

text: Hi I know why feature scaling is preferred in SVM, I have two questions: 1-does anyone know of legit articles of books explaining it. I am writing my thesis and I need references. It doesnt have to be for SVM specific, even if it is for Neural Networks or other supervised learning algorithms it will be appreciated. 2-Other than the referencing part does anyone have another source for more explanations on the scaling or rescaling of inputs. Thanks


title: Frequencies calculated over 6 categories

text: I have 60 patients with squint classified according to age group into 6 categories and according to prognosis on their treatment into improved and same. If I want to see if there is a difference among the 6 age groups in regard to prognosis what test statistic should I use? I tried to use chi square but couldn't apply it, first because I have 6 age groups and 2 outcomes which are improved or same and this means it is 6 by 2 table while ordinary chi square is 2 by 2 or 2 by 3 table, and second because I got some expected values less than 5.   * What test statistic should use?    * Should it be Fisher's exact test or goodness of fit test or any other test?    * Would correlation work?    * What is the formula for the right test please?


title: How to interpret odds ratio?

text: Suppose we have a baseline exposure group and 2 other exposure groups for a case control study. Suppose the odds ratio for the first exposure is $1.5$ and the odds ratio for the second exposure is $1.8$. Does this mean that cases are $1.5$ times as likely to have exposure 1 than the controls? Also does this mean that cases are $1.8$ times as likely to have exposure 2 than the controls?


title: "system is exactly singular" in R function BoxCox.ar

text: I'm trying to perform a Box-Cox transformation on some financial data (SPY). The BoxCox.ar function (in the package `TSA`) gives me the following error:               Error in solve.default(res$hessian * length(x)) :        Lapack routine dgesv: system is exactly singular: U[2,2] = 0      Here's the code I'm using:               data = ts(read.csv("http://ichart.finance.yahoo.com/table.csv?s=SPY&d=3&e=21&f=2013&g=d&a=0&b=29&c=1993&ignore=.csv")[,7])     # (Normally I would reverse the data here!)     plot(data)     BoxCox.ar(data)      How can I fix this? I tried smoothing the data with a moving average but the error persisted. Is it possible to perturb the data to avoid the singular system?


title: Maximum expected difference between players when using 4d6 drop lowest

text: Players of a certain TRPG have characters with 6 ability scores, each ability score ranging from 3-18. One method of generating those is by rolling 4d6 drop lowest for each of the scores. That means four six-faced-dice are rolled, and the three highest results are added. What's the expected highest difference in the sum of the ability scores between 2 of those 5 players? The related question here shows how to get the distribution of 4d3 drop lowest, but how do I get from there to an answer of my question above? A good answer would explain the result in a way that a statistics novice can follow.


title: Clustering distributions

text: I have several distributions (10 distributions in the figure below). ![distributions](http://i.stack.imgur.com/Tcz9b.jpg) In fact these are histograms: there are 70 values on the x-axis which are the sizes of some particles in a solution and for each value of x the corresponding value of y is the proportion of particles whose size is around the value of x. I would like to cluster these distributions. Currently I use a hierarchical clustering with the Euclidean distance for example. I am not satisfied by the choice of the distance. I have tried information-theoretic distance such as Kullback-Leibler but there are many zeros in the data and this causes difficulties. Do you have a proposal of an appropriate distance and/or another clustering method ?


title: Best way to perform multiclass SVM

text: I know that the SVM is a binary classifier. I would like to extend it to multi-class SVM. Which is the best, and maybe the easiest, way to perform it? code: in MATLAB                  u=unique(TrainLabel);          N=length(u);          if(N>2)                 itr=1;                 classes=0;                while((classes~=1)&&(itr<=length(u)))                    c1=(TrainLabel==u(itr));                     newClass=double(c1);                  tst = double((TestLabel == itr));                 model = svmtrain(newClass, TrainVec, '-c 1 -g 0.00154');                   [predict_label, accuracy, dec_values] = svmpredict(tst, TestVec, model);                     itr=itr+1;                end             itr=itr-1;         end      How can this be improved?


title: learning R-understanding calculations for specific variable

text: Clarify why the prediction for a female in the example below in related topics (taken from a question by @MsSnowy) do we use the new calculations and not the original lm: `fitted.model <\- lm(spending ~ sex + status + income, data=spending)`, my results were as follows:               Coefficients:                            Estimate  Std. Error t value   Pr(>|t|)         (Intercept)    22.55565   17.19680   1.312   0.1968         sex         **-22.11833**  8.21111  -2.694   0.0101 *       status          0.05223    0.28111   0.186   0.8535         income          4.96198    1.02539   4.839 1.79e-05 ***     verbal         -2.95949    2.17215  -1.362   0.1803       Now, the new model is sex and all other predictors constant in lm model `mydata<-lm(spending ~ sex, data=spending)` was               Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   29.775      5.498   5.415 2.28e-06 ***     **sex        **-25.909****    8.648  -2.996  0.00444 **       **Question: Why isn't the prediction that females spend less than 22.118 (from the 1st lm) than males but rather the new lm coefficient value of less than -25.909?** Someone please clarify because I would think its the first lm for prediction.


title: Gini coefficient and error bounds

text: I have a time series of data with N=14 counts at each time point, and I want to calculate the Gini coefficient and a standard error for this estimate at each time point. Since I have only N=14 counts at each time point I proceeded by calculating the jackknife variance, i.e. $\operatorname{var}(G) = \frac{n-1}{n} \times \sum_{k=1}^n (G(n,k)-\bar{G}(n))^2$ from equation 7 of Tomson Ogwang 'A convenient method of computing the Gini index and its' standard error'. Where $G(n,k)$ is the Gini coefficient of the N values without element $k$ and $\bar{G}(x)$ is the mean of the $G(n,k)$. Direct naive implementation of the above formula for Variance.               calc.Gini.variance <- function(x) {       N <- length(x)       # using jacknifing as suggested by Tomson Ogwang - equation 7       # in the Oxford Bulletin of Economics and Statistics, 62, 1 (2000)       # ((n-1)/n) \times \sum_{k=1}^n (G(n,k)-\bar{G}(n))^2       gini.bar <- Gini(x)            gini.tmp <- vector(mode='numeric', length=N)       for (k in 1:N) {         gini.tmp[k] <- Gini(x[-k])       }       gini.bar <- mean(gini.tmp)       sum((gini.tmp-gini.bar)^2)*(N-1)/N      }      calc.Gini.variance(c(1,2,2,3,4,99))       # [1] 0.1696173      Gini(c(1,2,2,3,4,99))      # [1] 0.7462462      Is this a reasonable approach for a small N? Any other suggestions?


title: Proof for a binomial equation

text: I'm trying to formally write out the problem: Light bulb color problem when I come up with the following equation $\sum_{i=0}^{k-1}{C}_{2k-1}^{2k-1-i}{p}^{2k-1-i}{(1-p)}^{i}=\sum_{i=0}^{k-1}{C}_{2k}^{2k-i}{p}^{2k-i}{(1-p)}^{i}+\frac{1}{2}{C}_{2k}^{k}{p}^{k}{(1-p)}^{k}$ where it is supposed to hold for all $k\geq1$. Let's forget the original light bulb problem for a while. If we just want to prove the above equation, how shall we do it? I tried induction and it seemed difficult and baffling... Is there any idea or perhaps tricks one can use to establish it? What is the general approach for establishing these kind of complex combinatorial equations? (By the way, I know from the light bulb problem that it must be true for $0\leq p\leq1$, but does it also hold for any $p$ beyond that domain?) * * * Edit: I vaguely remember there's a theorem or something that says if the above equation holds for any $0\leq p\leq1$ then it must hold for any $p\in\mathbb{R}$. Does any one know what that theorem is?


title: Sampling an interpolated model with MCMC

text: Is it safe to run a MCMC by interpolating in tabulated data of a model? For background, I have output of a model that involves a set of coupled non-linear differential equations. Calculating models on the fly is slow, so I've tried using a MCMC simulation where the function call is an interpolation in a grid of pre-calculated output, rather than the actual output for a given set of parameters. My question is, basically, has this kind of application been studied before? I've been Googling for scholarly material for a while and looking in some textbooks for clues but I can't find anything that discusses this. My gut instinct is that the results of the MCMC calculation are as reliable as the interpolation scheme but I'm not sure if there are any formal breakdowns in the subtle details of the chains. I have a reasonable background in maths and physics but I'm not a master statistician. I apologize if I've overlooked an obvious answer!


title: Expected distribution of random draws

text: I have a two part question; First Part: I have an urn with 20 balls, 2 of those balls are purple, and I pull out 6 balls at random. I witness 100 realizations of this process. Given the observed frequency at which I drew purple balls, how do I determine if I am really pulling balls out at random? Also given that there are 2 purple balls, I have a hunch that if purple balls are pulled out disproprionately, I expect that both purple balls would be pulled out (i.e. I'm more interested in seeing if 2 purple balls are pulled out disproportionately than I am if 1 purple ball is pulled out more frequently than expected). Second Part: I have an urn with a variable number of balls, a variable number of purple balls within that urn, and a variable number of draws. I witness 100 realizations of this process, and I observe in each realization how many balls there were, how many of those balls were purple, and how many balls I drew from the urn. Same questions; Given the observed frequency at which I drew purple balls, how do I determine if I am really pulling balls out at random? Again I'm more interested to see if higher frequencies of purple balls are disproportionately drawn than I am to see if 1 purple ball being drawn happens more than expected by chance. (I'm open to suggestions for title of question and tags) **Edit:** Srikant suggested I may need to make distributional assumptions about my variables, which I am willing to do. Lets say the number of balls in the urn is uniform between 20 and 30, the number of purple balls is uniform between 0 and 4, and the number of draws is uniform between 6 and 12. See my answer that describes my motivation for asking this question.


title: Combining two confidence intervals/point estimates

text: Suppose one has two independent samples from the same population, and different methods were used on the two samples to derive point estimate and confidence intervals. In trivial cases a sensible person would just pool the two samples and use one method to do the analysis, but let's suppose for the moment that different method has to be used due to limitation of one of the sample such as missing data. These two separate analyses would generate independent, equally valid estimates for the population attribute of interest. Intuitively I think there should be a way to properly combine these two estimates, both in terms of point estimate and confidence interval, resulting in a better estimation procedure. My question is what should be the best way to do it? I can imagine a weighted mean of some sort according to the information/sample size in each sample, but what about the confidence intervals?


title: Adding coefficients to obtain interaction effects - what to do with SEs?

text: I have a multivariate regression, which includes interactions. For example, to get the estimate of the treatment effect for the poorest quintile I need to add the coefficients from the treatment regressor to the coefficient from the interaction variable (which interacts treatment and quintile 1). When adding two coefficients from a regression, how does one obtain standard errors? Is it possible to add the standard errors from the two coefficients? What about the t-stats? Is it possible to add these as well? I'm guessing not but I can't find any guidance on this. Thanks so much in advance for your help!


title: Examples of Bayesian and frequentist approach giving different answers

text: ### Note: I _am_ aware of philosophical differences between Bayesian and frequentist statistics. For example "what is the probability that the coin on the table is heads" doesn't make sense in frequentist statistics, since it has either already landed heads or tails -- there is nothing probabilistic about it. So the question has no answer in frequentist terms. _But such a difference is specifically **not** the kind of difference I'm asking about._ Rather, I would like to know how their predictions for _well-formed_ questions actually _differ_ in the real world, _excluding_ any theoretical/philosophical differences such as the example I mentioned above. So in other words: ### What's an example of a question, answerable in _both_ frequentist _and_ Bayesian statistics, whose answer is different between the two? (e.g. Perhaps one of them answers "1/2" to a particular question, and the other answers "2/3".) Are there any such differences?   * If so, what are some examples?   * If not, then when does it actually ever make a _difference_ whether I use Bayesian or frequentist statistics when solving a particular problem?   Why would I avoid one in favor of the other?


title: Is there any book on statistical distributions?

text: I am a civil engineering student and my research area is traffic engineering. I have encountered different types of statistical distributions e.g. Poisson, binomial, Gumbel, Gaussian, Pearson type III, etc in traffic engineering literature. I understand some of the distributions and what kind of data they fit. But there are many distributions (e.g. Pearson distributions, Gumbel distribution, etc.) which are completely alien to me: I can't find a thorough explanation online or in any book which I could understand and apply in my research. Is there any book or online resource covering all statistical distributions with their applications?


title: Getting Bayes Prior

text: We want to predict how many points per game LeBron James is going to get per game. Assume there is some underlying theta that does not change game to game, that predicts how many points he will get. For previous information, we have NBA league averages, and we have Lebron James points per game for previous years and games. How do you choose a prior here? It would be easy to get say a league average of points per game on all starters who play roughly the same minutes as Lebron. And it is easy to get Lebron's point per games for previous games. The main question is, if I am using a beta distribution to describe the prior, how do I decide what level of sureness to assign the prior? How do we combine Lebron's previous years of points, the league average, and a beta distribution into an accurate description of a prior for Lebron and how sure we are of it (without making arbitrary decisions). Thanks for any thoughts!!


title: How do I sample without replacement using a sampling-with-replacement function?

text: I vaguely recall from grad school that the following is a valid approach to do a weighted sampling without replacement:   1. Start with an initially empty "sampled set".   2. Draw a (single) weighted sample with replacement with whatever method you have.    3. Check whether you have already picked it.   4. If you did, ignore it and move to the next sample. If you didn't, add it to your sampled set.    5. Repeat 2-4 until the sampled set is of the desired size. Will this give me a valid weighted sample? (Context: for very large samples, R's `sample(x,n,replace=FALSE,p)` takes a really long time, much longer than implementing the strategy above.) **UPDATE 2012-01-05** Thank you everyone for your comments and responses! Regarding code, here's an example that faithfully reproduces what I'm trying to do (although I'm actually calling R from Python with RPy, but it seems the runtime characteristics are unchanged). First, define the following function, which assumes a sufficiently long sequence of with-replacement samples as input:               get_rejection_sample <- function(replace_sample, n) {         top = length(replace_sample)         bottom = 1         mid = (top+bottom)/2         u = unique(replace_sample[1:mid])         while(length(u) != n) {             if (length(u) < n) {                 bottom = mid             }             else {                 top = mid             }             mid = (top+bottom)/2             u = unique(replace_sample[1:mid])             if (mid == top || mid == bottom) {                 break             }         }         u     }      It simply does a binary search for the number of non-unique samples needed to get enough unique samples. Now, in the interpreter:               > x = 1:10**7     > w = rexp(length(x))     > system.time(y <- sample(x, 500, FALSE, w))        user  system elapsed       12.820   0.048  12.902      > system.time(y <- get_rejection_sample(sample(x, 1000, TRUE, w), 500))        user  system elapsed        0.311   0.066   0.378      Warning message:     In sample(x, 1000, TRUE, w) :       Walker's alias method used: results are different from R < 2.2.0      As you can see, even for just 500 elements the difference is dramatic. I'm actually trying to sample $10^6$ indices from a list of about $1.2 \times 10^7$, which takes hours using `replace=FALSE`.


title: At which extent can standardized moment be trusted in comparison to plotting data?

text: ## Motivation I am interested on using R to do data analysis on a considerable amount of data. Having to plot each time I want to observe if it fits any particular distribution by eyeball (which can be misleading), and recall in memory how the shape looks like (or compare, say, on a book). As a programmer, I would like to be able to have the program instead looking this to me. Standardized moment was something I studied in an undergraduate course and seems that Kurtosis is of particular interest of the community according to my search. ## Questions   1. Wikipedia -- Kurtosis seems to have a considerable value on them that is very precise on saying the distribution simple based on the number. How far can I trust this instead of eyeball the plotted data? If this is not a yes/no answer, how can I claim that after a given value it is acceptable to apply, say for example Linear Regression assuming Kurtosis indicated a normal distribution (given the other constraints are correct)?   2. The same as above, but for moment 3, Skewness. I would also appreciate if someone would point how to understand if it is left or right skewed or symmetric based on the values of 'the method'.   3. When I mention 'the method' above I noticed wikipedia suggest few of them. So the third question is, is there currently a chosen one or depends on the context? If depends on the context, how can I go about analyzing? I would appreciate if someone could point an R function/package that does this as well. Thank you.


title: Is this approach using the central-limit theorem applicable for reasoning about open source project data?

text: ## Definitions **Open Source Project** : > In short (and roughly described, only to the purpose of clarification) open > source projects allow me to have access to the code of a certain 'program'. ## Motivation and Background (so I can understand your answer) I have a undergoing majors in CS, expect a very high level course in statistics and self-learning using books. Do also consider I am trying to be practical here as I will do most of the calculations in R (that doesn't means your explanation needs to be in R), but its just a pointer that I am concerned on understand what I am doing without going to every single detail. ## Question   * It might have flaws on the theory along the path, much I salvaged from books on my own. Please feel free to point my mistakes. According to **Blalock Jr, Social Statistics, Mc Gran Hill** we have that **The Central Limit Theorem** is as follows: > If repeated random samples of size N are drawn from any population (of what- > ever form) having a mean $\mu$ and a variance $\sigma^2$, then as N becomes > large, the sampling distribution of sample means approaches normality, with > mean $\mu$ and variance $\sigma^2/N$. Now, one of my interests is analyzing open source data using some statistics rigor. I wanted to use the **Central Limit Theorem** to reason about data, as one of the most trouble some issues I have been dealing with is that they do not follow a normal distribution. My background course in statistics used as textbook **Probability and Statistics by Paul Meyer**. There and probably anywhere else we have the definition of random variables, however with a much great emphasis. The first part of the approach I wanted to use starts with a doubt at this point: The first book approaches binomial distributions, t student distribution, chi square, ANOVA, and others after talking about hypothesis testing. Thus, as defined by the book my concern is:   * To define and follow the basic rules of an experiment and given the assumptions of one of these distributions.   * Consider one of them to be my hypothesis and the remaining I must assume to be true (which constitutes my model).  The second book however is more concerned on observing each distribution as a random variable and does not follow the reasoning of choosing a distribution to verify a hypothesis when introducing the distributions. In fact, hypothesis testing is introduced much later.   * **Question 1** : Are these two different approaches to analyze data or am I being mislead by the order of subjects being taught on different textbooks? * * * Let me assume at this point despite question 1 that I decided to define the experiment first as wanting to test my hypothesis on an open source data. **Pang, on Introduction to Data Mining** book suggest that there are two types of data: >   * **Observational Data** : Data was already available and we are trying to > understand something out of it. >   * **Experimental Data** : An experiment is settled so we can control each > of the variables during the experiment, and thus based on my description of > the book above of choosing one distribution assumption as hypothesis and the > remaining as true seems appropriated. >   * **Question 2** : Using hypothesis testing, given this definition would be inappropriate to approach open source projects given that they fall under observational data or is this tolerable?  * * * Assuming that anything at this point poses a serious violation, I moved on the approach to believe that if observing, for example, the number of open defects or closed defects on a given open source project doesn't follow a normal distribution, I could obtain one out of the distribution of the mean of the samples and use it instead. After this I would try to see which distribution seems more appropriated based on its assumptions from the parametric 'family', since the normal distribution constraint which is the most troublesome one would be gone due to the **Central Limit Theorem**. The first book also suggests that N>100 is enough for a given sample to be considered a large N, so for a open source project this is requirement is easily match.   * **Question 3** : When we go about using the distribution of the means, can we actually use it for creating hypothesis or I would need to have the 'original data' to follow the normal distribution? (In this example open defects and closed defects).   * **Question 4** : Which would be better: Using the whole open source project to this point as one single huge sample or extract a smaller (which size?) sample out of it? * * * Question 3 is actually of the most particular of interest of mine since I am unable to use any parametric distribution for hypothesis testing with the 'raw data'.   * **Question 5** : How much do I lose on being able to infer only from the distribution of the samples instead of, say, the real distribution of open defects and closed defects?   * **Question 6** : Would be better to attempt data analysis from the open source code using non-parametric methods? Two possibilities that I have salvaged from books seemed possible here:     * Attempt to convert a certain observation from the data to Binomial distribution (leaving to a binary mapping) so that the normal distribution constraint does not apply.      * Going for ordinal scales using ranking, in exchange losing statistical inference power (How much, in this case does usually is?).  Thank you.


title: Marginalizing a conditional distribution

text: Suppose I have a conditional discrete distribution $P_{M|A}(m,a)$. If I take sum over $m$, i.e. $\sum_{m}P_{M|A}(m,a)$, do I get $f_A(a)$, where $f_A(a)$ is the pmf of the random variable A?


title: Modeling seasonality in sales time series

text: Suppose I have 5 years of weekly sales data for a particular product. I also have other variables such as weekly unemployment rate, weekly coupon rates, and percentages of marketing amount spent on internet, TV, and mail advertising. Suppose I want to determine whether season has an effect on sales. Would it be best to treat season as a dummy variable where season takes values: Winter, Fall, Spring Summer, or would it be better to perform spectral analysis on the sales time series and determine seasonality that way? The current setup I have is: $sales_{t} = week_{t}+ unemploy_{t} + coupon_{t} + internetprop_{t} + TVprop_{t} + mailprop_{t} + season_{t} + w_{t}$ where $w_{t}$ is white noise. I run a linear regression first and depending on the ACF and PACF charts, I choose a model for the error term. Is the correct approach? I am unsure of how to model seasonality.


title: Is a rank means stat appropriate?

text: I have attached an image of a subset of data consisting of snow measurements (D1-D6) taken at unique locations (Site1-Site25) over the course of several years. Each measurement was taken at the same location each time, once a month for 2 years. I'm trying to convert these measurements of snow depth into a single variable for each site to represent which sites consistently have deeper snow. At this point, I don't care about differences in means or anything like that, I just would like a value for each site representing it's "relative snow depth" compared to all other sites. So far I've essentially used a Kruskal-Wallis H-test to determine a mean rank for each site (without then conducting the test) and I'm wondering if that appropriately represents the relative snow depth I'm trying to calculate. If not, does anyone know of a test that may be more appropriate for this data? Thanks!! ![Snow data frame](http://i.stack.imgur.com/WTweU.png)


title: How to explain smoothing functions in the logistic regression model

text: I fitted a logistic regression model with some smoothing functions, and the software made beautiful plots for them. Here is one example: ![enter image description here](http://i.stack.imgur.com/caGUb.png) My main concern is that there is no reference level in the continuous variable, so converting the estimates by an exponential function to odds ratio seems not quite right for smoothing functions. Hence, I am wondering how to explain smoothing functions appropriately? Here is a part of output of the above smoothing function:               age estimate  95% CI     18  0.802419 (0.708936, 0.90305)     19  0.744268 (0.666983, 0.823497)     20  0.688122 (0.625178, 0.755262)     21  0.634746 (0.57693, 0.69252)     22  0.584887 (0.530877, 0.638431)     23  0.538855 (0.485444, 0.589725)     24  0.496467 (0.44487, 0.550205)     25  0.457521 (0.405356, 0.511471)     26  0.421809 (0.370094, 0.47534)     27  0.389193 (0.340622, 0.441149)      Any suggestion will be appreciated!


title: Neuro fuzzy system for classification

text: I am applying Neuro fuzzy system for classifiying objects. I found only a single paper **"Rough-Neuro-Fuzzy System with MICOG Defuzzification"** and it is very hard to follow how the output which is the class is defuzzified and the interpretation of the defuzzified class. For example: If object{color} is Green And object{size} is small And object{shape} is Square Then Object is GreenBox. Hence, the Situations are the class labels {GreenBox, RedCup, Table}. The application is that these represent situations where the shopper is at a particular counter. That is the shopper is shopping for boxes or cups or furniture. I have questions regarding the defuzzification step of neuro fuzzy system used in classification. Neurons in layer 1 provide input values to the network. Neurons in layer 2 represent fuzzy sets defined on the input variables: each neuron receives the input value and computes a membership value through a Gaussian function. Each neuron in layer 3 corresponds to a rule; it is connected to the neurons in layer 2 which implement the fuzzy sets and computes the activation strenght as product of the membership values. Layer 4 is the class variable. Say there are 3 classes C1: [0 1 0]; C2: [0 0 1]; C3: [0 1 1]; and there are 5 input nodes. According to fuzzy set theory, after the variables are fuzzified there should be a module for defuzzification. So, how is the class label defuzzified and then we can say with a confidence that 85% probability that the shopper is shopping for furniture if the input patterns are labelled under Class 3. Shall be greateful for insights. Thank you