title: What statistical tests can be used on a distribution of distributions?

text: I have two experimental conditions, each with 20 reps (40 experiments total). In all 40 experiments a non-fixed number of sounds were recorded, and I've measured each of their durations. I would now like to test the null hypothesis that, on average, the sound durations in condition A are the same as in condition B. However, I'm a little unsure how to do this in a way that captures all of the data. ***I should note that I have the same data format for a number of pairs of conditions A and B, and sometimes the data are paired and sometimes they are not paired. I'm struggling with both data types. My current ideas are:   1. Pool all of the measurements from each condition into two giant distributions and do a t-test (I'm going to ignore the question of parametric vs. non-parametric testing for simplicity's sake). This seems reasonable, but in the paired case it ignores the pattern between the two conditions across individuals.   2. Reduce each distribution to a measure of central tendency and then do a t-test on the distributions of MoCTs. This seems a bit fishy, as it doesn't really take into account most of the data in the raw data distributions, but it does work equally well in both the paired and the unpaired data conditions.   3. In the paired data condition, do individual t-tests between each pair of paired data, then do a t-test to see if the resulting distribution of t-stats is different from 0 (which would indicate that in some pairs the sound durations were longer in condition A and in some cases they were longer in condition B, but across the population of pairs there isn't a strong pattern in one direction or the other). I feel that this scheme does the best at taking all of the data into account, but I don't know how to do something equivalent with the unpaired data.


title: How can I estimate parameters of a convolution of an exponential and gamma?

text: Ideally, I would input a one-dimensional array of data, and output the estimates for the three parameters. I'm not very familiar with any statistics software but I have MATLAB (and all of the free software available, I suppose).


title: Plotting ratios from x:1 to 1:x

text: I'm comparing the ratio of two measured variables as either 1:x or x:1, depending on which is higher. Absolute difference isn't so important, rather the proportion, that's why I'm trying this method. On the plot's y axis, I don't want any decimal numbers below 1 (being 1:1). Instead, I would like marks on y at say 3:1, 2:1, 1:1, 1:2, 1:3, all with equal spacing. some raw data as an example               Date        Control       Treatment     24/10/2013  38.16702859   50.98325719     5/11/2013   40.40145245   39.76689213     13/11/2013  25.79553081   54.50968963      proportions are to be graphed as               1 : 1.335793198     1.015957001 : 1     1 : 2.113144716      So a line joining the three should start above 1:1, then dips just below 1:1, then back above. Plotting these in Excel includes the area between -1 and <1 where values cannot exist, so I need to fold this space out of the graph and have 1:1 as the origin (so to speak), which it seems Excel can't do so I'm after other software. I imagine R might be suitable here. Firstly, is there a name for this particular type of graph so I can find more info? Secondly, is there some existing software that can work with such number pairs and easily plot data in this way?


title: Generating (forcing) confidence percentages in RapidMiner

text: I have a dependent variable (my 'label' in RapidMiner terms), that is a binary classification expressed as 'WIN' or 'NOTWIN'. I know 'NOTWIN' is exhibited in about 90% of all observations. When I try to run a K Nearest Neighbors approach, the generated results do not show a confidence percentage. I think this has got to do with the fact that 'NOTWIN' is overwhelmingly the correct classification. Nevertheless, what is important to me is the confidence percentage assigned to 'NOTWIN' (however large it may be), so as to gain an understanding of observations that fall into the 'WIN' classification. How do I 'force' RapidMiner to show me a confidence percentage given this scenario?


title: terminology question: Is there a generic name for the set of probabilities corresponding to quantiles

text: Given a set $P_i$ of probabilities and a Cumulative Distribution Function of $X$ you can find the set $X_i$ of values of $X$ corresponding to $P_i$. The $X_i$ are the "quantiles". Is there an equivalent generic term for referring to the $P_i$ corresponding to a given set of $X_i$? You could refer to them as the CDF but in standard use that refers to the $P_i$ corresponding to _all_ the $X_i$ rather than some designated subset of interest. In the absence of a pre-existing name I'll probably use some awful phrase like "Corresponding Cumulative Probabilities" (abbreviated as CCP and pronounced "soviets").


title: Is there a verifiable way to generate discrete random variables on the basis of an external event?

text: Suppose that we want to generate a draw from the following distribution: $P(X=0) = 0.5$ $P(X=1) = 0.5$ There are two constraints though: (a) The draw has to be on the basis of an external event. (b) Related to (a), the draw must be verifiable by a third party. In other words, a third party should be able to verify that my draw was in fact $X = 0$ (say). Qn 1: Can such a system be devised and if so how? Qn 2: Can the system be extended to discrete variables with more than 2 possible outcomes (like the roll of a dice)? Qn 3: Similarly, can it be extended to continuous variables (e.g., the normal)?


title: Metric for ordinal, non-normal distribution in Markov models

text: I am working with Markov models with ordinal, non-normal distributions of probabilities. Ultimately, I would like to create a metric for determining the probability of a specific path occurring. At any given point along the path, there is a range of possible moves with associated probabilities. However, these are not normal distributions which makes it difficult to discuss their statistical likelihood. I would like to create some type of metric that would measure the amount of "deviation" from more probable paths. Because the data varies so drastically, it is difficult to get a good idea of how likely each event was. A sampling distribution doesn't say much because variation from one path to the next path could be a difference of up to 60%. Also, at each state there might be one possible move or 30. This high degree of variation between states makes it difficult to measure.


title: Collection of continuous variables with >70% more zeros

text: This is dataset that is going to be data mined for factors that affect an output that of interest A large Part of my dataset (150 of 300 potential inputs) has a heavy skew of Zero values. usually this is >90% and always >70% zeros (true zeros) This has risen as the input variables are usually OR options. -Say i have 10 input variables of the above in a wider dataset. of these 10 any one data line will have approximately 1-4 of these 10 populated with non zeros. In real world terms if you use input variable 1 you may not need variable 2 to 10 the above statement isnt absolute. you may use a some of input 1 and some of 2 and none of 3 to 10. i was tempted to transform these into binary catagorical 0 = No 1 = Yes but i would loose the scale on the non zeros. -If input variable 1 IS NOT 0 the values may range between 2 orders of magnitude


title: What do you call a graph where you color a grid based on some measure?

text: I mean something like this chart showing the frequencies of a user's public contributions on GitHub. ![enter image description here](http://i.stack.imgur.com/X9SRH.png) I really like this kind of chart, as I think that they display the intended data in a really nice and easy to see way, but I have no clue what they are called, so I can't look for a library to help me generate them. Can someone tell me the name of this type of graph / chart?


title: about the definition of bayesian network

text: In this PDF http://people.csail.mit.edu/yks/documents/classes/mlbook/pdf/chapter2.pdf page 5 says: Given a set of functions $f(x_i,pa(x_i))$ non-negative and sum to 1, we define a joint probability as follow: $$p(x_1,\dots, x_n)=\prod_i f(x_i,pa(x_i))$$ [...] what is the relationship between the conditional probability $p(x_i|pa(x_i))$ and $f(x_i,pa(x_i))$, a function with has the properties of a conditional function but it is arbitrary? As will see, this function is in the fact one and same. That is, under the definition of joint probability, the function $f(x_i,pa(x_i))$ is $p(x_i|pa(x_i))$. How can it go from $f(x_i,pa(x_i))$ to $p(x_i|pa(x_i))$ and state that it is the only possible function that can be assigned to $f(x_i,pa(x_i))$?


title: HPD interval for the mean

text: Suppose we have iid observation with the following model $ Y_t \sim \mathcal{N}(\mu,1/\mu) , t=1,2,..T$ The question is " Assuming a flat prior on $(0 ,\infty )$ find a 95% HPD interval for $\mu"$ The aditional information that the problem gives is: $T=6; \sum{y_i}^{2}=76;\sum{y_i}=18$ * * * this is my first try, $p(\mu)=\mu^{-1}1_{y>0}$ then $ \begin{eqnarray*} p(\mu\vert Y) & \propto & p(Y \vert \mu) p(\mu) & \propto \mu^{T}\exp[\sum{(-\mu/2){}(y_i-\mu)}^{2} ] \mu^{-1}1_{y>0} \end{eqnarray*} $ $Pr(\mu ϵ C│y)=\int_{a}^{b}p(\mu\vert Y)d\mu=\int_{a}^{b} \mu^{T-1}\exp[\sum{(-\mu/2){}(y_i-\mu)}^{2} ] 1_{y>0}d\mu=0.95$ * * * I think this problem doesn´t have an analytical solution, so the solution is numeric. The question is which algorithm ,step by step, can I use to solve it . * * * Update $\sum{(y_i-\mu)}^{2}= \sum{y_i}^{2}-2\mu\sum{y_i}+T\mu^{2}=76-36\mu+6\mu^{2}$ $\ (-\mu/2) \sum{(y_i-\mu)}^{2}= -38\mu+18\mu^{2}-3\mu^{3}$ $Pr(\mu ϵ C│y)=\int_{a}^{b}p(\mu\vert Y)d\mu=\int_{a}^{b} \mu^{5}\exp[-38\mu+18\mu^{2}-3\mu^{3} ] 1_{y>0}d\mu=0.95$                   %%Matlab Code         a=0.0001;b=0.0001;step=0.0001;         x=[a:step:1];         fx= x.^(5).*exp(-38.*x+18.*(x.^2)-3.*(x.^3)); plot(x,fx)      ![Density](http://i.stack.imgur.com/RqE7r.jpg) How can I find (a,b) that give me a 0.95 credible region with the minimum distiance between a and b


title: Create grid and color cells with average values of scatterplot using ggplot2

text: Given a numeric dataset $\\{(x_i, y_i, z_i)\\}_{i = 1}^N$, one can create a scatterplot by drawing a point $P_i = (x_i, y_i)$ for each $i = 1, \dots N$ and color each point with an intensity depending on the value of $z_i$. Assume I have created this scatterplot. What I would like to do is to divide the plane with a grid (rectangular, hexagonal, triangular, ... doesn't matter) and color each cell of the grid with the average intensity of all the points that fall within the cell. Additionally, suppose I have created two such plots $Plot_A$ and $Plot_B$ (with identical boundaries) for two different datasets $A$ and $B$. Let $c_i^k$ be the $i$-th cell of plot $k$. I want to create a third plot such that $c_i^3 = c_i^1 * c_i^2$ for every $i$. Thank you.


title: Does an exponential model fit my data?

text: I am measuring accumulation of a fluorescent-tagged protein at a particular location within a cell over time. In previous experiments that I have performed, I see a standard exponential distribution where the fluorescence intensity reaches a plateau, however in my current experiment, I see a distribution as shown below: ![enter image description here](http://i.stack.imgur.com/dJFX7.png) What is the best model to use for this data? Should I use two separate exponential models, one for the increase in intensity up to the peak and one for the decay phase, or is there another statistical model for this type of distribution. Thanks!


title: How to determine the demographics of consumers of a product, given demographic information and consumption rates

text: Okay, so I'm going to admit right away that I'm very new to statistical analysis, so please go easy on me. I have a 50+ different groups of people and for each group, I have the number of people using a product, and the total number of people in the sample. So, for group number 1 of 50+, 353 out of 26719 use this product. Then, for each of these groups I have demographic information. So, for example, 4078 of the same 26719 people are white males between the ages of 35 and 55, or whatever. And I have this for 50+ groups and many demographic categories. How do I use this information to determine who it is that is using our product? Who isn't using our product? All that jazz.


title: How can I counteract the effect of a degenerate classifier in an OVA Model?

text: Suppose I build a OVA classification model for classification with more than 2 possible classes (a model of sub-models, where each submodel predicts the probability of a data point belonging in a single class). Occasionally, one or more of the submodels might not perform very well. For instance, the model for predicting membership in class 1 might be something silly like $P(y=1|x) = 1$. Since this is the highest possible value, the OVA model is going to assign everything to class one, regardless of the quality of the other sub-models. When dealing with many classes, this seems to be a fairly frequent event (i.e. at least one of the classes is hard to predict, so we get a degenerate model). I can see a number of ad-hoc ways of correcting this. The simplest is to generate a confidence factor $\frac{tp}{tp+fp}$, which is essentially the probability that a positive prediction made by the sub-classifier is correct. If I multiply P(y=1|x) by the confidence factor, then a degenerate classifier's probability estimate is re-weighted to just be the percentage of the time it would be correct, which seems reasonable. Likewise, a perfect classifier is not reweighted at all. I imagine other people have encountered this problem, but don't know what the name of this problem is, or what a standard approach to fixing it might look like. Any insights are much appreciated.


title: Usage of rma.glmm() for random effect single group study

text: Below are first few rows of my data that need to be analysed. I would need to understand the relation between number of years and survival rate. My statistical teacher advised on using rma.glmm() function in R language to do poission regression. However I find glmm() is being used only for two group of variables and not for single group. Can anyone advise the correct formula for glmm() for my case where the objective is to obtain correlation between number of years and survival rate.               Study        num_years  sample_size     survived    failed     Mohl, 1988      2         33               28         5     Aboush, 2001    2         57               49         8     Probster , 1990 2         48               47         1     Clyde , 1988    2.5       122             109        13     Wahad, 2004     3         21               19         2     Rashid, 2003    3         70               66         4     Isidor, 1992    4         20               20         0


In the example from the web site I was trying out this problem in **page 107** about Goodness of fit test for a normal distribution.Question is about analysis of fat content of hamburgers. I understand that we need to have the two intervals **-$\infty$ < x< 26** and **$40 < x<$$\infty$** to find the class probabilities of intervals 26< x <28 and 38< x<40.   But why should the **expected frequencies of the two intervals -$\infty$ < x< 26** and **$40 < x<$$\infty$** be also **calculated and used for chi Squared value**. These two categories are **not** there in the given data set. Is it necessary to break our intervals such that they start from -infinity and end from infinity as this is a normal distribution?


title: Hamiltonian Monte Carlo: why is reparameterizing needed?

text: In the Stan user's manual (Version 2.0.1, page 157), it says > A hierarchical model such as the above will suffer from the same kind of > inefficiencies... [for a Hamiltonian Monte Carlo method] because the values > of `beta`, `mu_beta` and `sigma_beta` are highly correlated in the posterior They then recommend reparameterizing the model to avoid this issue. I'm not sure I have good intuition about why correlations would be a problem. I thought that the whole point of Hamiltonian Monte Carlo was that it was invariant to rotation. So it shouldn't matter whether the posterior is aligned with the axes or if correlations put the major axis at an angle. The best guess I can come up with is that, if the correlations get to be extreme enough, the posterior could end up very long and narrow, and that something about the difference in scales causes problems. But I'm not sure that this is a reasonable interpretation.


title: How to Reduce Error Term

text: My question is "What could you do if you wanted to reduce the error term (e)? I know the error term is basically the distance between the line and the point but I don't know how you would reduce it. The problem is: > The Postal Service is attempting to reduce the number of complaints made by > the public against its workers. In order to facilitate this task, a staff > analyst for the Service regresses the number of complaints lodged against an > employee last year (Y) on the hourly wage of the employee for the year (X). > She obtains the following results: Y = 10.2 - 1.6X + e r2 = .73 Syx = 1.1 sb > = .87 n = 348


title: Bootstrapping in R using the boot {boot} and Boot {car}

text: I'm trying my hand at resampling techniques with a dataset I have, and I think either I'm missing a conceptual point with bootstrapping, or I'm doing something incorrectly in `R`. Basically, I'm trying to use it in a correlation/regression framework, and I'm able to get the original coefficients, the bootstrap bias, and the bootstrap coefficients but I can't find a way to have `R` easily display the bootstrap model $R^2$ (when I'm working with several predictors), the Pearson $r$, or the $p$-values for individual regression coefficients. (I'm using the `Boot` function in the `car` package). A secondary question...the more general function `boot` in the `boot` package requires defining a function to use as an argument. The function must include an argument for the original data set, and a second argument which is a set of indices, frequencies, or weights for the bootstrap sample. I'm a little confused by this. What conceptually are these indices I am specifying, and how do I specify them syntactically within my function?


title: Generalized linear models dispersion analysis

text: I was wondering, in relation to model selection for some count data I am looking at using poisson GLMs, would a high residual deviance (measure of dispersion) warrant switching to a quasi-poisson model before or after dropping unimportant explanatory variables?


title: Power analysis in a linear model - pooled vs individual samples

text: I'm interested in how many samples are required to see a relationship between allele frequencies (proportions of alleles) and environmental gradients, such as temperatures and water quality. Environmental gradients are either continuous or categorised into 5-9 levels. I understand that we can use pwr.2p.test(h = , n = , sig.level =, power = ) to estimate required sample size when comparing two proportions (e.g. between control and a treatment). But I want to see a linear relationship. I found that the following R code can be used to estimate sample size in a linear model:               pwr.f2.test(u =, v = , f2 = , sig.level = , power = )      where u and v are the numerator and denominator degrees of freedom. We use f2 as the effect size measure. The code above requires numerator and denominator degrees of freedom. I understand that numerator degrees of freedom is k - 1, denominator degrees of freedom is N - k, where k is the number of treatments and N is the total sample size. My allele frequency data is pooled data. I have a total of 18 populations (i.e. N=18) which consist of 40~50 individuals each. And let's say k (=the number of environmental levels) is 5. Cohen suggests f2 values of 0.02, 0.15, and 0.35 represent small, medium and large effect sizes. I picked small effect size. I input the information in the code: > pwr.f2.test(u = 4, v = 13, f2 = 0.02, sig.level = .05, power = )                Multiple regression power calculation                     u = 4               v = 13              f2 = 0.02       sig.level = 0.05           power = 0.06317734      a power is very very low. If I use the number of individual samples (N=~1000: this is the number of samples that I used. i.e. not pooled.), I get a power of >0.95. Maybe this is not a right way to estimate sample size, is it? Please let me know what you think. Is there other way to estimate required sample size when a linear model is used?


title: Help with R Zeroinfl model

text: I am carrying out a zero-inflated negative binomial GLM on some insect count data in R. My problem is how to get R to read my species data as one stacked column so as to preserve the zero inflation. If I subtotal and import it into R as a single row titled Abundance, I loose the zeros and the model doesn't work. Already, I have tried to: stack the data myself (there are 80 columns * 47 rows) so with 3760 rows after stacking manually you can imagine how slow R gets when using the pscl zeroinfl() command (It takes 20mins on my computer!, It still worked) The next problem concerns a spatial correlation. Certain samplers sampled from the same medium so as to violate independence. Can I just put medium in as a factor in the model? platypezid


title: Can a negative binomial distribution be used to model a continuous distribution?

text: I have a data set which is a set of continuous distances from some origin. I originally modeled this as a negative binomial distribution by rounding the data and using it as an input in the Matlab function, nbinfit:               nbReg = nbinfit(round(data));      Using the output from nbinfit, I created the following CDF: ![enter image description here](http://i.stack.imgur.com/FFYjs.png) Here are the corresponding QQ Plots ![enter image description here](http://i.stack.imgur.com/Qlaz4.png) **Questions**   1. Is it valid to use a negative binomial distribution to model a continuous distribution? Is there any documentation in the literature where negative binomial was used to model continuous data?    2. It is possible to have a similarly distributed continuous distribution   (you can remove the x == round(x) condition in nbinpdf to get a continuous distribution)?   3. The CDF fits very well, but is there a "better" distribution you would recommend?**


title: Mixed effects vs linear regression: accuracy of beta estimation

text: Let's say you have a response variable and an independent variable. Your data is measured across several levels of a categorical independent variable. One approach in analysing these data would be to use linear regression to estimate a slope at each level of the categorical independent variable. This is the approach I've used here, using `sleepstudy` dataset from the `R` `lme4` package (I've stored the betas from each model in `lmBetas`):               library(lme4); library(plyr); library(ggplot2)     lmBetas <- daply(sleepstudy, .(Subject), function(x) coef(lm(Reaction ~ Days, data=x))["Days"])      Another approach in analysing these data would be to use a mixed effects model to estimate slopes for each level of the categorical independent variable, which in this case is `Subject`. This is the approach I've taken here (I've stored the betas from the model in `lmerBetas`):               lmerBetas <- coef(lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy))$Subject[,"Days"]      I have learned that a single mixed effects model, as implemented through the `lmer` function in R, is more accurate at estimating slopes than a multiple linear regression model applied to multilevel data. This can be demonstrated with this plot of betas from the above models.               betas <- data.frame(method.betas = c(lmerBetas, lmBetas))     betas$method <- c(rep("lmer", 18), rep("lm", 18))          ggplot(betas, aes(method.betas)) +        geom_histogram() +       facet_grid(method ~ .)      The top histogram shows betas estimated using linear regression, and the bottom histogram shows betas estimated using mixed effects. You can see betas estimated using linear regression are more widely spread than those estimated through the mixed effects model. ![enter image description here](http://i.stack.imgur.com/YtB22.jpg) So finally, my questions:   1. Is a mixed effects model's higher accuracy in betas estimation connected with the fact that it models intercepts and slopes for each level of the categorical independent variable under a joint probability distribution?   2. Generally speaking, why is a mixed effects model more accurate in its betas estimation?


title: Is there a well-defined class of ensemble methods?

text: Ensemble methodology's main aim is to somehow aggregate or summarize estimates from multiple models. In some cases this is aggregating different bootstrap estimates or Monte Carlo estimates, but others involve combining the estimates from starkly different models. When you look up "ensemble methods", especially in ESL, they mention boosting, bagging, bumping, random forests, and the EM algorithm. Question is: is this all of them? Or more generally, are all so- called ensemble methods somehow contained within these classes? Or is there a more general term that I'm missing? The list so far:   * Boosting   * Bagging   * Bumping   * Random forests   * EM algorithm Can you add to this?


title: Plot Pareto tails in QQ-plot for log-normal distributions

text: I'm working on samples that I'm trying to fit into log-normal distributions. In some cases, Kolmogorov-Smirnov test statistics is something like D = 0.0056 with an associated p-value of 0. Hence, my sample shows very small departures from a theoretical log-normal distribution, but looking at the p-value I reject the null hypothesis that that my sample is drawn from the reference distribution (log-normal). KS-test is performed through R code:               sample.z <- std(sample) # I standardize data to allow for comparisons     LN <- rlnorm(1e5, 0, 1) # theoretical lognormal with mean = 0 and sigma = 1     ks.test(sample.z, LN)      Looking at the QQ-plot I see significant departures in the tails of my distributions. Thus, I was thinking that maybe I'm getting these results because of Pareto tails in my approximately log-normal distributions. Indeed,               library(igraph)     power.law.fit(sample)      confirms this hypothesis by identifying a Pareto tail after a certain lower bound (xmin) with a certain alpha. Now, I would like to show the fit of my Pareto tail in the QQ-plot. How can I do that? Can you suggest other data-visualization methods to stress out the presence of Pareto tails in log-normal distributions? ![PDFs](http://i.stack.imgur.com/56yhc.jpg) ![CDFs](http://i.stack.imgur.com/BKGnm.jpg) ![fitdist output](http://i.stack.imgur.com/7PgM3.jpg)


title: Why do Friedman and Wilcoxon post-hoc contradict each other for my data?

text: I conducted an experiment, testing three conditions with 27 participants in a repeated measures design. Since the data is not normally distributed I used a Friedman test that reported a p-value of 0.1, a slight disappointment given how much time I spent in the lab. Being disappointed but curious, I still performed Wilcoxon's signed rank test as post-hoc analyses for the three pairs. The interesting thing is that these tests reported p-values 0.002, 0.004 and 0.25, respectively. Even after a - conservative - Bonferroni correction two of the p's are still significant. Now, besides the ethical/scientific issues surrounding digging around in the data like that: why are the p-values of the Friedman and the Wilcoxon tests so far apart? And: how do I interpret these findings? They seem to contradict each other.


title: How would I go about analysing MMA statistics and data for future predictions?

text: My maths knowledge is pretty poor but I'm looking to improve on a new project I'm interested in. I have data from thousands of MMA fights/fighters including weights, height, reach etc. and I'm trying to make a calculation to predict future fights. Now I know this is probably a huge task, but I'm thinking of it more as an exercise of how I'd go about doing it, rather than expecting accurate results. How would I go about analyzing data like this? What should I be reading into, and what methods should I be using to get to find patterns in the data etc.? Any resources to start me off would be hugely appreciated. I can remember very little of mathematics from school, so the more simple the better.


title: Hidden Markov Model: Predict observation sequence from state sequence

text: Given a transition matrix, starting probability, means and covariances Is it possible to predict the most likely observed sequence for a given state sequence using the above details? If yes, how? Thanks in advance.


title: Defining contrasts in ANOVA to compare mean differences

text: I would like to run an analysis of variance (parametric or not) with two factors; each factor having two levels. So if x, y is are the levels of the first factor and 1, 2 the levels of the second I would like to calculate whether the difference x1-x2 is significantly different to the difference y1-y2 (with $n_{x1} \neq n_{x2} \neq n_{y1} \neq n_{y2} $), using either SPSS or Matlab. In SPSS I need to specify the LMATRIX coefficients but I have no idea what the coefficients need to be in this case. I am not sure if it possible to do this in Matlab at all. Any ideas or suggestions will be really appreciated. This is related to a question I asked a while ago: Comparison of differences between pairs of samples of unequal size


title: Remove effect of a factor on continuous proportion data using regression in R

text: I have a data set of continuous proportions which depend on a fixed-effect factor, e.g.:               df <- data.frame(        x = c(0.1,0.2,0.12,0.3,0.2,0.35,0.5,0.3,0.45,0.6,0.54,0.77),        f = as.factor(c(rep("A", 6), rep("B", 6)))     )      I'm looking for a method to remove the effect of the factor, e.g. if the residuals were normally distributed, you could achieve this with:               fit <- lm(x ~ f, df)     fit$residuals + fit$coef[1]      What would be the appropriate method for continuous proportion data? I would like the output (after removal of the estimated effect of the factor) to be on the same scale (0-1).


title: Sampling machine learning output to calculate confidence intervals

text: I need to classify large numbers of short answer, free response data from a study with a between-group design. In order to reduce the amount of manual labor costs, I was thinking of manually coding a small sample set, running the rest of the responses through an SVM classifier, and then coding a random sample of the SVM classifier to obtain a classical statistical measures of the automatically coded data-set. The original, and overly verbose, title to this question was, "Is applying random sampling to output from a machine learning classifier a statistically valid way to calculate confidence intervals?" I have already done a conceptual sanity check with a friend of mine who worked with machine-learning algorithms and atmospheric modeling, but I wanted to run it past some real statisticians before I start basing my workflow around this. Thanks!


title: Variable importance from SVM

text: How to obtain a variable (attribute) importance using SVM?


title: Fast computation/estimation of a full-rank linear system

text: This is somewhat related to this question. I have two full rank matrices $A_1,A_2$ each of dimension $p \times p$ and a p vector $y$. These matrices are closely related in the sense that matrix $A_2$ has $p-1$ row in common with matrix $A_1$. Now, I'm interested in the vector $\beta_2$, where $\beta_2=(A_2'A_2)^{-1}(A_2'y)$. My question: Is there is a fast way to get $\beta_2|(\beta_1,A_1,A_2,y)$ in R/matlab?


title: How to get a confidence interval on parameters that were fitted using multiple functions and datasets at once?

text: I use matlab a lot to fit functions to datasets. I then need to determine some kind of error estimate on the fitted parameters. Say I have a dataset of $N$ points which is given by the points $(x_i,y_i),\ i=1,\ldots,N$. I need to fit a function $f(x,P)$ where $P$ is the parameter matrix $$P = \left[\begin{array}{cc} a_1\\\ a_2\\\\\vdots\\\a_k\end{array}\right]$$ where $k$ is the number of parameters to be fitted. I use the _root-mean- squared_ method to find estimates for parameters $a_1,\ldots,a_k$. I build an error function: $$\varepsilon(P,x,y)= \frac{1}{N}\left(\sum\limits_{i=1}^N (f(x,P)-y)^2\right)^{\frac{1}{2}}$$ and then use the matlab function `fminsearch` to get the minimum of this function by optimizing $P$. (I know there are better matlab functions, but bare with me). To get the estimates I calculate the Jacobian matrix $J_f$: $$J_f = \left[\begin{array}{c}\frac{\partial f}{\partial a_1}(P',x) \\\ \frac{\partial f}{\partial a_2}(P',x)\\\ \vdots\end{array}\right]$$ where $P'$ denotes the fitted parameter matrix (best values to minimize the error function $\varepsilon$). I also calculate the variance on the residuals: $$ \sigma_r = \frac{1}{N-k}\sum\limits_{i=1}^N r_i^2$$ with $$ r_i = y_i - f(P',x_i), $$ the residuals. I can get the covariance matrix now: $$ \operatorname{cov}(P) = \sigma_r(J_f^TJ_f)^{-1} $$ where the upper script $T$ stands for transposing and upper script $-1$ stands for the inverse matrix. TO get the error estimates I then need to get the point-wise square root of the covariance matrix: $$s_P = \left[ \begin{array}{c} s_{a_1} \\\ s_{a_2} \\\ \vdots \end{array}\right] = \sqrt{\operatorname{cov}(P)}$$ ## The Question Given the above method is at least a bit correct (I assume it is). Is there a similar approach for the case where not one function $f$ is fitted, but when there are multiple functions $f_1$ and $f_2$? I'll explain what I mean by this. Let's say I have two data sets $(x_i,y_i)$ and $(x'_i,y'_i)$ which comply respectively to functions $f(x,P)$ and $f'(x',P)$ where $P$ is the parameter matrix. While $f$ and $f'$ are different functions, $P$ should be the same for both fits. Eg $f$ might be: $$f(x,P)=a_1(1-exp(-x/a_2)exp(x/a_3)$$ while $f'$ is: $$f'(x',P)=a_4(1+exp(x'/a_2))exp(x'/a_4)$$ I would therefore fit simultaneously with the error function: $$\varepsilon(P,x,y) = \frac{1}{N}\left(\sum\limits_{i=1}^N (f(x,P)-y)^2\right)^{\frac{1}{2}} + \frac{1}{N'}\left(\sum\limits_{i=1}^N (f'(x',P)-y')^2\right)^{\frac{1}{2}}$$ and using `fminsearch` again. This gives me the best estimate for $P$ using all data available (I actually have 4 datasets and 2 functions). How would I estimate the error/certainty bounds on the parameters in this case? (Preferrably using Jacobian again). My english isn't top-notch, nor is my knowledge of statistics and therefore I may be mixing nomenclature a bit. Forgive me for this.


I have started my PhD in statistics this year, and I am looking for your best- practices, advice and (meta-advises) regarding how to grow and become a good academic researcher in the fields of statistics/ML. General thoughts and links are welcomed, but in order to start the ball rolling, here are a bunch of questions gathered from Michael Steele's great article "Advice For Graduate Students in Statistics" (if I am missing important questions, or if some of the questions are meaningless - please also comment on it):   * Papers vs Thesis - how much should one focus on publishing papers during his PhD work? How many papers should one realistically aspire to write?   * In what journals should one strive to get published in? (relevant questions link1, link2)   * How many hours a day should one spend on research (developing/dealing with your research question), and on learning (reading new papers/ attending courses)   * Where does one go to find "hot topic", or even better - a "soon to be hot topic"? (link1, link2)   * Once a "hot topic is found" how should one balance learning the basics of many aspect of the problem, with focusing on one aspect? Obviously these questions are VERY general, and there are many angles for thinking/answering them - I hope to read your perspective on how to think about these general issues. Thanks in advance!


title: How to predict the model fit of a particular % of a sample?

text: I have a large number of projects (about 300). Each project contains a set of test cases. Each test case can result in a success or a failure. The dataset looks like this. Original:               project, passed, total, time, score     p1,      121,    210,  13,    (121/210)     p2,      100,    120,  14,    (100/210)      p3,      100,    120,  11,    (100/120)     p4,      110,    310,  10,    (110/310)     ....p300..      I compute the project score as passed/total Now, I have a method of selecting test cases, which can reduce the time taken for the experiment. The dataset from that after selection would look like this. Selected:               project, passed, total, time     p1,      61,    110,  13     p2,      50,    60,   14     p3,      40,    50,   11     p4,      110,   210,  10     ....p300..      I have a model to relate the original score to the selected score: `original.score = \beta_0 + selected.score * \beta_1` Which gives me a certain `R^2` value. Now, I want to compare this with the case I would have if I had randomly picked 50% of the tests to run. How do I predict the `R^2` if I had selected 50% test cases randomly, and used that in the model? How can I do this in R?


title: Using t-test for comparing Likert responses before and after intervention

text: I am testing a hypothesis that states that participants after an intervention will be more likely to perceive themselves as being at risk for disease. Two questionnaires will be used. One will consist of 5 questions and will measure perceived susceptibility and the the other will consist of 7 questions and will measure perceived disease severity. I think the best way though uncertain to analyze these two questionnaires which are rated on 5-point Likert scales is to use a $t$-test but I am not sure of how or the best way to do it. I am looking for suggestions.


title: Testing that a sample corresponds to an arbitrary distribution

text: I'm writing a set of unit test for a function that draws random numbers from an arbitrary distribution, defined as a PDF (example of such a function can be seen here). I am confused about the proper way of obtaining a good confidence that the function output indeed follows that specified distribution. I know that proving equality is much harder that disproving it, but good approximation will do. An answer on this site suggests using Pearson's chi-squared test for similar purposes in the case of discrete variable. Here, I'm dealing with continuous variable, for which the PDF is defined using paired vectors $X$ and $P$ of final sizes.


title: Suspiciously high pseudo r^2 values in gaussian & ordered regressions (R)

text: I am running 9 regressions with 5 predictors each. The dependent variables are ordered, but have been treated before as count (and have a Gaussian distribution). I ran the 9 regressions using the MASS package, using both 1) glm with family=gaussian, and 2) polr The estimates of the models are very similar. I want to determine pseudo $R^2$ now, and use the `pR2`function of the PSCL package. Again, $R^2$ values are basically identical for Gaussian and ordered models However, the pseudo $R^2$ values are drastically higher than I had expected. m1 0.13 0.27 0.32 0.42 m2 0.17 0.27 0.33 0.43 m3 0.06 0.19 0.36 0.40 m4 0.11 0.21 0.38 0.42 m5 0.14 0.26 0.39 0.45 m6 0.18 0.33 0.37 0.49 m7 0.10 0.28 0.29 0.41 m8 0.04 0.23 0.13 0.28 m9 0.07 0.31 0.12 0.36 The first column is $R^2$ obtained using linear regression models, the second McFadden pseudo $R^2$ , the third Maximum likelihood pseudo $r^2$, the fourth Cragg and Uhler's pseudo $r^2$, which is similar to Nagelkerke. As far as I understand, Cragg and Uhler's pseudo $r^2$ has a range of 0-1, in contrast to the other pseudo $r^2$ values, and should be the one I should use. The values are between 3 and 7 times higher than the $r^2$ values obtained from linear regression models. Now I understand that linear regression doesn't fit very well because of the ordered nature of the variable, but a increase of 3-7 times sounds a bit too much. Also my dependent variables are self reported single questionnaire items of a screening instrument for mental disorders, the predictors covariates like gender and personality traits. I am just skeptical that 40 or 50% variance can be explained by my covariates. Maybe you have input on what might have gone wrong.


title: How to measure usage data across cohorts/segments?

text: Netflix and Google among others use customer usage (number of videos streamed or unique searches conducted) as a proxy for value, or for determining which version of a test is best. I am curious how they would measure this but I'm not great at statistics. If they split their customer base into even groups, the usage within the group will still follow a distribution. So just comparing the mean usage in each group might be misleading if one group has a Bill Gates in it (or someone whose usage dominates the distribution). **What are some proxies I can use to help measure "this segment's usage is higher than the other segment's"?** The distribution I am trying to measure is much more likely to have extreme usage numbers than Netflix or Google whose usage numbers are limited by the number of hours in a day. Would a Q-Q plot be appropriate?


title: Does this quantity related to independence have a name?

text: Obviously events A and B are independent iff Pr$(A\cap B)$ = Pr$(A)$Pr$(B)$. Let's define a related quantity Q: $Q\equiv\frac{\mathrm{Pr}(A\cap B)}{\mathrm{Pr}(A)\mathrm{Pr}(B)}$ So A and B are independent iff Q = 1 (assuming the denominator is nonzero). Does Q actually have a name though? I feel like it refers to some elementary concept that is escaping me right now and that I will feel quite silly for even asking this.


title: Calculating total score for a scale including polytomous and dichotomous items

text: I have a variable that is measured using a 13-item, 4-point Likert scale with the exception of one item (yes $=1$, no $=2)$. How can I calculate the total score for this variable?


title: Determine ratio constant and formula for an equitable score in my algorithm

text: I am working on a algorithm for a software module but decided I was hung up on the mathematical formula so I posted here instead of Programmers. Each scenario will have a team of people competing. Each team can be thought as a unique set of players each of random variable size. Minimum team size is 2. > Eg. Team2 - 9 players Each player can be thought as a binary decision... for simplicity sake lets say that they will randomly grab a hat to wear without looking, either a black hat or a white hat. If every person on the team happens to have a BLACK hat then they will be awarded points. If so much as one person happens to be wearing a white hat then they will not be awarded any points. The above is to give some background and context to the actual problem because it is an easy binary decision in my algorithm. This is incidental to the formula so far. _The important thing to note is that a team can have a random number of players and point awards are all or nothing._ I am trying to determine a formula to determine _how many points to award equitably_. Clearly it is statistically more likely that small teams will have better odds thus the point award should be smaller. Larger teams are riskier so they should get much more bountiful rewards. The starting point is to assume that the smallest possible team (2 players) should receive 10 points. _Let t = 2, p = 10_. I am trying to determine the most statistically equitable formula and an appropriate constant ratio that can be applied in this calculation to determine the most equitable point reward. **EDIT:** In response to a request for more information:   1. `what precisely does it mean to award points "equitably"?` I apologize for the loose term, but equitably in my case would be one where if I were to have a sufficient sample size of Teams, that on average any team regardless of size will have a statistically similar amount of points. I hope this is clear: (Eg. Team1 4 players - 53 points over a year, Team 2 9 players, 51 points over the same year)   2. `Are these shared by team members?` Points are awarded to a team, not to team members. There can be any number of events that a number of teams will participate in. Many teams will not get points. Some will get points. If points are awarded depends on the random player criteria listed above. The amount of points rewarded should depend on the statistical probability that every team member flips a coin and gets heads.   3. `Are they used to determine which team "wins" in the competition?` It is not really a competition and the ones with the most and least points are outliers and not very interesting to me.   4. `Can individual members appear on more than one team?` No. Teams are more interesting than players anyway, at least in this scenario. Teams have a variable number of players and for each event they all essentially flip a coin hoping that everybody gets heads. If the population of events is in the thousands, then I should be able to get a sufficient sample size to determine that teams with smaller amounts of players DO NOT have a statistical advantage at earning points than teams with MORE players. This is what I mean by equitable.


title: Statistical hypothesis testing, one-sided t-test

text: These are the values:               18.2    17.6                                                     ------------------------------------------------------------------------------------------------------------     18.9    18.6                                                     ------------------------------------------------------------------------------------------------------------     20  19.7    19.8    19.6                                             ------------------------------------------------------------------------------------------------------------     21.4    20.6    21                                               ------------------------------------------------------------------------------------------------------------     21.9    22.1    22.2    22.2    21.6                                         ------------------------------------------------------------------------------------------------------------     23.1    22.6    23.2    23  22.9    23.2    23.2    23  22.5    22.5                     ------------------------------------------------------------------------------------------------------------     23.8    24.1    23.8    23.6    24.3    23.7    24  24.2    24.4    23.8    23.5                 ------------------------------------------------------------------------------------------------------------     24.7    24.9    25  25.4    25.1    25.3    25.3    25.2    25.2    24.7    24.6    24.5    24.5         ------------------------------------------------------------------------------------------------------------     25.6    26  25.6    25.6    25.6    25.7    25.6    26.3    26.3    26.2    26  25.9    26.2    26     ------------------------------------------------------------------------------------------------------------     27  27.2    27  27.2    26.7    27  26.6    27.3    26.8    26.6    26.9    27  26.5         ------------------------------------------------------------------------------------------------------------     28.2    27.7    28.1    27.9    27.6    27.7    28.1    28.2    27.5                         ------------------------------------------------------------------------------------------------------------     28.6    29.1    28.7    29.1                                             ------------------------------------------------------------------------------------------------------------     30.3    29.6    29.7    29.5                                             ------------------------------------------------------------------------------------------------------------     30.9    31  31  30.5                                             ------------------------------------------------------------------------------------------------------------     32.3    31.6     ------------------------------------------------------------------------------------------------------------         These are the groups:                - 17.5 -   18.5      - 18.5 -   19.5      - 19.5 -   20.5      - 20.5 -   21.5      - 21.5 -   22.5      - 22.5 -   23.5      - 23.5 -   24.5      - 24.5 -   25.5      - 25.5 -   26.5      - 26.5 -   27.5      - 27.5 -   28.5      - 28.5 -   29.5      - 29.5 -   30.5      - 30.5 -   31.5      - 31.5 -   32.5      $$ \bar{x}= \frac{1}{100} \sum_{i=1}^{100} y_i= \frac{1}{100} \sum_{g=1}^{15} n_g \, y_g = 25.37 $$ $$ \hat{\sigma} = \sqrt{\frac{1}{99} \sum_{g=1}^{15} n_g (y_g - \bar{x})^2 } = 3.1031 $$ Now the EU said that the standard is:   a) at least 25 MPa   b) at least 24 MPa. Decide if the concrete complies the standard with a reliability of 99%.   a) The hypothesis H1 is that the mean is at least 25 MPa.   b) The hypothesis H2 is that the mean is at least 24 MPa. > The results are   >  a) We cannot deny H1; we don't know if the concrete complies the standard > 25 MPa.   >  b) The standard 24 MPa is met with a reliability of 99%.] I know how to calculate the t-score: it's the difference between the measured mean (25.37) and the requested mean (25) divided by standard deviation(3.103) and the whole fraction is multiplied by $\sqrt{n}$ (n=100)... I get:   a) t = 1.192   b) t = 4.41 In the table for the Student's t-distribution for df=99 and $\alpha=.01$, I found t* = 2.365 a) 1.192 < 2.365, so I cannot decide if H1 should be denied or not.   b) 4.41 > 2.365, so I can say that H2 cannot be deined and the standard of 24 MPa is met with a reliability of 99%. So, my question is : When do I deny a hypothesis and when can I be sure that the standard is met?


title: When do we use the Durbin-Watson test?

text: Isn't this test for the determination of auto-correlation of residuals only necessary when time is some sort of a factor in the observed variables? As it is I had a data-set that had one dependent variable (sales) and three independent variables, (unemployment rate, population size and advertizing expense). I find it difficult to imagine how the residuals in this situation can be serially-related without a sequence of observations being specified. My group-mate was unable to explain this to me as I was very skeptical that the Durbin Watson test belonged in our report. Please help me understand


title: Repeated measures ANOVA for an experiment with missing values

text: I have an experiment where several subjects (subjects $= S_1,S_2,...,S_m $) were asked to perform a set of tasks (tasks $= T_1, T_2, T_3,...,T_n$) using both their left ($L$) and right ($R$) arms. Each task for each arm was repeated $r$ times. The response is measured in a variable called 'measure'. Unfortunately, there were some tasks and repetitions missing. I tried using the aov function in R to perform a repeated measures ANOVA analysis, but later found out that this is not appropriate for unbalanced designs. Here is the sample of what I did:               > summary(aov(measure ~ arm*task + Error(subject/(arm*task)), data=all_data))          Error: subject               Df Sum Sq Mean Sq F value Pr(>F)     arm        1 0.3240  0.3240   0.398  0.573     task       4 0.1426  0.0357   0.044  0.994     Residuals  3 2.4397  0.8132                         Error: subject:arm               Df Sum Sq Mean Sq F value Pr(>F)       arm        1 0.0023 0.00234   0.074 0.8027       task       4 0.9972 0.24931   7.941 0.0601 .     arm:task   1 0.0112 0.01117   0.356 0.5928       Residuals  3 0.0942 0.03139                      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Error: subject:task               Df Sum Sq Mean Sq F value Pr(>F)     task       5 0.3898 0.07795   1.652  0.172     arm:task   4 0.1482 0.03706   0.785  0.542     Residuals 35 1.6511 0.04718                         Error: subject:arm:task               Df Sum Sq  Mean Sq F value Pr(>F)     arm:task   5 0.0352 0.007032   0.351  0.878     Residuals 35 0.7013 0.020036                         Error: Within                Df Sum Sq  Mean Sq F value Pr(>F)     Residuals 203  1.288 0.006345                    Warning message:     In aov(measure ~ arm * task + Error(subject/(arm * task)),  :       Error() model is singular      Why am I getting the warning message: "Error() model is singular"? I also found that for unbalanced design it better to use the function anova() from the "car" package. I tried to search for the documentation for anova(), but I am very confused about how to use this function. Which model must I use for testing the following hypothesis: 1\. Within-subjects there is no difference between the left and right arms. 2\. Within-subjects there is no difference between the different tasks. 3\. There is no interaction effect between arm and task within a subject. It would also be very helpful if someone could point me to suitable books that can help me learn about these concepts.


title: Valid way to highlight trends in growing data set

text: I have a summer research project that I'm currently working on that involves face recognition during a video capture session. I currently have an algorithm that returns a ranking of the closest matching results along with a similarity measurement that indicates how closely the input face image was matched. These similarity measurements are based on a comparison to gallery face images on which the algorithm has been trained. As of right now, the application simply continually captures a frame from the video that contains a face image and then queries the algorithm until we get a top result that meets some similarity threshold that the user can set. However, we would like the application to have a "memory" so that it can capture information about the results that are consistently in the top rankings. For example, consider the case where we have a threshold of 65% and we get the following results(only showing top 3 rankings): **Result 1**   1. Jim--62%    2. Mark--61%   3. Jane--53% **Result 2**   1. Mark--62%   2. Brian--55%   3. Sarah--52% **Result 3**   1. Jane--56%   2. Mark--53%   3. Brian--49% **Result 4**   1. Brian--66%   2. Mark--63.5%   3. Sarah--57% The inconsistent top rankings seem to be a result of an idiosyncracy of the algorithm which requires that the orientation/expression of the face closely match its associated training image. This is a crude example, but hopefully I've relayed the concept. In this situation, we would most likely want 'Mark' to be the label for the face since he is consistently placed in the top rankings, but the returned label for the face would be 'Brian'. The main problem is that if we simply aim for a top result with a similarity that meets our baseline target, we fail to take into account any consistent and relevant trends in the history of results, and that is the problem that I was hoping to get some assistance with. The naive approach that I've come up with is simply to continually compute the average of these similarity values, and simply assign the to the face the label with the highest associated current average. Unfortunately, this approach fails to take into account the ranking of the results, which seems like valuable information. One further thing to note is that we are using an eye detection algorithm to preprocess the images, which sometimes fails. If it does fail, we simply guess the locations of the eyes since they are usually in the same general locations in the detected face. If the eye detection does fail, we'd like those results to have less "weight" than those results for which the eye detection was successful. I'll admit that this is a somewhat new domain to me, so I was hoping someone could suggest to me a good statistical/machine learning technique to tackle this problem. Thanks for your time!


title: Error with syntax for BayesGLM in R with Gaussian family

text: I am trying to fit a Bayes glm model with an inverse Gaussian family and `link="identity"`. However, I'm getting the following messsage:               Error in if (iter > 1 & abs(state$dev - devold)/(0.1 + abs(state$dev)) <  :        missing value where TRUE/FALSE needed      My code is as follows:               age<-bayesglm(formula=y~factor(x1)+factor(x2)+factor(x3)+factor(x4)+factor(x5),       family="inverse.gaussian"(link="identity"),prior.scale=Inf,prior.df=Inf,data=data)       Can anyone assist?


title: Why haven't robust (and resistant) statistics replaced classical techniques?

text: When solving business problems using data, it's common that at least one key assumption that under-pins classical statistics is invalid. Most of the time, no one bothers to check those assumptions so you never actually know. For instance, that so many of the common web metrics are "long-tailed" (relative to the normal distribution) is, by now, so well documented that we take it for granted. Another example, online communities--even in communities with thousands of members, it's well-documented that by far the largest share of contribution to/participation in many of these community is attributable to a minuscule group of 'super-contributors.' (E.g., a few months ago, just after the SO API was made available in beta, a _StackOverflow_ member published a brief analysis from data he collected through the API; his conclusion-- _less than one percent of the SO members account for most of the activity on SO_ (presumably asking questions, and answering them), another 1-2% accounted for the rest, and the overwhelming majority of the members do nothing). Distributions of that sort--again more often the rule rather than the exception--are often best modeled with a _power law_ density function. For these type of distributions, even the central limit theorem is problematic to apply. So given the abundance of populations like this of interest to analysts, and given that classical models perform demonstrably poorly on these data, and given that robust and resistant methods have been around for a while (at least 20 years, i believe)--why are they not used more often? (I am also wondering why _I_ don't use them more often, but that's not really a question for _CrossValidated_.) Yes i know that there are textbook chapters devoted entirely to robust statistics and i know there are (a few) R Packages ( _robustbase_ is the one i am familiar with and use), etc. And yet given the obvious advantages of these techniques, they are often clearly the better tools for the job-- ** _why are they not used much more often_**? Shouldn't we expect to see robust (and resistant) statistics used far more often (perhaps even presumptively) compared with the classical analogs? The only substantive (i.e., technical) explanation i have heard is that robust techniques (likewise for resistant methods) lack the power/sensitivity of classical techniques. I don't know if this is indeed true in some cases, but i do know it is not true in many cases. A final word of preemption: yes I know this question does not have a single demonstrably correct answer; very few questions on this Site do. Moreover, this question is a genuine inquiry; it's not a pretext to advance a point of view--i don't have a point of view here, just a question for which i am hoping for some insightful answers.


title: Density function from percentiles (P10, P25, P75, P90, mean and median)?

text: I have percentile data (P10, P25, P75 and P90) for a variable.   I also have the mean and median for each group:               group    mean   median  P10     P25     P75     P90     1        30100  26200   19900   22500   32800   44200     2        38700  36600   28000   31500   44000   52100      **How do I:**   1. Create a probability density function based on these variables.   2. Use that function to give me the % in specific step intervals? (I.e. answering the question: How many out of 100 are in the 30000-31000 interval for group 2?) Thanks.


title: Generating from Dirichlet distribution with the differences in a sequence of ordered uniform

text: First, let's assume that we want to generated from a Dirichlet(1,1,1,1) distribution. Would the following method be correct?:   * generate three variates from a Uniform(0,1). Call them $x_1$, $x_2$, $x_3$.   * then, order these such that $0 \leq x_{(1)} \leq x_{(2)} \leq x_{(3)} \leq 1$   * then, return the differences as our Dirichlet variate: $(x_{(1)}, x_{(2)}-x_{(1)}, x_{(3)}-x_{(2)}, 1-x_{(3)})$ Is this correct? I have a feeling it is correct, but I'm not sure and this does not seem to be the same as either method described on Wikipedia or any other search I did. Maybe it is slow or has other problems, but I'm curious if it is correct. Assuming this is correct, can it be extended to non-uniform Dirichlets, such as Dirichlet(a,b,c,d)? **Extra note** : I am _not_ simply asking how to generate a Dirichlet; there is plenty of information about that already. I'm just curious to see if the method for uniforms can be extended. Is there a more general method that involves drawing from a distribution, then ordering those numbers, then using the gaps?


title: Segmentation of employees

text: I am working on a determining why certain employees cause errors in a company's process and why others do not. I have the employee information, information about the errors they have made and the teams that the employees are in. All the employees in all the teams have caused errors at one point or another. What i want to find out is if there are certain characteristcs about certain employees that causes them to make certain types or number of errors. The sample size is a few thousand over a six months period. What do you guys think is the best approach, i.e. clustering method/ general data mining method?


title: Credible interval for Bayesian posterior of variance and mean, and posterior predictive of normal

text: I am trying to do a "round trip" to verify that my code works with regards to recovering parameters of normal distributions. I'm using this document regarding the normal-gamma distribution: http://www.cs.ubc.ca/~murphyk/Teaching/CS340-Fall07/reading/NG.pdf The basic flow is to generate 100,000 mean and variance parameters from the normal-gamma with prior parameters given. Then for each mean-variance parameter pair, randomly generate 100 samples. Next, using the prior parameters and the data, develop 95% credible intervals for each of the 100,000 mean-variance pairs. And finally create 95% credible intervals for the posterior predictive distribution and compare these to a randomly generated observation. What's interesting is that the mean and variance CIs never seem to contain the actual parameter 95% of the time, but the posteriror predictive interval seems fine. Is this expected? Here is the R code:               rm(list = ls(all = TRUE))          set.seed(115)     no_of_examples = 100000      examples = seq(1:no_of_examples)          #priors     m_0 = 0     k_0 = 6     a_0 = 10     b_0 = 840          means = rnorm(no_of_examples, m_0, sqrt(k_0))     vars = 1/rgamma(no_of_examples, a_0, b_0)          samples = 100      data = matrix(NA, no_of_examples, samples)          for (example in examples){         data[example, ] = rnorm(samples,means[example], sqrt(vars[example]))     }                 # now recover parameters          m_n = rep(NA, no_of_examples)     k_n = rep(NA, no_of_examples)     a_n = rep(NA, no_of_examples)     b_n = rep(NA, no_of_examples)          for (example in examples){         n = samples         x_bar = mean(data[example,])         sample_mean = rep(mean(data[example,]),samples)              m_n[example] = (k_0 * m_0 + n * x_bar)/(k_0 + n)         k_n[example] = (k_0 + n)         a_n[example] = a_0 + n/2         b_n[example] = b_0 + 1/2*(sum((data[example,]-sample_mean)^2))+(k_0 * n * (sample_mean[1] - m_0)^2)/(2 * (k_0 + n))     }          # check posterior of means         conf = matrix(NA,no_of_examples,2)     in_range = rep(NA,no_of_examples)          for (example in examples){         conf[example,] = qt(c(0.025,0.975),2*a_n[example]) * sqrt(b_n[example]/(a_n[example]*k_n[example])) + m_n[example]     }          in_range = ifelse(means > conf[,1],ifelse(means < conf[,2],1,0), 0)          print("% of means within 95% CI:")     print(sum(in_range)/length(in_range))          # check posterior of variances     ci_v = matrix(NA, no_of_examples, 2)          for (example in examples){         ci_v[example,] = (1/qgamma(c(0.975,0.025), a_n[example], b_n[example]))      }          in_range_v = ifelse(vars > ci_v[,1],ifelse(vars < ci_v[,2],1,0), 0)          print("% of variance within 95% CI:")     print(sum(in_range_v)/length(in_range_v))               # check posterior predictive distribution          y_p = rep(NA, no_of_examples)     ci_p = matrix(NA, no_of_examples, 2)          for(example in examples){         # generate an observation         y_p[example] = rnorm(1, means[example], sqrt(vars[example]))         # find 95% CI for posterior predictive         ci_p[example,] = qt(c(0.025,0.975),2*a_n[example]) * sqrt((b_n[example]*(k_n[example] + 1))/(a_n[example]*k_n[example])) + m_n[example]     }          in_range_p = ifelse(y_p > ci_p[,1],ifelse(y_p < ci_p[,2],1,0), 0)     print("% of observations within posterior predictive CI:")     print(sum(in_range_p)/length(in_range_p))      My results are:               [1] "% of means within 95% CI:"     [1] 0.94122     [1] "% of variance within 95% CI:"     [1] 0.93993     [1] "% of observations within posterior predictive CI:"     [1] 0.95123      Is there an error in my code? Or am I completely misunderstanding credible intervals? Update: I edited a small error in the b_n update, and got a little bit closer on the variance CI, but it still seems to be consistently below 95% (when I change the seed). The percentage of means and variances falling within the 95% credible interval both seem to be right around 94%.


title: Quadratic term and variance inflation factor in OLS estimation

text: One question I have carried around with me for a while is related to including quadratic terms for model specifications. I wonder why it is considered ok to include linear and quadratic terms into OLS analysis, despite the fact that the Variance Inflation Factor (VIF) gets high. A model specification that I was taught to be fine to use is the following: $y = c + x + x^2 + \epsilon$ But this model has quite high numerical linear correlation between $x$ and $x^2$ and I do not understand why in this case the VIF does not matter. An example for what I mean using R is the following:               # Load packages     require(car)     requre(lmtest)          # Create artificial data     dat <- data.frame(x = 1:1000, x2 = (1:1000)^2)     dat$y <- dat$x2 + rnorm(1000, 0, 10)          # Run model with linear term     dat.lm1 <- lm(y ~ x, dat)          # Run model with linear and quadratic term     dat.lm2 <- lm(y ~ x + x2, dat)      In this setup, the RESET test suggests using quadratic terms but the VIF is high.               > reset(dat.lm1)     RESET test     data:  dat.lm1      RESET = 27344126617, df1 = 2, df2 = 996, p-value < 2.2e-16          > reset(dat.lm2)     RESET test     data:  dat.lm2      RESET = 0.4404, df1 = 2, df2 = 995, p-value = 0.6439          > vif(dat.lm2)     x       x2      16.03008 16.03008       So the question is: Is this really a model specification that is valid and if yes, why does the VIF not matter in that case?


title: Use fitdist parameters in variables

text: I'm sure this is a simple problem, but I'm not sure how to search and find the answer. I have a set of parameters that I've estimated using 'fitdist':               Fitting of the distribution ' gamma ' by matching moments      Parameters:            estimate     shape 0.2018062     rate  3.1255336      now I want to pull those parameters and store them in variables. But I can't figure out how to separate them. I've tried:               > gamma.fit$estimate         shape      rate      0.2018062 3.1255336       but how do I pull the specific elements? (e.g. shape = 0.2018062 and rate = 3.1255336) I've tried using normal data frame tools but I get an 'incorrect number of dimensions' error. Any ideas would be greatly appreciated, sorry if this is covering old territory. Nick


title: CI for a difference based on independent CIs

text: I'm looking to check my logic here. Say you measure a quantity in group A, and find the mean is 2 and your 95% confidence interval ranges from 1 to 3. Then you measure the same quantity in group B and find a mean of 4 with a 95% confidence interval that ranges from 3.5 to 4.5. Assuming that A & B are independent, what is the 95% confidence interval for the difference between the groups? Presumably you can compute this using standard t-statistics, but I'd like to know if it's also possible to compute an estimate based on the CI's alone. I reason that the lower bound of the CI of the difference should be the minimum credible difference between A & B; that is, the lower bound of the interval for B (3.5) minus upper bound of the interval for A (3), which yields a lower bound for the difference of 0.5. Similarly, the upper bound of the CI of the difference should be the maximum credible difference between A & B; that is, the upper bound of the interval for B (4.5) minus lower bound of the interval for A (1), which yields a lower bound for the difference of 3.5. This reasoning thus yields a confidence interval for the difference that ranges from 0.5 to 3.5. Does that make sense, or is this a case where logic and statistics diverge?


title: Is alpha cumulated in ordinal regression?

text: My outcome is an ordinal variable having 4 levels. I did an ordinal regression and get one universal beta-coefficient describing each level change. If I do 3 logistic regressions on each level change, I would have an alpha- error of 3 * 5% (multiple testing) My question: Is there such a cumulation of alpha error in ordinal regression models (proportional odds or continuation ratio), too?


title: Test whether count of events is significantly greater in a feature list than the general feature population

text: I begin with a matrix of features (~30000 rows) vs cases (columns). I have binary - TRUE or FALSE - data in each cell indicating whether there is an event spanning the feature for that case. Next, I sum each row to get the total number of cases from the sample set that have an event for each feature. I will call this **_c_** for each feature. Next, I retrieve a list of features. How can I test whether **_c_** (the total number of cases with an event for a given feature) is significantly greater than **_c_** among the feature population in general? Furthermore, is there some way of calculating some uncertainty measure of my result based on the fact that only known features associated with a particular process are in the list. NB. My attempt (probably wrong and overly complicated):   1. Bootstrap with replacement 100000 times to get 100000 samples each of size 1023 from the 30000 features. The 1023 features (for example) are excluded from the population from which the features for each sample are drawn.   2. Use `var.test` across **_c_** for each sample _vs_ the 1023 feature list. If the test **does not** produce a significant value, that sample can be used in step 3. Samples that produce a significant value are discarded. This is done to satisfy the assumption of the Wilcoxon test that the variances between the two test samples are similar.   3. Do a Kolmogorov–Smirnov test (`ks.test`) for each sample brought forward _vs_ the 1023 feature list. If the test **does not** produce a significant value, that sample can be used in step 4. Samples that produce a significant value are discarded. This is done to satisfy the assumption of the Wilcoxon test that the distributions between the two test samples are similar.   4. Do a Wilcoxon test (`wilcox.test`) for each sample brought forward _vs_ the 1023 feature list. If the mean (or median?) p-value for these tests is significant, I can say that **_c_** is significantly greater for this set of features than it is for features that do not belong to this set.


title: How to interpret result from mvpart object in R?

text: I used R mvpart package to create a multivariate regression tree. This is part of the output:                         CP               nsplit         rel error    xerror       xstd     1       0.02717093      0         1.0000000 1.0005358 0.03481409     2       0.01302184      2         0.9456581 0.9521266 0.03306820          Node number 1: 3479 observations,    complexity param=0.02717093     Means=12.94,0.5749,9.375,0.72,1.611,0.973,2.153,0.6209,3.307,3.702,2.422,0.3837,1.499,     Summed MSE=1305.19      left son=2 (992 obs) right son=3 (2487 obs)     Primary splits:       Dag           splits as  RRLLRRR, improve=0.02478172, (0 missing)       Hoofdberoep    splits as  RLRLLRRLLLLRLL, improve=0.02313676, (0 missing)       Ploegenstelsel splits as  RRLRRRL, improve=0.02191660, (0 missing)       Werksituatie   splits as  LRLR, improve=0.02179270, (0 missing)       Werkuren       splits as  RLRRRR, improve=0.02130351, (0 missing)      How do I have to interpret the different 'improve values' (0.02478172,0.02313676,...) and how are they related to the complexity parameter (0.02717093)?


title: Logistic regression sub-group size parameters

text: Someone in my lab has a sample of 500 older kids and he wants to investigate what factors are related to the probability that they will bully. Groups:                       bullied  not     bully       22    28     don't      200   250      They created an interaction variable using age*not bullied/bullied and added 14 other "predictors" to the mode. Interaction of age and bullying is sig but the prototypical plot looks weird, the group who were bullied have a 0 slope by age and the line starts and ends at 0. My guess is the the bullied/bullying group is way too small for all the info in the model. IF so, what are the group size parameters or ratio to total group parameters for logistic regression? Thank you very much! ![enter image description here](http://i.stack.imgur.com/zw6Gd.jpg)


title: Normalization prior to cross-validation

text: Does normalizing data (to have zero mean and unity standard deviation) prior to performing a repeated k-fold cross-validation have any negative conquences such as overfitting? Note: this is for a situation where #cases > total #features I am transforming some of my data using a log transform, then normalizing all data as above. I am then performing feature selection. Next I apply the selected features and normalized data to a repeated 10-fold cross-validation to try and estimate generalized classifier performance and am concerned that using all data to normalize may not be appropriate. Should I normalize the test data for each fold using normalizing data obtained from the training data for that fold? Any opinions gratefully received! Apologies if this question seems obvious. **Edit:** On testing this (in line with suggestions below) I found that normalization prior to CV did not make much difference performance-wise when compared to normalization within CV.


title: Where are the t-values in my pairwise comparisons?

text: I have data from a study with three conditions. The independent variable was within-subjects and within-items. I did an F1 and F2 analysis (= twice a repeated measures ANOVA) on the averages, by subject (the participants) and by item (the words). I checked/corrected for sphericity and the means do differ significantly. Now I want to know how the three conditions differ from one another. I have a "Pairwise Comparisons" box and p-values (I checked "Bonferroni" when doing the analyses). But I do not have t-values. How can I know my corresponding t-values? Or should I do some paired-samples t-tests? If so, what does this "Pairwise Comparisons" box tell me? Screenshot of the "PC" box:   ![enter image description here](http://i.stack.imgur.com/cM1Zn.jpg)


title: Finding the instantaneous period of a signal (or calculating period degradation)

text: I have a signal, of varying amplitude, and after a certain point in time, I expect to see the period of this signal lengthen. I am looking for a way to measure the lengthening of this period. Is this a way I can calculate the instantaneous period of the signal below, at any given point? For example, at index 220 I would expect the period to be roughly 40 units, whereas at index 50 I expect it to be 24 units. I imagine the solution relates to either Fourier Transforms or Match Filtering, but I have little practical experience of either of these, so would appreciate any guidance. I am currently using R, but willing to consider other tools if they can help solve this problem. ![Signal](http://i.stack.imgur.com/zFxbz.png)


title: Accelerometer normalization

text: I'm working with data from different accelerometer hardware. Each hardware has a different maximum range, different resolution at which data changes, and different minimum delay after which a new value is emitted. For example for one hardware, the maximum value it can obtain is 10240, emitted after a minimum of 10000 milliseconds, with an increment of 1.0 each time. For another hardware, the maximum value is 2, emitted after 0 milliseconds and incremented by a minimum of 0.038 each time. How cab I possibly normalize the data? Thanks


title: Compare samples with different structure in AB testing

text: Consider two samples in AB testing with different structure. For example for sample 1 we have:   * city1 - 55%   * city2 - 45% (55% of website visitors, who see design1, are from city1) For sample 2 we have:   * city1 - 47%   * city2 - 53% The idea is to "normalize" both samples so they have same structure. I try just to multiply: $x_1/0.47*0.55$ but I think its a bad idea. * * * I try to search by keywords "sample structure" and "sample composition" but with no result. * * * So the question is:   1. Can I get better result (more stable), if I take into account structure of samples?   2. If yes - how to do it?


title: Weighted multiple regression in R with prespecified weights

text: I would like to run a regression of the following form:               Y ~ B1*predictor1 + B2*predictor2 + B3*predictor3      I would like to specify `B1,B2,B3` to hold the values: `0.4,0.2,0.1` Is there way to weight these three predictors by the values of `B1,B2,B3` using lm? So for example I would like to do:               lm(y ~ 0.4*predictor1 + 0.2*predictor2 + 0.1*predictor3)      I know there is the argument weights but that only seems to work with one predictor.


title: High level interpretation of Cramer-Rao bound and Fisher information matrix

text: I am reviewing a manuscript and am struggling to understand why some statistical techniques were chosen, i.e. what information they can give. The paper looks at the effect of predicting a variable such as temperature based on noisy measurements from different sensors. Can you please correct any errors in the following interpretation of the paper's statistical tools:   1. For the temperature prediction, the authors obtain a Cramer Rao Bound (CRB) which defines the smallest error that the prediction might have. A small CRB is desirable, since it means the prediction has a smaller error.   2. The authors decompose the CRB into a numerator and a divisor. The divisor is the Fisher Information Matrix (FIM). A large FIM is desirable, since it means the prediction has a smaller error.    3. The authors lists several situations that would result in a singular FIM, which is undesirable because it means that the temperature prediction isn't valid.


title: Find one slope on dataset with two y-intercepts?

text: I have some analysis that compares relative velocity shifts between regions. I have a representative sample below (in reality there are more regions -- this is for clarity). Assume that I have a dataset with 4 regions -- and the possibility of a relative velocity shift between all four. The only method that I have for determining this relative shift involves a process that forces one of the regions to have a shift of "zero" -- and the other values to shift relative to that. This means that the y-intercept is arbitrary, but the relative y-values between data points is still important. The plot below shows an example of what I'm doing: ![two datasets](http://i.stack.imgur.com/pvcI5.png) This plot shows the result of finding the relative shifts by alternately setting the first region to zero (red); and the last region to zero (blue). As you can tell, the slope is the same between the two measurements. I also have an error estimate for each point not set to be zero. I want to fit for the slope across all of the regions. My question is: is there a statistically robust way to combine this into one fit? Or to fit two linear regressions with independent y-intercepts but tied slopes? Is there another way to think about this entirely?


title: "Build a mean equation for serially correlated data". Should I restrict the search to AR(p), MA(q) or beyond?

text: A problem asks us to build a mean equation for a log return series of a stock: "if there are serial correlations in the log return, build a mean equation for the log returns." Does the mention of serial correlations mean we should restrict our search to AR(p) models? Or should I be exploring whether this is an MA(q) process as well? ![enter image description here](http://i.stack.imgur.com/fTSyh.png)


title: How do I compare multiple timeseries?

text: I'm working with data on litter (fallen leaves, flowers, etc) accumulating under different genera of trees. I have a year's timeseries for each genus, and I'm really more interested in comparisons among genera than trying to fit regressions to seasonality, etc. I'm thinking to use some sort of timeseries tests because subsequent measurements are not fully independent -- this is the standing stock of litter, not the flux of fresh inputs. Data files: all samples, and average by genus/date. In the following plot, points are individual samples, and lines connect genus averages. ![Points are individual samples; lines connect genus averages](http://i.stack.imgur.com/rCkn9.png) There's enough overlap that I'm not going to try to do all pairwise comparisons. Two that are of a priori interest: Genus C is essentially a control, and as you see its litter is much less, especially in the fall. Genus E suffered an insect outbreak in the summer and dropped a lot of green leaves. I'm interested in testing if either of these is statistically significant. I'm working in R. Thank you in advance!


title: Expected peak to trough calculation

text: I am having trouble coding up in R equation 3.1 in the following paper: http://www.intelligenthedgefundinvesting.com/pubs/rb-cb.pdf My function:               expected.drawdown<-function(mn,sd,t){       ret<-mn - ((sd^2)/2) #return R            a<- (ret * sqrt(t))/sd       b<- (ret * t) + (((sd^2) * t)/2)       c<- -(sqrt(t) * (ret + (sd^2))) / sd            x<-(2*ret*pnorm(a)) + (2 * exp(b) * (ret * (sd^2)) * pnorm(c))            x<-x/((2*ret)+ (sd^2))            exp.dd<- 1 - x            return(exp.dd)     }      In page four with the graph, I tried plugging in the values               sd<-0.3162 #/sqrt(252)     mu<-(sd^2) * 0.5     t<-60 #years     expected.drawdown(mu,sd,t)     > 1      while the graph approximates to 0.73 at year 60. Anyone has any idea where I am making mistake?


title: Feature selection and latent variables

text: I would like to know if it is useful (or maybe dangerous) to reduce the number of attributes (by selecting the most informative ones among thousands) before seeking for latent variables or not (in an exploratory perspective). A subsidiary question: in the same case, would it be beneficial to select the most important features for each categorie of features (these can be compressed using an entity-attribute-value model which is not really suitable for data mining) before detecting the latent variables?


title: Class activities/experiments to teach statistical concepts?

text: I intend to give a one hour lecture to teenagers about statistics. I will probably see them only once. This scenario might happen over and over. **I would like to give them some activity to make them experience statistics.** But I am forced to do it with people who do not know anything about probability, statistical inference, exploratory analysis etc. My thought was to go through some simple visualization "tricks" that the media sometimes uses, and debunk that a bit. (please do not give me a link to "how to lie with statistics" :) ) The other idea is to (also) give them an assignment for running an experiment to discover something. For example: discovering if they can detect the difference between coca cola and RC cola. **I am looking for any suggestions for what to do with them, or resources with relevant materials.** Thanks :) (of course, making this a community wiki will be nice)


title: Is it appropriate to use the term "bits" to discuss a log-base-2 likelihood ratio?

text: I'm quite enamoured with likelihood ratios as a means of quantifying relative evidence in scientific endeavours. However, in practice I find that the raw likelihood ratio can get unprintably large, so I've taken to log-transforming them, which has the nice side-benefit of representing evidence for/against the denominator in a symmetric fashion (i.e. the absolute value of the log likelihood ratio represents the strength of evidence and the sign indicates which model, the numerator or denominator, is the supported model). Now, what choice of logarithm base? Most likelihood metrics use log-base-e, but this strikes me as a not very intuition-friendly base. For a while I used log- base-10, which apparently was dubbed the "ban" scale by Alan Turing and has the nice property that one can easily discern relative orders of magnitude of evidence. It recently occurred to me that it might be useful also to employ log-base-2, in which case I thought it might be appropriate to use the term "bit" to refer to the resulting values. For example, a raw likelihood ratio of 16 would transform to 4 bits of evidence for the denominator relative to the numerator. However, I wonder if this use of the term "bit" violates its conventional information theoretic sense. Any thoughts?


title: Choosing models with similar AIC values

text: I'm using a multinomial logistic regression analysis to examine deer behavioural responses to camera traps in terms of 7 predictors (both singly and their interactions). I have found that the model with the lowest AIC value is season and vegetation (AIC = 1005.023). However, the AIC value for season as predictor variable on its own is 1005.103. From what I have read I gather that this means that these 2 models are nominally equivalent in indicating the processes influencing reactions as the AIC value difference is less than 2 or 3. I'm wondering if it is okay to report both AIC models and say that the model season and vegetation indicates that these variables when combined have a strong influence on the underlying processes which influence behavioural responses to camera traps even though the AIC values are nearly the same? The AIC value for vegetation as an isolated predictor variable is 1008. 289. Also, when I run various models based on my hypothses, the best models seem to centre around the 1008 AIC value mark. As the difference in AIC value between the model season and vegetation and these models is approximately 3, does this means that the models with an AIC value of ~1008 are equally as good at indicating the underlying processes influencing behavioural responses as the model season and vegetation? Thanks a mil!


title: Probabilistic latent semantic analysis - objective function

text: In pLSA we are trying to maximize a likelihood function. The likelihood function can be given as: $L = \Pi_i^N \Pi_j^M P(d_i,w_j)^{n(d_i,w_j)}$ Where, M is the number of words and N is the number of documents. I didnt understand this objective function. Here are we trying to maximize the likelihood $P(w|d)$? Could you please explain this objective function (derivation). I didnt understand the power $n(d_i,w_j)$ in the equation. How it came there?


title: Why could centering independent variables change the main effects with moderation?

text: I have a question related to multiple regression and interaction, inspired by this CV thread: Interaction term using centered variables hierarchical regression analysis? What variables should we center? When checking for a moderation effect I do center my independent variables and multiply the centered variables in order to calculate my interaction term. Then I run my regression analysis and check for main and interaction effects, which may show the moderation. If I redo the analysis without centering, apparently the coefficient of determination ($R^2$) does not change but the regression coefficients ($\beta$s) do. That seems clear and logical. What I do not understand: The p-values of the main effects change substantially with centering, although the interaction does not (which is right). So my interpretation of main effects could change dramatically - just determined by centering or not. (It's is still the same data, in both analyses!) Can somebody clarify? - Because that would mean that the option to center my variables would be mandatory and everybody should do it in order to get the same results with the same data. * * * Thanks a lot for distributing to that problem and your comprehensive explanations. Be assured that your help is very much appreciated! For me, the biggest advantage of centering is to avoid multicollinearity. It is still quite confusing to establish a rule, whether to center or not. My impression is, that most of the resources suggest to center, although there are some "risks" when doing it. Again I want to put out the fact, that 2 researchers dealing with the same material and data might conclude different results, because one does centering and the other does not. I just read some part of a book by Bortz (he was a Professor and kind of a Statistics Star in Germany and Europe), and he does not even mention that technique; just points out to be careful in interpreting main effects of variables when they are involved in interactions. After all, when you conduct a regression with one IV, one moderator (or second IV) and a DV, would you recommend to center or not?


title: How do you generate data from the multivariate Marshall Olkin distributed data?

text: I recall having been told that to generated $p$ variate Marshall Olkin distribution, one proceeds as follow: For $p=2$, generate $X_1,X_2,X_3\sim exp(\lambda)$, then $(Y_1,Y_2)\sim MMO$ where $Y_1=\max(X_1,X_2)$ and $Y_2=\max(X_1,X_3)$. But I must have gotten the definition wrong because this doesn't look like a 'MMO' (because there are many observations along $Y_1=Y_2$, which I don't recall was the case):               library(matrixStats)     n<-100     p<-2     a1<-matrix(rexp(n*(p+1),2),n,p+1)     a2<-cbind(rowMaxs(a1[,1:2]),rowMaxs(a1[,c(1,3)]))     plot(a2)      Can anyone help with the algorithm for generating MMO's?


title: Determine a subset of random variables that are most correlated with the whole

text: I have an index variable that is calculated as the sum of $N$ random variables: $Y = \sum_{i=1}^N x_i$. I would like to come up with a subset, $\mathcal{S}$, of the $x_i$ such that $Z \equiv \sum_{i\in\mathcal{S}} x_i$ is as highly correlated with $Y$ as possible. More specifically, I would like to be able to either (a) choose the cardinality of $\mathcal{S}$ and determine the highest possible correlation between $Y$ and $Z$; or (b) choose a desired correlation, and determine how big the set must be. (I don't mean to do this theoretically for a general case; I have data and I'm just interested in finding the subset of variables for the particular variables I have.) If $N$ was small, I could do this with brute force: just calculate all possible subsets and correlations. Alas, in my setting, $N$ is about 4000. I'm looking for suggestions on how to solve this. I wouldn't be surprised if it's been solved in some other setting, but I don't know where. It sounds somewhat related to principal component analysis, so I wouldn't be surprised if the solution involves using eigenvectors.


title: How can I rediscretize my data?

text: Related to my previous question, I have a dataset of 2D points with an associated label (this label can take 6 different values). As suggested in the answers to my other question, this can be modeled as a marked point process (or 6 different point processes), allowing to apply standard tools to study this dataset. I would like to take the approach that I first suggested in my first question, and try to apply PCA on this dataset, to see if the different types of points are correlated or not (i.e. are some types always happening together?). Here's how I want to do it:   1. Split my 2D space in a grid   2. For each cell of that grid, count the number of points of each type. For one cell, this gives me a point in $\mathbb{R}^6 : x_i = (N_1(A_i), N_2(A_i), N_3(A_i), N_4(A_i), N_5(A_i), N_6(A_i))$, where $N_k(A_i)$ is the number of points of the $k^{th}$ point process (corresponding to points of type $k$) in the cell $A_i$   3. Combine all the $x_i$ into a matrix $X \in \mathbb{R}^{6 \times M}$ and apply PCA to this matrix. My question is the following: how do I build the grid? In other words, how do I rediscretize this dataset? Indeed, the intensities of each process are not equal: some types appear more than others. If I just use a regular grid (all cells have the same area), the resulting points will have one or two components that dominates the others. I was thinking of building my grid such that each cell has at most $N$ points, thus bounding the norm of the data points, but I don't think this will solve my "balance" problem. Any suggestion, or pointer to litterature, are appreciated.


title: Maximal depth of a decision tree

text: Probably this is a very basic question, but is it possible to compute the full depth of a decision tree given the sample size of the training set, number of its categories, and the feature dimension?


title: Positive estimator

text: Suppose that one can construct an unbiased estimator $X$ of the quantity $E$, is there a way of getting an unbiased and positive estimator of $E^2$? Indeed, if $X_1$ and $X_2$ are two independent copies of the estimator $X$, the quantity $X_1 X_2$ is an unbiased estimator of $E^2$, but is not positive in general. _Motivation_ : when interested in the marginal distribution $\pi(\theta) = \int_{X} \pi(\theta, x) \, dx$ of an intractable distribution $\pi(\theta,x)$, one can use an MCMC algorithm if one can compute $\pi(\theta)$. This is generally not the case. Nevertheless, it is now well-known that this can also be done if one can only construct an unbiased and positive estimator of $\pi(\theta)$, hence the question.


title: Cost functions for contextual bandits

text: I'm using vowpal wabbit to solve a contextual-bandit problem. I'm showing ads to users, and I have a fair bit of information about the context in which the ad is shown (e.g. who the user is, what site they're on, etc.). This seems to be a pretty classic contextual bandit problem, as described by John Langford. In my situation, there are 2 main responses a user can have to an ad: clicking (possibly multiple times) or not clicking. I have about 1,000 ads I can choose between. Vowpal Wabbit requires a target variable in the form of `action:cost:probability` for each context. In my case, `action` and `probability` are easy to figure out: `action` is the ad I chose to display, and `probability` is the likelihood of choosing that ad given my current policy for showing ads. However, I'm having trouble coming up with a good way to map my payoffs (clicks) to costs. Clicks are obviously good, and multiple clicks on the same ad are also better than single clicks on the same ad. However, not clicking on an ad is neutral: it doesn't actually cost me anything other than the missed opportunity for a click (I'm working in an odd advertising context). Some ideas I've had are:   1. cost = -1 * sign(clicks) + 0 * (not clicked)   2. cost = -1 * clicks + 0 * (not clicked)   3. cost = -1 * sign(clicks) + 0.01 * (not clicked)   4. cost = -1 * clicks + 0.01 * (not clicked) In the case of an action vector of `(0, 1, 5, 0)` the costs from these 4 functions would be:   1. `(0, -1, -1, 0)`   2. `(0, -1, -5, 0)`   3. `(0.01, -1, -1, 0.01)`   4. `(0.01, -1, -5, 0.01)` There are obviously many other ways to represent that `clicks=good` and `no clicks=bad.` In general, how should I be modeling costs for contextual bandit problems in vowpal wabbit? Is it ok to represent benefits as negative costs, or should I re-scale everything such that all costs are positive? Is it ok for relatively neutral actions to have a zero cost, or should I give them a small positive cost to push the model towards the positive actions?


title: Linkage research using longitudial data

text: I have data from three companies that span 10 years. Two sets of variables (mean value of 300 answers pr. year pr. Company), that have the acceptable cronbach's alpha: a) the average employee satisfaction (computed from 3 questions likert from 1-5) b) the average customer satisfaction (computed from 3 questions 0-100)in addition to   C) profit pr. year. Is it possible for me to statistcally test the link between employee satisfaction-customer satisfaction and profit? What is the best way to analyse this data using either SPSS or Excel? Thank you for your answer in advance Marrgie