title: How to interpret underrepresented low p-values in Multiple Testing?

text: This is, I hope, quite a simple question. Being new to the concepts and methods of Multiple Testing, I've been histogramming some of the p-values I've been computing in my genomics research. What I have been seeing is that the distribution of p-values seems to be depleted at low values. Of course, this is horrid for the various FDR, etc methods. But, I'm simply curious as to what this implies about my data? Why would I be getting fewer low p-values than expected from the null hypothesis? The sample size for the test is quite small. I suspect that has something to do with it. These p-values are computed for Cox regression of a clinical variable with a genomic variable. ![histogram of p-values for all tests](http://i.stack.imgur.com/FYbOO.png)


title: Generating quartiles from an imputed variable

text: I am using Multiple Imputation to impute a continuous variable (`X`). I have a question regarding the generation of a new variable, starting from this imputed variable (what, in Stata jargon, would be called a _passive variable_ ). In particular, my aim is to divide this variable `X` into, say, quartiles. Let's suppose that I imputed this variable $n$ times (`X_1, ..., X_n`). What am I supposed to do? Should I create the quartiles "independently" for each of the $n$ imputations? In this case the cut-offs of the quartiles (i.e. the 25th, 50th and 75th percentiles) wouldn't be the same across the $n$ _passive_ variables I created. Should I pool the $n$ imputations, calculate the 25th, 50th and 75th percentiles and use those numbers to divide my $n$ imputed variables into quartiles?


title: How to interpret the variance inflation factors?

text: I have run a Poisson glm and I want to test for multicollinearity in my data. I have used the vif in R and obtained the following result. How can I interpret this?                GVIF Df GVIF^(1/(2*Df))          make    1.083927  1        1.041118     cc      1.093900  1        1.045897     area    1.113460  3        1.018073     insgen  1.029887  1        1.014834     insage  2.484052  2        1.255423     ageveh  1.067675  1        1.033284     marital 2.369840  1        1.539428


title: Linear statistical model with two-class variables

text: Suppose I have a set of (discrete) variables, say $X_1,\dots,X_n$. Each $i$ belongs to either class A or class B. When it belongs to A the contribution is $Y_{A,i}f(X_i)$, and when it belongs to B the contribution is $Y_{B,i}g(X_i)$, so that the final outcome is $Z=\sum_{i=1}^N Y_{A,i}f(X_i) I_{\{i \in A\}} + Y_{B,i}g(X_i) I_{\{i \in B\}}$, where $I_{\{.\}}$ denotes the indicator- function. The values for $X_1,\dots,X_n$ are given, as are the functions $f$ and $g$. Now suppose I have $M$ obervations for $Z$, what method do you advice to (i) Identify whether $i$ is in $A$ or $B$ (ii) Estimate the values for $Y_{A,i}$ (if $i\in A)$ and $Y_{B,i}$ (if $i \in B$). I tried to estimate the coefficients $Y_{A,i}$ and $Y_{B,i}$ using standard linear regression, assuming that $i\in A$ for all $i$ or $i \in B$ for all $i$, but that gives bad results. (iii) What if I have $M$ different inputs $X_{j,1},\dots,X_{j,N}$, for $j=1,\dots,M$ and one observation for $Z$ for each input? Edit: As an example, take the following R-code:               N=50     M=10     X=rbinom(N*M,1,0.5)-rbinom(N*M,1,0.5) #Generate -1, 0, 1 variables     dim(X) <- c(N,M)     pA=1/2 #pB=1-pA     XA = X #f(x)=x     XB = -X^2 + 1 #g(x)=-x^2+1     isA = as.logical(rbinom(M,1,pA))     YA = rnorm(M,1,5) #generate coefficients     YB = rnorm(M,1,5) #generate coefficients     Z = XA[,isA] %*% YA[isA] + XB[,!isA] %*% YB[!isA] #Calculate Z without noise     Z = Z+rnorm(N,sd=sqrt(var(Z))) #add noise       Thus given the observations `Z` and the input data `XA, XB` (these are known) I want to estimate `isA, YA[isA]` and `YB[!isA]`.


title: Random initialization with k-means clustering

text: I read on my machine learning course (on coursera) that random initialization performed several times and then taking the cluster with the lowest cose could help when the number of clusters is "small", but didn't help much for K>>10. Why so?


title: Confusion regarding bias/variance and feature selection

text: I have come across this confusion regarding why we need to add more features when our learning model is bias. I am referring to this lecture from Andrew Ng related to machine learning. I mean suppose I have data point such that they can be modeled better by a quadratic function but instead I choose a linear one So my model is highly biased or underfitted. Suppose my inputs x--> R. I mean one dimensional. Now suppose I use more features and x --> R^n. Now it is possible that even though my inputs cannot be modeled by a linear line in one dimension. But in higher dimension they could fit a high dimensional plane. So, that could fix some issues. This is just an example that I thought of and is no better than using a quadratic function in the low dimensional space. I just need some good examples. I hope I am able to clarify my situation


title: Why are svm training labels +/- 1?

text: I'm new to svm and I've been looking for a svm to use. From all the ones that I've seen, the training label vector is basically a m by 1 vector of 1 and -1. I don't understand why this is so. I was under the assumption that each row of the training vector should be a unique number that labels its respective trained instance. So a m by n matrix, of m training instances each with n features should have m labels that labels each instance. That seems logical to me. But it seems like I'm not quite understanding this. Can someone explain?


title: Does JMP aggregate data in graph-building?

text: I am using JMP to find a relationship between drought index values and yearly corn yields for a 30 year period. I have drought data for seven different indeces and each index ranges from -6 (severe drought) to +6 (extremely high precipitation). The drought data contains monthly drought index values for March-September for every year from 1981 to 2011. The corn yield data contains one number (bushels/acre) for every year from 1981 to 2011. I built two side-by-side graphs in JMP, with one showing corn yield data over the 30 year period and the other showing drought data from one particular index over the same time period. Both are smooth curves and visually line up and seem to have a relationship. My problem is that I have no idea how JMP is processing the drought data. There are eight separate drought index values per year, every year, for each drought index, and only one value per year for the corn yield data. I assumed that JMP was smoothing out the drought curve by taking an average of the eight index values for every year, then using that number in the graph construction but I haven't been able to find any literature that can tell me this for certain.


title: MCMC on a bounded parameter space?

text: I am trying to apply MCMC on a problem, but my priors(in my case they are $\alpha\in[0,1],\beta\in[0,1]$)) are restricted to an area? Can I use normal MCMC and ignore the samples that fall outside of the restricted zone(which in my case is [0,1]^2), i.e. reuse transition function when the new transition falls out of restricted(constrained) area?


title: MCMC for an explicitly uncomputable prior?

text: I am trying to sample from a posterior distribution and I only have an explicit formula for likelihood but I can sample from the prior distribution. How can I sample from the posterior distribution with such a restriction. Is there any specific method? After seeing the answers I've decided to write my exact question to clarify stuff: Its about learning hyper-parameters $\alpha$ and $\beta$ and parameters $\theta_i$ in the following case: $\alpha$ and $\beta$ are uniformly chosen from the perimeter of a square by following vertices: $(0,0),(0,1),(1,0),(1,1)$. Now $\theta_i$ is uniformly chosen from this line. $\theta_i$ it self is the parameter for data $y_i\sim\text{Bin}(n_i,\theta_i)$. In my first attempt, and maybe being foolish I wrote a neat vectored algorithm which would sample from $p(\alpha,\beta,\theta)$ where $\theta=(\theta_1,\theta_2,...)$. But afterwards I realised that it is hardly related to sampling from $p(y|\alpha,\beta,\theta)$ maybe as a result of the answers here. So what I am doing now is that I ignored the whole sampling algorithm I had for the joint priors. To solve the problem is to make a MC random-walk on parameter space $(\alpha,\beta)$ and sub-sampling from it (according to discussion on another question of mine in each step), then sampling from $p(\theta|\alpha,\beta)$ and then calculating the likelihood and then test the new sample according to Metropolis Hasting! I am not even sure this is correct but after my studies, this is the what I can think of!


title: Did anyone know how to fit factor analytic covriance struture either in R or SAS?

text: I would like to make meaningful interpretation from a two-way interaction data using factor analytic covariance structure. I have a genotype x environment matrix and I would like to know which environments account common source of variation from large set of environments. I would like to use mixed effect model and I found factor analytic covariance structure model is promising, but I'm not sure how to proceed. Thanks


title: What is the most appropriate test for a multi-year different group/yr experiment?

text: I've been working on a research project for close to five years now. For my thesis I have to show how "well" did my approach improve things. **Setup** : Every year we use a tool A to brainstorm and negotiate software requirements. The tool was wiki-based and had very low participation from technical and non-technical (client) stakeholders. (Yes we have real living and breathing clients for our class :). I looked at the state of affairs and saw that perhaps social networking based/influenced technologies could help increase participation. So I created tool B to replace tool A. However, tool A was used for the first 2 years and tool B for the latter 3 years. **Environment** : The students changed every year but the overall composition of the class was relatively same (i.e., we had the similar amount of awesome, average and underperforming teams/students). We strive for selecting client projects which are similar level of difficulty and are doable in the duration of the class. They can all be considered projects of the same 'class' (class as in category - same level of complexity etc). Here's my **hypothesis** : Tool B will increase stakeholder participation as compared to Tool A (it was initially in past tense, then changed to present and now to future. Not sure what's right. Keep getting corrections for tense from advisors.) Here are my **measurements** :   1. Client (non-technical) participation using Tool B vs Tool A - via access logs   2. Team (student) participation using Tool B vs Tool A - via access logs + observational data   3. Number of requirements captured/negotiated (new/updated) in Tool A vs Tool B   4. Client surveys (for ascertaining the usefulness of tool A for capturing/negotiating requirements. Already have for Tool B.) One of my advisors suggests I use MANCOVA with 1-3 above as DVs and covariates capturing "sense of complexity of projects", "average number of use-cases per project" and "some metric for team composition/makeup" (which are pretty much the same across the years). Another advisor thinks that simple t-tests would work just fine: That is, I compare the average performance of groups across the years (Group 1 = Tool A; Group 2 = Tool B) and it should be sufficient or maybe an ANOVA at most. Another advior says not to do anything since the data itself is highly skewed i.e., using Tool B has really increased each of the above! He said that doing a statistical test is only to increase the perceived success of the tool and just makes a pompous show of the rigor in analysis. I'm not really sure what would be a good approach here? I'm familiar with t-tests but have never done a MANCOVA ever and am afraid that I may just crunch the numbers and falsify underlying assumptions. What would be an appropriate test for such an experimental design, which is done across multiple years, with different groups, keeping the environment relatively constant? I have many such hypotheses w.r.t. tool B since there are many things that it enables from a process standpoint than what was doable before. It's really confusing with 3 advisors giving different advice and I not being a statistician to be able to decide.


title: When is the differential entropy negative?

text: The definition of entropy for a continuous signal is: $h[f] = \operatorname{E}[-\ln (f(x))] = -\int\limits_{-\infty}^{\infty} f(x) \ln (f(x))\, dx$ According to Wikipedia, it can be negative. When would that happen? As far as I understand, $f(x)$ is always $\in[0,1]$ so $f(x)\cdot ln(f(x))$ can only be negative. What am I missing ?


title: Support vector regression on skewed/high kurtosis data

text: I'm using support vector regression to model some fairly skewed data (with high kurtosis). I've tried modeling the data directly but I'm getting erroneous predictions I think mainly due to the distribution of the data, which is right skewed with very fat tails. I'm pretty sure a few outliers (which are legitimate data points) are affecting the SVR training, and perhaps also in the cross validation, where at the moment I'm optimizing the hyperparameters by minimizing mean-squared error. I've tried to scale my data before applying SVR (e.g. using a sqrt function to reduce the outliers) as well as use a different hyperparameter minimization function (e.g. absolute error), which seems to give better results, but still not very good. I'm curious if anyone has encountered similar problems and how they approached it? Any suggestions and/or alternate methods most welcome.


title: Variance for different values of p

text: > There are two programs in a high school. Boys are a majority (65%) in > program A, and a minority (45%) in program B. There is an equal number of > classes in each of the two programs. > > You enter a class at random, and observe that 55% of the students are boys. > What is your best guess-does the class belong to program A or to program B? Most people chose program A, but program B is slightly more likely. Why is that?


title: Spearman's correlation useful with data that has an upper limit?

text: I am analyzing weather data from 2 sites, and trying to determine if they are correlated. For the most part, I am able to run linear regressions and Pearson's correlations since the data are relatively normal. However, when looking at relative humidity, it obviously has an upper limit in place (100%; it would have a lower limit of 0 as well, but this is not applicable, as I'm only looking at maximum RH). One site reaches 100% RH frequently; the other occasionally does. The Question: Is Spearman's `\rho` appropriate in the case of data with upper limits?


title: Standard errors in weighted least squares

text: My data is a panel of countries by year. Suppose my main RHS variable is a country's GDP and my main LHS variable classifies countries by whether the country is a democracy. Is it desirably to weight each observation by the country's population in a particular year? Since the world population has generally increased over time, does that put too much weight on more recent time periods? To account for that, should I instead weight each country's population by say its population in the initial for each of that country's observations? Also, isn't one purpose of weighting to correct for heteroskedasticity? Should I still adjust standard errors by using robust standard errors? Or, does using robust standard errors obviate the need to weight? Should the weighted and unweighted regressions yield similar results?


title: Monte Carlo variation in bootstrapping?

text: I am obtaining confidence intervals using different bootstrap methods. When comparing them, what is does it mean to use Monte Carlo variation to look at how much the upper and lower limits vary for each bootstrap?


title: Model generation during nested cross validation

text: I'm a little confused here. Appreciate if anyone can help me out. During a _k-fold_ nested CV, I understand for each combination of training fold and testing fold , the training fold will be further split into _k_ subsets and a small CV will be carried out to determine the optimal hyper- parameters. My question is: are the models generated by each training fold the same? If not, which model should be used when apply to productive environment. I know there is a similar thread 'Nested cross validation for model selection'. But I still don't quite get it after reading for several times. Please bear with me as I just step into the field. Thanks in advance. **[Note]** After reading some other threads, I can understand the mechanism below:   1. Create hyper-parameter matrix;   2. Use **sub-training set** to fit model;   3. Use fitted model and **validation set** to select model;   4. Use **test fold** to evaluate the performance of the chosen model. However, my question is triggered by the book 'Data Science for Business'. When it describes 'Nested Cross Validation', it says: > ... before building the model for each fold, we take the training set and > first run an experiment: we run another entire cross-validation on just that > training set to find the value of _C_ estimated to give the best accuracy. > The result of that experiment is used only to set the value of _C_ to build > the actual model for **that** fold of the cross-validation... In my understanding, each fold may get a different value of _C_ which leads to different models. Then which model should I use for productive operation? Please help point out anything wrong in above understanding or if I have any misunderstanding of the book. Thanks a lot.


title: SVM prediction accuracy drops when using Test data

text: I am using the Kaggle Scikit data to learn R. I am using the R e1071 SVM function to predict classes. When I use:               svm(train, trainLabels, scale = TRUE, type = NULL, kernel = "polynomial")      I obtain this level of accuracy on a sample of the Train data:               > table(pred, trainLabels)         trainLabels     pred   0   1        0 478   8     1   12 502      which I interpret as being 98% accurate (8+12) / (478+8+12+502). Though when I use the same prediction model on the Test data, Kaggle returns a **0.82** score, based on classification accuracy. Can you explain why I can get such a different accuracy level?


title: How to create recommender system that integrates both collaborative filtering and content features?

text: I am creating a Recommender System and want to incorporate both the ratings of "similar" users and the features of the items. The output is a predicted rating [0-1].I am considering a Neural Network (to start with). So, the inputs are a combination of the features of the items and the ratings of each user. For item A and user 1, the system could be trained on the combined data, A1. This would be one training example. What if user 1 also rated movie B? Then would the data B1 also be a training example? Is there a problem with repeating the training with user 1's features in this way? Do you have any suggestions about a better way to approach the problem?


title: Regression Statistics model predict

text: Need some help answering this please! The price and the year of manufacture for 124 randomly selected used Mazda cars were recorded in 1992. It is desired to be able to predict the 1992 price of a used Mazda car from its year of manufacture. In order to use the method of modelling using a least-squares line, it was found necessary to first transform the price variable by taking the (natural) logarithm. The relationship between ln price and year of manufacture was found to be suitable for modelling using the least-squares line. The Excel regression output is shown below. (You are reminded that log means natural logarithm, ln, unless otherwise stated.)                   Regression Statistics         Multiple R                    0.901         R Square                      0.811         Adjusted R Square             0.809         Standard Error                0.367         Observations                  124                               Coefficients       Standard Error        t Stat          P-value         Intercept          -314.144             14.121            -22.246           0.000     Year of Manufacturer     0.163               0.007             22.879           0.000      Suppose another used Mazda car was manufactured in 1989. What, to the nearest whole dollar, does this model predict to be its 1992 price?


title: How to test whether the variance of two distributions is different if the distributions are not normal

text: I'm studying two geographically-isolated populations of the same species. Inspecting the distributions, I see that both are bimodal (there's some seasonality to their occurrence), but the peaks in one population are much higher and much narrower (i.e., the variance of the local peaks is smaller). What sort of statistical test would be appropriate to determine whether these differences are significant? To clarify, my y-axis is the number of individuals identified in a trap on a particular day, and the x-axis is Julian day.


title: Variance of product of 2 independent random vector

text: This should be a fairly straight forward but I still couldn't quite get it. Let $X$ be $n$ by $1$ random vector, $Y$ be another $n$ by $1$ and that $X$ and $Y$ are independent. What is $Var(X^tY)=:V_{XY}$? I tried generalizing this formula from the univariate case $Var(xy)=E(x)^2Var(y)+E(y)^2Var(x)+Var(x)Var(y)$ except I hit a bit of a roadblock for the 3rd term above. That is, I get $$ V_{XY}=X^tV_YX+Y^tV_XY\\+E[(X^t-\mu_X^t)(Y-\mu_Y)(Y^t-\mu_Y^t)(X-\mu_X)] $$ So... what is it?


title: Showing spatial and temporal correlation on maps

text: I have data for a network of weather stations across the United States. This gives me a data frame that contains date, latitude, longitude, and some measured value. Assume that data are collected once per day and driven by regional-scale weather (no, we are not going to get into that discussion). I'd like to show graphically how simultaneously-measured values are correlated across time and space. My goal is to show the regional homogeneity (or lack thereof) of the value that is being investigated. ## Data set To start with, I took a group of stations in the region of Massachusetts and Maine. I selected sites by latitude and longitude from an index file that is available on NOAA's FTP site. ![enter image description here](http://i.stack.imgur.com/aZm4N.jpg) Straight away you see one problem: there are lots of sites that have similar identifiers or are very close. FWIW, I identify them using both the USAF and WBAN codes. Looking deeper in to the metadata I saw that they have different coordinates and elevations, and data stop at one site then start at another. So, because I don't know any better, I have to treat them as separate stations. This means the data contains pairs of stations that are very close to each other. ## Preliminary Analysis I tried grouping the data by calendar month and then calculating the ordinary least squares regression between different pairs of data. I then plot the correlation between all pairs as a line connecting the stations (below). The line color shows the value of R2 from the OLS fit. The figure then shows how the 30+ data points from January, February, etc. are correlated between different stations in the area of interest. ![correlation between daily data during each calendar month](http://i.stack.imgur.com/X4YZI.jpg) I've written the underlying codes so that the daily mean is only calculated if there are data points every 6-hour period, so data should be comparable across sites. ### Problems Unfortunately, there is simply too much data to make sense of on one plot. That can't be fixed by reducing the size of the lines. I've tried plotting the correlations between the nearest neighbors in the region, but that turns into a mess very quickly. The facets below show the network without correlation values, using $k$ nearest neighbors from a subset of the stations. This figure was just to test the concept. ![enter image description here](http://i.stack.imgur.com/NWzm2.jpg) The network appears to be too complex, so I think I need to figure out a way to reduce the complexity, or apply some kind of spatial kernel. I am also not sure what is the most appropriate metric to show correlation, but for the intended (non-technical) audience, the correlation coefficient from OLS might just be the simplest to explain. I may need to present some other information like the gradient or standard error as well. ### Questions I'm learning my way into this field and R at the same time, and would appreciate suggestions on:   1. What's the more formal name for what I'm trying to do? Are there some helpful terms that would let me find more literature? My searches are drawing blanks for what must be a common application.   2. Are there more appropriate methods to show the correlation between multiple data sets separated in space?   3. ... in particular, methods that are easy to show results from visually?   4. Are any of these implemented in R?   5. Do any of these approaches lend themselves to automation?


title: Difference of 'centers' of 2 non-normal samples with Mann-Whitney test

text: I have 2 non-normally distributed samples of different sizes (N1~=N2). To evaluate whether there is a significant difference between these samples, I used the Mann Whitney U test (`ranksum` in MATLAB). Now I want to evaluate by how much the populations differ. With normally distributed data, I would just use the difference of the means with a SEM confidence interval. However, for nonparamteric data, my Wikipedia research suggested using the Hodges–Lehmann estimator. My questions are:   1. Is the Hodges–Lehmann estimator indeed correct?   2. How do I calculate confidence intervals for the HL estimator (like the SEM)   3. How can I implement this in MATLAB? Thanks


title: Checking whether a density is exponential family

text: Trying to prove that this doesn't belong to exponential family. $f(y|a)=4\frac{(y+a)}{(1+4a)} ; 0 < y < 1 , a>0$ Here is my approach: $$f(y|a) = 4(y+a)e^{-log(1+4a)}$$ $$f(y|a) = (4y)(1+\frac{a}{y})e^{-log(1+4a)}$$ Comparing it with standard form, $h(y) = 4y$ and $g(a)$ which has to be a function of only $a$, can not be defined in terms of $a$ alone, as $y$ in the $1+\frac{a}{y}$ is inseparable. Is this enough to show that this distribution doesn't belong to exponential family. Please review my approach.


title: Normalize vector-based integer data

text: I'll preface my question by saying I have a very limited knowledge of statistics, and while I've put some thought into this problem, I'm a bit stuck! Onwards... I have a collection of fixed-size vectors of the form $\overline{v} = \{ c_0,\ldots,c_n \}$, where each $c_i \in \mathbb{Z}^*$ and relates to an observation of some unique phenomenon $i$. Each $c_i$ are independent of each other. From a 'training' set $V$ of many such vectors, I would like to compute:   * The average (expected?) count and standard deviation (error? confidence?) for each phenomenon $i$   * Given some new vector $\overline{v}$, some measure of 'distance' from what is 'expected', with respect to $V$. There's one troubling complication: _the observation magnitudes $c_i$ for each vector may be scaled differently_. I do not know what these scaling factors are for each vector. However, I do know that each $c_i$ in a vector have correct relative proportions to one another. I'm thinking that I need to:   1. Normalize the $c_i$ values in each vector $\overline{v}$ so that all vectors in $V$ are directly comparable to each other. I'm currently thinking I should either normalize each vector by standardizing to z-scores or by converting to unit vectors using the L2-norm.    2. Once normalized and comparable, take the average and standard deviation of all $c_i$ values for each $\overline{v} \in V$, giving the expected normalized average $u_i$ and standard deviation $\sigma_i$ for each phenomenon $i$. This gives rise to a vector of 'expected' values with standard deviations for each phenomenon, $\hat{v} = \{u_0 \pm \sigma_0, \ldots, u_n \pm \sigma_n \}$.    3. Use $\hat{v}$ together with the notion of vector distance (in the 'normalized' space) to determine how close/far some new vector $\overline{v}$ is from that is expected. Am I on the right track? Any help much appreciated!


title: Why does Natural Language Processing not fall under Machine Learning domain?

text: I encounter it in many books as well as web. Natural Language Processing and Machine Learning are said to be different subsets of Artificial Intelligence. Why is it? We can achieve results of Natural Language Processing by feeding sound patterns to Machine Learning algorithms. Then, what's the difference?


title: Bayes-factor for testing a null-hypothesis?

text: I heard somewhere, that I can directly test (or gather support for) a null- hypothesis using the Bayes-Factor. In my specific experiment, I hypothesize that an experimental manipulation does not have an impact on some variable but does selectively impact another one. Somehow, simply showing that a t-test gives non-significant results does not seem appropriate (because this can only "not-reject" the H0). My specific problem is as follows: I have two experimental conditions, say medication and control, and I would like to show that estimated parameter-values are the same in both conditions (I measure some data $X$ and have a model with parameters $\alpha,\beta$ giving the likelihood $P(X|\alpha,\beta)$. I measured $N$ subjects in both the control and medication condition (repeated measures). My first approach would be to set up a hierarchical model where $\alpha,\beta$ are distributed for each individual according to some group-level distribution (since both parameters have to be positive and are more likely to be small, I could use an exponential). The group-level distribution would have a uniform prior in a feasible range. I would use this model to sample (MCMC) from the posterior and I think I would be interested in comparing the group-level estimates. However it is unclear to me, how to integrate the two competing hypotheses (H0 being that $\alpha_{med}=\alpha_{cntrl}$ and H1 being that $\alpha_{med}\ne\alpha_{ctrl}$ and the same for $\beta$). So my question is: How can I go from this model to the Bayes-factor? Concretely, I am running my MCMC in the pymc package so any help with concrete code would help a lot!


title: Categorizing text into IPTC subjects

text: Does anyone know where I can find a corpus I can use to train a classifier into IPTC news categories (http://www.iptc.org/site/NewsCodes/) ? A google search was not very useful. Thank you


title: Meta-analysis in R with multiple probes of a microarray

text: Regarding a question posted here, I have the same problem and I can't figure out a solution. I started by doing my meta-analysis in PLINK, but I want the confidence intervals of the meta-analysis, which PLINK does not give. I then tried metagen() function in the meta package, but I'm struggling to write a loop, because the output is a meta class object, and I don't know how to extract the results. I would be happy with any of these solutions: -A way of getting the CI from PLINK or -A way of extracting the results from the meta object to a table. Thanks!


title: Should I perform multiple t-test or an ANOVA?

text: My causal comparative study seeks to three questions related to the hypothesis: a) is there a statistically significant difference between the academic achievement of 6th grade mathematics students based on non-band music participation status; b) is there a statistically significant difference between the academic achievement of 6th grade males based on non-band music participation; and c) is there a statistically significant difference between the academic achievement of 6th grade females based on non-band music participation? Am I correct in stating that I would perform three different t-test to compare means, or do I need to use an ANOVA? Additionally, I have a sample size of approximately 765 students, is this large enough?


title: Learning tasks where semi-supervised naive Bayes is beneficial

text: I am looking into semi-supervised naive Bayes, which appears to be most commonly applied to text classification problems, such as the 20 newsgroups problem. Can anyone recommend some other benchmark datasets (preferably unrelated to problems in text classification) where semi-supervised naive Bayes works well, that I can use in testing my implementation?


title: How to make sense of 6mo worth of weight loss data?

text: -I have ~6mo of data that contains daily weight/caloric intake & #carbs/fats/proteins. I'm trying to figure out a way to make sense of the data. It's not a perfect experiment so there are other factors I should have considered and taken into account some sort of maintenance metabolic conditions (calories/macros to stay at a given weight) but I'm trying to derive something meaningful out of this data. I'm currently a junior/senior in EE with 1 stats course under my belt. My goal is to use what I've learned in class to understand my data better. -What I've done so far: I've found the correlation coefficient between day n-1's calorie, carbs, fat or protein with day n - day n-1' weights, ie: how much a given macro or calorie amount correlates to the following day's change in weight. -The correlation coefficient seems to make sense. The highest positive correlation was with carbs (larger % of diet as carbs, gain weight) and the only negative correlation was with fat (larger % of diet as fat, lose weight).  -Any help would be greatly appreciated, thank you! (if suggestions are beyond the scope of an entry stats course please direct me to where I can learn more about said approach) (I can attach my matlab code if anyone would like)


title: Comparing Bernoulli means across subpopulations in which the number of observed successes may be zero

text: I've got a binary characteristic and a population $S$ with size $n$ and $P[X] = p$ such that $p$ may be small and $n$ is extremely large. Within this population are subpopulations of various sizes $S_0, S_1, \dots, S_k \subset S$. I'd like to be able to select each subpopulation in which $p_i < p$ with some concept of statistical significance. My first inclination is to observe that the standard error on each $p_i$ is $SE_i = \sqrt{\frac{\hat{p_i}(1-\hat{p_i})}{n}}$ and to compare upper bounds on confidence intervals. $\{S_i \; | \; \hat{p_i} + 3 \cdot SE_i < p\}$, for example. But when $\hat{p_i} = 0$, then $SE_i = 0$, and this upper bound is 0 even for the smallest subpopulations (like those where $n_i = 1$). Is there any way to express uncertainty in $p_i$ when $\hat{p_i} = 0$? Maybe through use of $p$ as a prior? **Edit:** It looks like the Jeffreys interval as described in Brown et al. is about what I'm after, though I'm not as-of-yet sure how to apply it.


title: construct the maximum likelihood estimator

text: Let $a_{1},a_{2},a_{3}$ be independent with a normal(0,1) distribution. Define $X_{1},X_{2},X_{3}$ by $X_{1}=a_{1}$, $X_{2}=\theta X_{1}+a_{2}$ and $X_{3}=\theta X_{2}+a_{3}$ Find the MLE for $\theta$ . My attempt: I get that $X_{1}$ has distribution $N(0,1)$, $X_{2}$ is $N(0,\theta^2+1)$ and $X_{3}$ is $N(0,\theta^4+\theta^2+1)$. Then I found the joint densisty and then the log likelihood function. Unfortunately, the expression is quite natsy and partial differentiation would involve simulataneous product, chain and quotient rule. This may be correct but can anyone confirm this is the right approach and the log likelihood function. EDIT; The joint density of $X_{1},X_{2}, X_{3}$ I found to be $\frac{1}{\sqrt{8\pi^3 (\theta^2+1)(\theta^4+\theta^2+1)}}exp(\frac{-1}{2}(x^2+\frac{x^2}{\theta^2+1}+\frac{x^2}{\theta^4+\theta^2+1})$


title: Histogram distance metric for extreme values only

text: I am interested in a histogram comparison method or histogram matching technique that takes into account only the tails of the distribution. Consider the following histograms: Histogram 1: ![enter image description here](http://i.stack.imgur.com/LFdFD.png) Histogram 2: ![enter image description here](http://i.stack.imgur.com/BlFw9.png) While both histograms seem to have a similar distribution, the "extreme" values of histogram 2 (values that fall between 2.4-2.7) are not in the range of histogram 1 extreme values (3-3.5). I am looking for a method (if such exists) that involves a minimal number of parameters, such as tuning percentile values, thresholds etc.


title: Google Prediction API Flaws

text: I have been testing the google prediction api https://developers.google.com/prediction/ It seems to be excellent, i gave it 5 features so it can predict a regression problem, the API fitts the data very well, however as soon as we want to predict on new data we see it diverging. as you will see in the figure below ![enter image description here](http://i.stack.imgur.com/Upkq2.png) So am i missing something here, or does the API just overfit ?


title: K-medoid clustering in python

text: How do I implement k-medoid clustering algorithms like PAM and CLARA in python 2.7? I am currently using Anaconda, and working with ipython 2.7. I have tried scipy.clusters but they don't seem to have the above algorithms. Please help


title: Is there an extension of multinomial process models that can model several dependent variables?

text: In the setting of psychological experiments, where a categorical response has to be given, the same response can, in theory, be generated by different latent processes. For example, in an experiment where participants have to indicate whether a particular word was part of a word list they learned earlier, some correct responses may be based on actual recollection of the word, and some may be due to guessing. Based on the proportions of yes and no responses in different conditions, multinomial process models allow one recover the occurrence probabilities of these latent processes. However, different latent processes may also effect other dependent variables. For example, they could associated with different response latencies: Guessing might require more time than a response based on recollection. **Is there an extension of multinomial process models that takes several dependent variables into account, and if so, is there software for doing that?**


title: Probability integral transformation : Why doesn't this inverse affect inequality?

text: I am reading this problem from DeGroot's "Probability & Statistics" 2nd edition. ![enter image description here](http://i.stack.imgur.com/g9gak.png) ![enter image description here](http://i.stack.imgur.com/2eVGj.png) I can't understand why $$\Pr[G^{-1}[F(X)]\leq z]=\Pr[F(X)\leq G(z)]$$, as nowhere is strict monotonicity assumed. I don't get why $G^{-1}[F(X)]\leq z\iff F(X)\leq G(z)$.


title: R package TSA: how to interpret the IO coefficients output of the arimax function

text: I was playing with the TSA package in R and wanted to test the `arimax` function to the solution provided in Pankratz's _Forecasting with Dynamic Regression Models_ , chapter 8. The savings rate and the function seems to provide similar results as the ones in the book except for the IO weights which are quite different. I bet there is a transformation that I might be missing. Any help on understanding why IO coefficients are so different would be appreciated... the solution states               AO @ t=82,43,89     LS @ t=99     IO @ t=62,55      with Parameters estimates               C = 6.1635     w82 = 2.3346     w99 = -1.5114     w43 = 1.1378     w62 = 1.4574     w55 = -1.4915     w89 = -1.0702     AR1 = 0.7976     MA2 = -0.3762      To fit the model in R, I used (`saving` is the data)               arimax(saving, order = c(1,0,2), fixed=c(NA,0,NA,NA,NA,NA,NA,NA,NA,NA), io=c(55,62),             xreg=data.frame(AO82=1*(seq(saving)==82),                            AO43=1*(seq(saving)==43),                            AO89=1*(seq(saving)==89),                            LS99=1*(seq(saving)>=99)),            method='ML')      The savings rate data is (100 points) 4.9 5.2 5.7 5.7 6.2 6.7 6.9 7.1 6.6 7 6.9 6.4 6.6 6.4 7 7.3 6 6.3 4.8 5.3 5.4 4.7 4.9 4.4 5.1 5.3 6 5.9 5.9 5.6 5.3 4.5 4.7 4.6 4.3 5 5.2 6.2 5.8 6.7 5.7 6.1 7.2 6.5 6.1 6.3 6.4 7 7.6 7.2 7.5 7.8 7.2 7.5 5.6 5.7 4.9 5.1 6.2 6 6.1 7.5 7.8 8 8 8.1 7.6 7.1 6.6 5.6 5.9 6.6 6.8 7.8 7.9 8.7 7.7 7.3 6.7 7.5 6.4 9.7 7.5 7.1 6.4 6 5.7 5 4.2 5.1 5.4 5.1 5.3 5 4.8 4.7 5 5.4 4.3 3.5 here it is my output               > arimax(saving, order = c(1,0,2),fixed=c(NA,0,NA,NA,NA,NA,NA,NA,NA,NA),io=c(55,62),xreg=data.frame(AO82=1*(seq(saving)==82),     + AO43=1*(seq(saving)==43),AO89=1*(seq(saving)==89),LS99=1*(seq(saving)>=99)),method='ML')          Call:     arimax(x = saving, order = c(1, 0, 2), xreg = data.frame(AO82 = 1 * (seq(saving) ==          82), AO43 = 1 * (seq(saving) == 43), AO89 = 1 * (seq(saving) ==          89), LS99 = 1 * (seq(saving) >= 99)), fixed = c(NA, 0, NA, NA, NA, NA,          NA, NA, NA, NA), method = "ML", io = c(55, 62))          Coefficients:              ar1  ma1     ma2  intercept    AO82    AO43     AO89     LS99    IO-55   IO-62           0.7918    0  0.3406     6.0628  2.3800  1.1297  -1.0466  -1.4885  -0.5958  0.5517     s.e.  0.0674    0  0.1060     0.3209  0.3969  0.3780   0.3835   0.5150   0.4044  0.3772          sigma^2 estimated as 0.2611:  log likelihood = -75.57,  aic = 169.14


title: Proportions of different types of nuts in bags of mixed nuts

text: Say I wanted to investigate the proportion of different types of nuts in bags of mixed nuts. So I weigh the nuts by type. Would I be able to use the information on weights and the Chi-squared test to analyze the nuts?


title: Bradley-Terry models by using R?

text: I want your help to Apply the Bradley-Terry models on any particular sport or any other idea of an application by using R. I am not familiar with R, can anyone help me ?


title: How to interpret inconsistent Beta values in different steps of hierarchical regression analysis?

text: I did hierarchical regression analysis on my data due to having moderation effects in my research model. ![enter image description here](http://i.stack.imgur.com/imZJq.jpg) R2 increased from .695 in model1 (main effect only) to .734 in model2 (main &interaction effects)(sig. F change = .000). All the assumptions for the regression analysis have been met. I have two problems with the "coefficients" table: ![enter image description here](http://i.stack.imgur.com/teR03.jpg)   1. As u can see in the table, the insignificant beta value of **ZSC** in model1 became significant in model2! Is it ok? I'm confused! Which value should i consider to reject/accept the related hypothesis? B of model1 (which rejects the hypothesis) or 2 (which confirms it!!)?   2. Although the Beta value for **ZSC_X_CS** is significant, its positive sign is against the hypothesis! it's supposed to have a negative sign according to the literature & also logic! How should i treat this hypothesis? Accept? Reject? Partially accept?!!! TQ. ![enter image description here](http://i.stack.imgur.com/dENFs.jpg)


title: Variance of X/Y

text: Suppose you have two normally distributed, independent random varianbles X and Y, where X has mean $\mu _x$ and variance $\sigma^2 _x$ and Y has mean $\mu _y$ and variance $\sigma^2 _y$. For each, $\mu \gg \sigma^2$ and $\mu \gg 0$. We know from sampling that the distribution of $X/Y$ is approximately normal with mean $\mu _x/\mu _ y$. The question is: what is the approximate variance of this distribution? Below is the sampling histogram: ![enter image description here](http://i.stack.imgur.com/KuArN.png)


title: Is visualization sufficient rationale for transforming data?

text: **Problem** I would like to plot the variance explained by each of 30 parameters, for example as a barplot with a different bar for each parameter, and variance on the y axis: ![alt text](http://i.stack.imgur.com/AERls.png) However, the variances are are strongly skewed toward small values, including 0, as can be seen in the histogram below: ![alt text](http://i.stack.imgur.com/tVWYB.png) If I transform them by $\log(x+1)$, it will be easier to see differences among the small values (histogram and barplot below): ![alt text](http://i.stack.imgur.com/imlTM.png)![alt text](http://i.stack.imgur.com/s8rEU.png) **Question** Plotting on a log-scale is common, but is plotting $\log(x+1)$ similarly reasonable?


title: Analysing questionnaire data

text: I want to test if there is a rivalry among two siblings in a family. I have 15 questions in my study and I let my 100 respondents ( distributed equally to two siblings) ranked them 1 to 15. How should I analyse this data?


title: Time series cross section forecasting with R

text: I have a (I suspect) simple question. I have time series cross section data on voting behaviour in the Council of the European Union (the monthly number of yes, no and abstentions for each member state from 1999 to 2007). So basically the variables are counts, thus a Poisson/negative binomial regression would be appropriate, possibly with lagged dependent variables on the right hand side to control for time dependencies. I have seen papers with people using such negative binomial models to forecast, for instance the number of monthly legislative acts adopted in the future, and I have three questions in this regard:   1. How can i run a negative binomial regression on panel data without making any inferential mistakes?   2. How can I use a negative binomial model with lags to forecast future values of the dependent variable.   3. Can this be done in R? Thomas


title: Proc Calis (or TCalis) and p-values

text: Anyone know why when you run a SEM model in SAS using proc calis (or tcalis) you do not get p-values for the parameter estimates? It does supply a t-value however. Two popular SEM packages in R, 'sem' and 'lavaan', both give p-values for the estimates but they use a Z test statistic.


title: Graph theory -- analysis and visualization

text: I am not sure the subject enters into the CrossValidated interest. You'll tell me. I have to study a graph (from the graph theory) ie. I have a certain number of dots that are connected. I have a table with all the dots and the dots each one is dependant on. (I have also another table with the implications) My questions are:   Is there a good software (or a R package) to study that easily?   Is there an easy way to display the graph? Thanks fRed


title: Low-pass Filter

text: I collected some data with an instruments with 1Hz sampling clock, now I want to low-pass filter the data to separate the mean and fluctuation part (Reynolds decomposition). How can I design a low-pass filter with a cutoff period of 20 minutes?


title: Automatic identification of distribution of data

text: Some statistical tools offer methods of automatic identification of the distribution of data, such as the one shown in this post. Generally, is this approach reliable? If not, what's a better alternative to the problem of working with data from unknown distributions?


title: Is it appropriate to run the Error Correction Model on data which are not I(1)?

text: I have intraday data (frequency = 1 min.) for 6 stocks and 1950 observations per each time series. I checked stationarity for the level data and first difference and it appears that:   1. 5 stocks's level data are non-stationary and 1 is stationary,   2. all 6 stocks's first differences are stationary, meaning that that 1 stock is not I(1) and that 5 stocks are I(1). **My question is:** Is it appropriate to run Engle-Granger methodology, ending with Error Correction Model, on the data pair which includes 1 time series which is I(1) and the other one which is not I(1)?


title: variance of doubly truncated binomial

text: I am in need of variance of doubly truncated binomial distribution (equation number 3.69 on page number 137 of third edition of Johnson, Kemp and Kotz Discrete Probability Distributions). Thanks. Anwer


title: How to deal with RAM limitations when working with big datasets in R?

text: I am currently playing around with the MNIST dataset (http://yann.lecun.com/exdb/mnist/) in R. The training set size is 60000x748 and it seems to drain all my memory even when constructing simple models like logistic regression. My question is: how do you guys usually deal with big datasets in R? And tangent: is it feasible to break the dataset into smaller chunks, construct a model on each, then perform a weighted average on the predicted values?


title: Is it allowed to use a 5% trimmed mean for analyzing data from a creativity task (quantity of ideas)?

text: For my research I conducted a creativity test and measured the quantity of ideas subjects had. Some people are extreme outliers as they have a lot of ideas or only 1 or 2 ideas. Intuitively I wanted to 5% trim my data to obtain a more robust estimation for my regression analyses.   1. Is this justifyable?   2. Is there a paper that can be cited confirming that 5% trimming is a good approach given skewed data?


title: How to create a loop without command "for" in Bugs?

text: I'll try to make an example to make it easier to understand the question.               model{         constant1=(number)     constan2=(number)         for(j in 1:k){ #k defined in data                 a[j] ~ dgamma(c[j],d[j])            b[j]~ dgamma(e[j], f[j])            c[j] ~ dgamma(constan1,constant2)            d[j] ~ dgamma(constan1,constant2)            e[j] ~ dgamma(constan1,constant2)            f[j] ~ dgamma(constan1,constant2)    }              l[i]<-           }      I want to make a likelihood ( l ) that is a sum of elements in 1:k (for each i => a simple loop would solve it..), but unlike R, in bugs I cannot make a for() function or a sum function to make l[i] (multiple definitions..). It's impossible to do something like this in bugs (it accepts only one specification).                for(i in N){ #N is a number defined on data, k also            l[i]<-0          for(j in 1:k){                 l[i]<-l[i]+ d[j]*log(c[j])-log(d[j])+loggam(c[i])+pow(a[j],b[j])+a[j]*t[i] +      delta[i]*a[j]            #it's just an example.. t and delta is in data             }     }      I tried to make this l[i] using functions sum and inprod, bugs doesn't allow to make log() of a function. Also, i keep getting errors when i use inprod with 2 "distributions", like c and d (inprod(d,c)).   And, if I make a l[i] specification directly with a sum all 1:k it gets a error of complex/too many constants. Any idea how to write it in a simple way and make a sum all at once of all 1:k? I am stucked, would be very usefull any comment. Thanks in advance.


title: Mean imputation/estimation of missing data

text: Could someone please refer me to papers that have imputed the mean to missing values of a continuous variable? (i.e. papers that have used mean imputation) I have imputed my missing IMD values using the mean value and I want to see examples of other **papers that have used mean imputation** as well.


title: Estimating abundance using non-normal count data

text: I have sample counts of $n=20$ or $n=7$ taken from right-skewed and zero- inflated populations. The challenge in each case is to use the sample to estimate the total count in that population. Each of the populations are constituted of 300 such counts. As the sample mean ($\times 300$) tends to be unrepresentative can anyone advise a more robust analysis? Some examples of the samples below. Many thanks. n=20 samples:               a <- c(0,    3,   31,    0,    0,  350,   32,  750,   63,  630,   91,   65,    0,              85,    0,  180,   14,  220,   75,  550)          b <- c(96,  210,   50,  195,  500,  450,   45,  180,  150,    0,   22,  250,  900,             190,  500,  110,  300,   25,  280,   40)          c <- c(6,    0,    3,    0,    3,   37,   52,   23,    8,    0,   23,    0,    0,                0,    3,   16,    6,   23,    1,    2)          d <- c(5400, 1800, 4500,  900, 3500, 3100, 5300, 1300, 4800, 4200, 7000,  777,    0,                   5,   42,  850,  700,  550,  500,   18)          e <- c(0,    0,  180,  135,  150,   90,    0,   22,   55,    0,  30,    0,  150,   45,               0,    8,    0,    0,    0,   40)      In addition there was a larger pilot sample of n=76:               p <- c(16,   0,   0,   0,   0,   0,   0,  76,  18,  12,   0, 121,   0,   0,   0,   0,             0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  30,   0,              78,  94,   0,  23,  43,   0,  46,  21, 112,   0,   0,   0,  20,   0,   0, 210,             453,  89, 596, 310, 602,  90, 254, 580,   7,   0,  10,   0,  64, 750,  18, 830,              45, 360,   0,  80,  52,  60,  93, 270,  35,  20,   0,  20)      nb.1. In case of interest these are counts of a plant species' flowers in $1\,\text{m}^2$-quadrats in stands of vegetation of area $300\,\text{m}^2$. The stands were sown with the same seed at the same intensity. Trailing zeroes are where the surveyor scaled up rather than count every individual flower. nb.2. My interpretation of these data is that the quadrats are confounding presence/absence with abundance-where-present. But it seems reasonable to assume that a species has a characteristic probability distribution of abundance-where-present, albeit a J-shaped one. My half-baked solution is to model that distribution (from the pilot sample) as a negative binomial, then combine it as a conjugate prior with the beta-distributed presence/absence of the 7- or 20- quadrat sample. If this approach is justified I would very much appreciate some demo R code.


title: Taking correlation before or after log-transformation of variables

text: Is there a general principle on whether one should compute pearson correlation for two random variables X and Y before taking their log transform or after? Is there a procedure to test which is more appropriate? They yield similar but different values, since log transform is non-linear. Does it depend on whether X or Y are closer to normality after log? If so, why does that matter? And does that mean that one should do a normality test on X and Y versus log(X) and log(Y) and based on that decide whether pearson(x,y) is more appropriate than pearson(log(x),log(y))?


title: Shanken (1992) correction for t-statistics

text: I have done a cross-sectional regression of time-series average returns on estimated Betas (over the same time horizon) to determine average premiums. So far so good. But I was told that the standard t-statistics can be biased, due to the fact that betas are estimated. There is a solution by: > Shanken (1992) - On the estimation of beta-pricing models [Review of > Financial Studies]. It does some small adjustments to the formula of the estimated covariance matrix. However I don't understand how to implement this in R. The paper is also very mathematical, but the solutions are supposed to be easy if you look e.g. at Cochrane Asset Pricing, chap 12 or http://www.uv.es/qf/06006.pdf. I cannot find anything close to that in the original paper though. I think the notation is very different. Does anyone know how to do it, or has done it already? I needed the adjusted formula in my context (in-sample regression), or even better the R-code.


title: MNIST and Vowpal Wabbit

text: I use Vowpal Wabbit 7.3 for 10-class MNIST classification, but can't get any reasonable results. My usage of vw:               ./vw -d mnist_data/mnist.train --oaa 10 -f mnist_data/mnist.model     ./vw -i mnist_data/mnist.model -t mnist_data/mnist.train -p mnist_data/mnist.res      I use 70k digits for training and then test on same data, but predictions always "10". How can I tune VW properly for MNIST?


title: Books about penalized regression

text: I'm looking for textbooks that discuss penalized regression in a methodical manner and present various model selection methods, such as AIC, BIC, Cp as special cases of penalized regression. The textbooks i am looking for don't have to be dedicated to penalized regression, rather this subject can be discussed in a chapter or a section of a more general book. An example of the kind of discussion that I am looking for can be seen in these slides by UCSD instructor Ery Arias-Castro. Books can be in English, French or German. Thanks.


title: Preconditioning the intercept for Logistic Regression

text: I am working on implementing a Logistic Regression model, using the newton-cg and lbfgs optimsers provided by scipy as the backend. I find the problems in which I fit the intercept, to be 50% slower than those in which I don't (which is equivalent to adding a column of ones on X). I am guessing that the reason behind this is because of the difference in between the other features in the training data and the column of ones. Is there are any way I can precondition/scale the intercept, so that this difference is nullified and the model becomes faster. Or is my guess that the fit intercept case is slower because of this reason completely wrong?


title: Why are confidence intervals computed by the delta method so uneven?

text: I plotted a treatment effect on some performance measure over the course of an experiment. Stata has a nice feature called margins and marginsplot to do this. So what I do is run a linear regression, where the time trend is a third degree polynomial, like               reg depvar control##c.time1##c.time1##c.time1 ..., vce(cluster clustvar)      $$y=\beta_1 control\times time+\beta_2 control\times time^2+\beta_3 control\times time^3+\beta_4 time+\beta_5 time^2+\beta_6 time^3+X'\gamma+\epsilon$$ where "control" is a dummy denoting the control group, and "c.time1" is the time variable in months (the ## denote the interactions in Stata). In other words, I fit a time trend for the treatment and one for the control group. The plot is obtained by first computing the average prediction over time via the margins command, and then plotting it: > margins, at(time1==(1(1)39)) over(control) > > marginsplot The end result is this: ![enter image description here](http://i.stack.imgur.com/8LezS.png) As you can see, the width of the confidence intervals looks like a continuous function of the month (x axis). But for example in month 4 or 21 the CI is _extremely_ narrow, whereas it is quite broad for, e.g., month 12. Why is that? It is weird looking, and I can't explain why width in one month is so much narrower than in the other. At first I though that maybe there are many more observations in month 21 that reduce standard errors (and hence narrow CIs), but this is not the case. Why else would the width of confidence intervals differ by so much? The default method to compute these CIs is the delta-method. If I use the option "vce(unconditional)" for margins, the confidence intervals are much broader, but roughly of equal width for each month. Which would be the better choice here?


title: Bayesian MLPs using the MCMC methods - any tricks of the trade?

text: Having used the NETLAB library for MATLAB to implement Bayesian Multi-Layer Perceptron (MLP) neural networks using MacKay's evidence framework, I am now experimenting with Markov Chain Monte Carlo (MCMC) methods using metrop() and hmc(). Does anyone have any advice for successful modelling using this approach (preferably using NETLAB), e.g. how to set the parameters to get good mixing and acceptably low rejection ratios etc? Has anybody used the Metrpolois-Hastings or HMC samplers for applications of MLP neural networks? I'm beginning to think that perhaps everybody just uses the evidence framework!


title: Does quasi-random number generator have a period?

text: I read somewhere, maybe incorrectly, that the Niederreiter quasi-random generator in MKL is 32 bit, and hence as a period of 2^32. This is pretty low, is this correct? This made me wonder if quasi-random number generator have such thing as a period.   My understanding is that, the low discrepancy sequence should go on forever, without repeating itself.   Is this correct? Or, say in another way, is the period of the Mersene-Twistter greater than the period of the Sobol/Niederreiter sequence?


title: How to make prediction models more interpretable?

text: It seems to me as if the go-to technique that is used to make any prediction model more "interpretable" is to reduce the number of input variables that are used in the model. I'm wondering if there are any _other_ big picture approaches that one can take to make these models more interpretable? By a model that is interpretable, I mean a model that is easier to use and understand.


title: Measurement error in dependent variable (type I or II error)

text: This source claims that measurement error in the dependent variable leads to more type I errors (page 2, 4th line of text). However I thought that higher variance of residuals --> larger elements of diagonal of regressor variance covariance matrix --> lower t-statistics / type II errors (failing to reject incorrect nulls).


title: local polynomial regression standard errors

text: I am attempting to find a reference which explains how one computes standard errors for local polynomial regression? Specifically, in R one can use the `loess` function to get a model object and then use the `predict` function to retrieve standard errors. Is there a reference somewhere to what is actually happening? What about in the case when there may be serial correlation in the residuals, one must adjust this using Newey-West type methods, is there a way to use the `sandwich` package to do this as you would for a regular OLS using `lm`?


title: Boosting for regression systems

text: I am quite a newbie in this area:   * What are the boosting methods for regression systems? I know about Gradient boosting; are there any other approaches?    * Are there textbooks or tutorials devoted to this area?


title: How can one evaluate Incremental Clustering Algorithms, in particular the goodness of the clusters formed?

text: I have been studying an incremental clustering algorithm for a large set of data that exhibit an inherent dynamic behavior (that is new data can get added over time and some older data may get deleted too based on the current situation). In this scenario, let us assume, some clusters are formed at one particular time instant. Now when new data comes in, the processing should be limited to only handling this new data and and after metadata extraction, they should be clustered into the already existing clusters (if they are similar, otherwise new cluster must be formed). Now, I am really stumped with the problem of analyzing this algorithm. What kind of cluster validity measures can I use to evaluate the goodness of the clusters formed here? I have only started working in the area and so forgive me if this problem has an obvious answer that I have missed. Thanks in advance. Really hope someone can help me with this.


title: The difference between MSE and MAPE

text: i was wondering what is the differences between Mean Squared Error (MSE) and Mean Absolute Percentage Error (MAPE) in determining the accuracy of a forecast? Which one is better? Thanks


title: Power analysis for binomial data when the null hypothesis is that $p = 0$

text: I'd like to do a power analysis for a single sample from binomial data, with $H_0: p = 0$, vs. $H_1: p = 0.001$, where $p$ is the proportion of successes in the population. If $0 < p <1$, I could use either the normal approximation to binomial, or $\chi^2$-test, but with $p =0$, these both fail. I'd love to know if there is a way to do this analysis. I'd very much appreciate any suggestions, comments, or references. Many thanks!


title: Compare two independent variables

text: I have two variables, one containing project delivery rates for a particular language, and the second containing team size that worked on these projects. I am particularly new to statistics and I want to test the two variables which may or may not be dependent. In other words I want to test the dependency of these two variables. I am using SPSS. Can someone guide me in this regard, which statistical tools can or should be employed for the above scenario?


title: Displaying relationships between nodes

text: I'm trying to create something similar to this. ![enter image description here](http://i.stack.imgur.com/QQVtd.png) So, 3 different Node classes, and a whole bunch of relationships between them. In my case, there should be roughly half of the number of nodes present at most. What I'm looking for is recommendations as to the best way to create a similar type of graph. Spent some time looking at R/GGplot2, but haven't found any solutions so far. I expect it's because I'm not using the correct vocabulary. The posted image was created using a proprietary app that I unfortunately am not able to leverage, otherwise I'd simply use that. Any suggestions/solutions would be fantastic! Thanks in advance.


title: Forecasting customer value as a function of digital engagement?

text: I want measure the value of a digital engagement, or a digitally engaged customer. Therefore I want to build a model to predict customer value as a function of digital engagement with a given digital asset such as my ipad app. Does anyone have thoughts on how to do this, or know of any literature that discusses appropriate regressors for measuring digital engagement? I am thinking my dependent variable is some measure of customer value (profit - cost to serve) and then my regressors are stages of digital engagement (downloaded App, login daily, login multiple times a day) in addition to demographic variables and I am hoping I can forecast customer value from this.


title: Unique or multiple maxima of log-likelihood function?

text: How can I find out if the log-likelihood function has only one global maximum or if it has multiple local maxima?


title: Linear regression 3 variables formula for slope coefficient estimates

text: Okay so I think I found a formula for the coefficient estimates but it is not very concise. It has like 6 sum of squares but it is in a single fraction so it is calculable. I was wondering what the simplest formula would be for estimating the coefficients for a linear regression. So say you had the regression equation: $$y = b_0 + b_1 (x-\overline{x}) + b_2 (z-\overline{z}) + \epsilon$$ where $\overline{x}$ and $\overline{z}$ are averages of $x$ and $z$. So (with no use of matrices) what is the simplest formula you can get that involves the first and second sample moments of $y$, $x$ and $z$ for the coefficient estimates ($\hat{b}_1$ and $\hat{b}_2$)? I just minimized the residual squared and used $s_{xx}$, $s_{xz}$ etc. notation to find the formulas but as I said they are not very nice.


title: Random forest on multi-level/hierarchical-structured data

text: I am quite new to machine learning, CART-techniques and the like, and I hope my naivete isn't too obvious. **How does Random Forest handle multi-level/hierarchical data structures (for example when cross-level interaction is of interest)?** That is, data sets with units of analysis at several hierarchical levels ( _e.g._ , students nested within schools, with data about both the students and the schools). Just as an example, consider a multi-level data set with individuals on the first level ( _e.g._ , with data on voting behavior, demographics etc.) nested within countries at the second level (with country-level data; _e.g._ , population): > `ID voted age female country population`   >  `1 1 19 1 1 53.01`   >  `2 1 23 0 1 53.01`   >  `3 0 43 1 1 53.01`   >  `4 1 27 1 1 53.01`   >  `5 0 67 0 1 53.01`   >  `6 1 34 1 2 47.54`   >  `7 0 54 1 2 47.54`   >  `8 0 22 1 2 47.54`   >  `9 0 78 0 2 47.54`   >  `10 1 52 0 2 47.54`   > Lets say that `voted` is the response/dependent variable and the others are predictor/independent variables. In these types of cases, margins and marginal effects of a variable (partial dependence) for some higher-level variable ( _e.g._ , `population`) for different individual-level variables, etc., could be very interesting. In a case similar to this, `glm` is of course more appropriate -- but when there are many variables, interactions and/or missing values, and/or very large-scale datasets etc., `glm` is not so reliable. Subquestions: _Can Random Forest explicitly handle this type of data structure in some way? If used regardless, what kind of bias does it introduce? If Random Forest is not appropriate, is there any other ensemble-type method that is?_ (Question Random forest on grouped data is perhaps similar, but doesn't really answer this.)


title: Symmetry of a distribution

text: I have received a sample, summarized by a frequency table with groups (i.e. 1-7, 7-10, 10-13, 13-16, etc...) I calculated the mean and median, both were around 15, the mean was 15.3, the median 15.6. The mode was 18. I have calculated the mean by using the middle point of each interval. The median was calculated using an histogram (there is also an identical formula, the one with L0 and L1 if you know which one I mean). The mode was calculated as the middle point of the group with the highest density. I need to say if this distribution is symmetric, symmetric with a bell curve, or not of the options. The first group is wider, and thus in the histogram, it looks only approximately symmetric, and not bell shaped. It looks a bit like left-skewed but just a bit. Is there a rule for the situation where mean=median smaller than mode ?


title: Analysing a rating scale with many predictor variables

text: My response variable is derived from a question in which respondents expressed a preference score (from 1 (do not like) to 10 (like alot)) for a range of scenarios. I would like to analyse this data in relation to a number of socio-economic explanatory variables to find out what factors influence these preferences. I'm not sure what test to use for this. The distribution of the averaged responses is non-normal. I'm wondering if a glm with poisson errors would be suitable? Please help!


title: SEM for Non-Normal Ordinal Data

text: In short, I conducted a satisfaction survey in which surveyees are required to answer on a satisfaction scale from 1 to 7: 285 observations, 37 satisfaction variables. Here is an example of what the scatterplot (jittered) between two variables from the data set looks like (I am working on R): https://drive.google.com/uc?export=download&id=0Bx2Sns2vaI9ycm1tV2pNSWUxQXc Therefore, the data set I am looking into is formed by ordinal data. Nonetheless, the data is non-normal with high skewness and kurtosis values (especially "multivariate"). Accordingly, the test statistics that would result from SEM are most likely all biased and need adjustments. Outlier deletion should be a last resort as there isn't a clear cut between outliers and non outliers (with the exception of a few extreme cases) as shown below: https://drive.google.com/uc?export=download&id=0Bx2Sns2vaI9yOXlYX2VodnVhNWM What would be suggested in this case in order to overcome the problem with non-normality in order to obtain correct statistical tests (whether for the model fit or the t-statistics for the significance of the resulting values)? Thank you very much, Deuterium


title: How to create appropriate number of data points that would be accurate enough for creating regression equations?

text: I found a similar question in this forum. As a rule of thumb, since there are 4 independent variables in my case, I need 4*10=40 data points. However, my question differs slightly, since I want to ask about data generation too. I want to develop a nonlinear regression equation as: $$y=c_{1} A^k B^l C^m D^n$$ by estimating regression parameters: $$c_{1},k, l, m, n$$ To simplify into linear regression, taking log $$\log(y)=log(c_{1})+ klog(A)+llog(B)+m log(C)+nlog(D)$$ where y=dependent variable and A,B,C,D are independent variables. The independent variables vary as follows: $$A=0.01-0.6$$ $$B=50-1000$$ $$C=0.001-10$$ $$D=0.1-10$$ What I can do is, for creating data points, I can take any arbitrary value of 3 out of 4 independent variables from this above range and vary the 4th one using arbitrary interval, and simulate each different value of "y" using my model. I can do this for each variable. But there can be an infinite number of combinations. What is the best way of creating data points of "y" by varying variables: A, B, C, D and how many points would be appropriate in this case? In other words, specifying, Margin of error at 99% confidence for each parameter, and how to choose my data points (A,B,C,D) within the range to sample such that the overall sample size is minimized? Note: There can be an infinite number of combinations of A, B, C, D using arbitrary intervals for each independent variable within their individual range.


title: Poisson Regression Model

text: Kindly explain me how to estimate poisson regression model using iterative weighted least square method. I know it can be easily estimated in any software like Stata, SAS, SPSP etc but I want to estimate it manually. I am unable to understand the weight matrix (W) in the following iterative estimation equation: $b =(X^\mathrm{T}WX)^{-1}X^\mathrm{T}Wz$


title: Why can't likelihood ratio tests be used for non-nested models?

text: More specifically, why do the likelihood ratio tests have asymptotically a $\chi^2$ distribution if the models are nested, but this is no longer the case for the not-nested models? I understand that this follows from the Wilks' theorem, but unfortunately, I don't understand its proof.


title: Enormous SEs in zero-inflated negative binomial regression

text: I have overdispersed count data where the outcome is events (occurrence of a rare disease) and the covariate of interest is season. The unit of analysis is the number of events occurring in a country-season combination. We have 16 countries and 4 seasons repeated across each country, thus 64 data points: ![enter image description here](http://i.stack.imgur.com/7zO6f.png) Since I was suspicious that there may also be an excess of zeroes, I ran several different regression models for comparison: **Negative binomial** ![enter image description here](http://i.stack.imgur.com/FT1l4.png) **Zero-inflated Poisson (ZIP)** ![enter image description here](http://i.stack.imgur.com/BzqEm.png) **Negative binomial hurdle (NBH)** ![enter image description here](http://i.stack.imgur.com/TUk0C.png) **Zero-inflated negative binomial (ZINB)** ![enter image description here](http://i.stack.imgur.com/W9oGx.png) The models yield similar results, except for one thing. **The SEs of ZINB's zero model are enormous.** The other three models have reasonable SEs. There is only one covariate (season) except for the offset term, so no collinearity. The residuals are asymmetric judging by the five-number summary in the output, but that's true for several of the models and it makes sense intuitively. What could be causing this? **EDIT #1** There doesn't seem to be perfect separation in the binomial part of the model. ![enter image description here](http://i.stack.imgur.com/7prR1.png) **EDIT #2** Here are some Pearson residual plots. Definitely not normal, and perhaps heteroscedastic (but the latter, at least, is to be expected). However, I really have no idea what residuals from a ZINB model "should" look like if the model fits. ![enter image description here](http://i.stack.imgur.com/BmB82.png)


title: Show that the random variable $V=X-a$ and $U=a-X$ have same distribution?

text: Given $X$ is a continuous random variable whose density is symmetric about a point $a$. Show that $V=X-a$ and $U=a-X$ have same distribution. * * * $$F_U(u) = P(U \leq u) = P(X-a \leq u) = F_X(a+u)$$ and similarly $$F_W(w) = 1 - F_X(a-w) \longrightarrow f_U(w) = f_X(a+w) = f_X(a-w)$$ by symmetry. Therefore, $f_U(u) = f_X(a+u)$ by changing variable $w$ to $u$, which shows $f_U(u)=f_W(u)$. Is this solution right? Thanks!


title: What is the RMSE of k-Fold Cross Validation?

text: I am testing a neural net to predict numeric values. For that i am using a Training,Validation and Test split. I made a manual 4-Fold CV, this means i am getting 4 RMSE error, each one is the error of the i-th Fold on the test data. How do i get global RMSE of all 4 Folds. Would it be (rmse_1 + rmse_2 + rmse_3 + rmse_4)/(number of all predictions)


title: Variable Transformation using Cumulative Distribution Function (CDF)

text: Consider two different data time-series, **_Data1_** and **_Data2_** , expressed using **inhomogeneous scales (units)**. Each of these two data series is itself a weighted-average of a bunch of **standardized** individual series. I would like to aggregate these two data series into one unique "index", by taking the equal-weight average of the two series at each point in time. However, before doing that, I must put these two series on a common scale. The predominant method is using standardization, i.e. by subtracting the sample mean from the raw values and dividing this difference by the sample standard deviation. However, this assumes that the variables are normally distributed. **_This assumption is clearly violated in the case of the series I am using._** One other method suggests to do the transformation of the raw values on the basis of their empirical cumulative distribution function (CDF) involving the computation of order statistics. I am not sure if this method is valid in my case for the following reason: - Let's take the first data series **_Data1_**. Imagine this series only has 5 data points: (2,-1,4,10,100). The transformed values using order statistic would be: (2/5, 1/5, 3/5, 4/5, 1). When you plot the raw values against time, it seems steady for the first 4 time periods, and there is a jump at t=5. However, when you plot the transformed values against time, there is no evident jump from t=4 to t=5. Therefore, the interpretation of the transformed values is not intuitive. ![enter image description here](http://i.stack.imgur.com/DOMAH.jpg) ![enter image description here](http://i.stack.imgur.com/x7SzD.jpg) Can I still use the order statistic to transform those two data series? If not, would there be any better way of aggregating them? Any help would be appreciated. Thank you!


title: A question about notation of Bayes' Theorem

text: When I studied Bayesian statistics, a question about the notation of Bayes' Theorem came to my mind. Below is the density function version of Bayes' Theorem, where $y$ is data vector and $\theta$ is the parameter vector: $$ p(\theta|y)=\frac{p(y|\theta)p(\theta)}{p(y)} $$ The numerator on the right handside can be written as: $$ p(y,\theta) $$ which is the joint probability distribution of $y$ and $\theta$, then Bayes' theorem could be written as: $$ p(\theta|y)=\frac{p(y,\theta)}{p(y)} $$ Furthermore, $$ p(\theta|y)\propto p(y,\theta) $$ Am I right on this? I think it does not look right. Because the posterior is the proportional to the joint density function. But where is the mistake?


title: Stochastic minimization of AIC criteria

text: I am trying to choose the relevants parameters for a logistic regression with a huge number of parameters. I don't know the business signification of most of them, but I still have to optimize my system for better predictions. I am currently using the function stepAIC of the package MASS. This function is determinist, and if I understood it correctly it works by, at each iteration, adding or removing the parameter that decrease the most the AIC. This function is not so bad, but is not enough for me for the following reasons :   * It is slow. I have more than 100 parameters, and each glm calculation is during 1s, so each iteration is one or two minutes long. It can run for hours, and I want to test more parameters. It is not multithreaded so I can't just use a bigger computer.   * It can stop at a local optima. This is not just a theorical issue, I have seen it stop at a lower AIC with removing some variables.   * It generate only one model. I would prefer to select the ones which have the properties I am looking for (Cramer's V, p-value) and to test them with cross-validation. I tried to do a kind of monte-carlo and to accept a modified model with a probability 2/(1+exp(-diffAIC))-1, but it doesn't work at all. I looked for papers, but I didn't found anything like that.


title: Two Sample Proportion Test - Finite Population

text: If there is a population that we consider fixed (finite) with size N - say the customers at a bank at the end of a given month - and we want to conduct an experiment on this population by randomly sampling m for treatment 1 and k for treatment2 (where m+k <=N), how can you test the difference in proportions between the two treatments? ![enter image description here](http://i.stack.imgur.com/NYEBd.jpg) 1) Traditional tests for two proportions (e.g. z test, chi square) all assume infinite populations in their calculation of p-values correct? Couple of thoughts - do any of these work in the case described?   * Fishers Exact test?   * Permutation tests?   * Bootstrap?


title: R package for feature set algorithm selection

text: I want to train a binary classification NN and part of this will require data pre-processing. However, I have a choice of which pre-processing algorithm to use. Of course I'd like to choose that one which is maximally informative for the purposes of training the NN. I feel sure that there might be some R package that I could us for this selection purpose, so could anyone point me in the right direction?


title: Why does the Cauchy distribution have no mean?

text: From the distribution density function we could identify a mean (=0) for Cauchy distribution just like the graph below shows. But why do we say Cauchy distribution has no mean? ![enter image description here](http://i.stack.imgur.com/zGTLU.png)


title: Difference in Differences: id switches between treatment and control group

text: In my difference in differences model firms $> x$ belong to the treatment group whereas firms$< x$ act as control. I have a two period model:   * In $t_1$ firm $i$ is $> x$ and thus belongs to the treatment group    * In $t_2$ the same firm $i$ is $< x$ and belongs to the control group Do you have any recommendations how to deal with this problem? Additional information: The treatment is a change in a law that only affects firms with more than $x$ employees. Firms smaller than $x$ are not affected by this law. However, a firm that is $<x$ in $t_1$ can certainly exceed the threshold in $t_2$. In this case it will be affected by the law change in $t_2$


title: Parameter fitting with STAN?

text: I have a model that produces data given a set of parameters. Now, given data, I'ld like to find out which parameters of the model are likely. I have an implementation in Matlab that uses Delayed Rejection Adaptive Metropolis for fitting (DRAM toolbox). Basically, DRAM samples parameter values, and tries to minimize an objective function. Not having used STAN before, I was wondering: (1) If and how such this could be implemented in STAN, and, in case it is possible, (2) if one could expect speed improvements as compared to using Matlab, since STAN code is compiled?


title: Matlab: Unable to plot partial autocorrelation plot

text: For a time series I wanted to plot separately the partial auto correlation. Below is the graph for a time series which shows PACF plot of the time series $x$ which I wanted to reproduce: ![acf_desired](http://i.stack.imgur.com/T53D6.png) This shows plot of two time series of two time series - one is a sequence of [-1,1] (variable s) and the other is just numeric numbers (x). The size of data = 256 points. X axis is the delay and Y axis is the PACF coefficients.   This is the plot of PACF that I got from the code: ![ACF which I got](http://i.stack.imgur.com/IFCGk.jpg) This is my code:               N=256;     x(1) = 0.1;     % generate the time series     for i =1 : N          x(i+1) = 4*x(i)*(1-x(i));     end          [Rx, lags] = xcorr(x, 'coeff');     plot (Rx)      where `data = x` time series variable.   I have used the following command to plot plot `x` and the graph shows the plot of two time series for diferent values of $x(1)$. Is this the correct way to do the plot? In Matlab there is a command               parcorr(x)      I even tried with this command, but the result is not any close to any of the two lines in the original plot the above picture. The axis are all different and appear weird. Can somebody please show how to plot PACF for time series?


title: Market / Customer Segmentation - Merging two different segmentations

text: I have a database where each observation is a person. They were questioned on their attitude towards the consumption of X category of product. I have being using K-means to segment this data. I have noticed that people under 19 years old tend to be quite different in their responses to those over 19. I was thinking of dividing the data in `<19` and `>=19` and producing two clusterings and then merging them so as to produce a single report. Does this make any sense from a data mining point-of-view? Is there precedent for doing so?


title: Efficient scale & shape parameter estimation for generalized secant hyperbolic distribution needed

text: the (symmetric) generalized secant hyperbolic distribution GSHD is very flexible but I found not much at all on how to estimate its 3 parameters. Given the location, I need to obtain scale & shape parameter. The generalized Gaussian (type 1) and the Student-t are similar, but the GGD1 is more peaky for the long tail case, and the Student-t does not treat short tail distributions; and even for them I hardly found anything generally applicable. For given degree of freedom or shape parameter it's of course much easier, but this comfortable situation is often not the case. So actually most important for me is the shape parameter estimation, but many methods like moment or L-moment fitting is not applicable here of course due to infinite moments of the GSHD. Bye Stephan


title: Weighting observations in binary logistic regression

text: I am working in program R. I am modeling the incidence of flight in a seabird in relation to distance to the nearest ship (potential disturbance, range = 0 to 74 km from the bird). 1= flight during observation, 0 = no flight. The bird does fly with some unknown probability when no ships are present or really "far" away. I am trying to find this really far distance and associated probability of flight using binary logistic regression. Model = Flight ~ ship distance. Other variables were explored but fell out with stepwise selection. During exploratory analysis I truncated the data down only looking at smaller distances from the ship (20, 15, 10 km). These models are highly significant and predict that as the ship gets closer the probability of flight increases. However when I include all the data (out to 74 km) the intercept is significant (and predicts the true % of observed flight events) but the slope term is non-significant. Can I use a weighting scheme to give more weight to observations when the ship was closer? Thanks. Edit: I am working through the suggestions made by @Scortchi and @Underminer. Here is a plot of a loess smooth on the observed data to better help visualize the pattern. ![Loess smooth of probability of flight as a function of distance to nearest ship](http://i.stack.imgur.com/ZabQh.jpg) The distance to the ship data does not discriminate between approaching ships and departing ships it is just a straight line measure to the nearest ship. The dip in the probability of flight at 8.5 I believe can be attributed to "unaffected" birds that did not fly as the ship passed by them. So as the ship departs and gets further from the observation site we were more like to be observing birds that for whatever reason did not fly when the ship passed and are less likely to fly for "naturally occurring" reason. As additional birds fill back into the observation area the "baseline" flushing rate is resumed and birds start to fly at "normal" probabilities.


title: Friedman test vs Wilcoxon test for binary IV

text: I'm trying to assess performance of a supervised machine learning classification algorithm. The observations fall into nominal classes (2 for the time being, however I'd like to generalize this to multi-class problems), drawn from a population of 99 subjects. One of the questions I'd like to be able to answer is, if the algorithm exhibits a significant difference in classification accuracy between the input classes. For the binary classification case I am comparing mean accuracy between the classes across subjects using a paired Wilcoxon test (since the underlying distribution is non-normal). In order to generalize this procedure to multi-class problems I inteded to use a Friedman test. However, the p values obtained by those two procedures in case of a binary IV vary wildly, with the Wilcoxon test yielding `p < .001` whereas `p = .25` for the Friedman test. This leads me to believe I have a fundamental misunderstanding of the structure of the Friedman test. Is it not appropriate to use a Friedman test in this case to compare the outcome of the repeated measures of the accuracy across all subjects? My R code to obtain those results (`subject` is the subject identifier, `acc` the accuracy DV and `expected` the observation class IV):               > head(subject.accuracy, n=10)        subject expected        acc     1       10     none 0.97826087     2       10     high 0.55319149     3      101     none 1.00000000     4      101     high 0.68085106     5      103     none 0.97826087     6      103     high 1.00000000     7      104     none 1.00000000     8      104     high 0.08510638     9      105     none 0.95121951     10     105     high 1.00000000     > ddply(subject.accuracy, .(expected), summarise, mean.acc = mean(acc), se.acc = sd(acc)/sqrt(length(acc)))       expected  mean.acc     se.acc     1     none 0.9750619 0.00317064     2     high 0.7571259 0.03491149     > wilcox.test(acc ~ expected, subject.accuracy, paired=T)              Wilcoxon signed rank test with continuity correction          data:  acc by expected     V = 3125.5, p-value = 0.0003101     alternative hypothesis: true location shift is not equal to 0          > friedman.test(acc ~ expected | subject, subject.accuracy)              Friedman rank sum test          data:  acc and expected and subject     Friedman chi-squared = 1.3011, df = 1, p-value = 0.254


title: How to produce the minimum forecast error using R?

text: Considering that we want to use optimize() on the interval [0,1] how can I write an R code for finding the value of β that produces the minimum forecast error without using external packages like `forecast`? ![enter image description here](http://i.stack.imgur.com/DE1Hn.png) ![enter image description here](http://i.stack.imgur.com/gb7Bh.png) For simplicity assume that: ![enter image description here](http://i.stack.imgur.com/B9c2A.png) I want to use the following package:               > require(datasets)     > str(nhtemp)      Time-Series [1:60] from 1912 to 1971: 49.9 52.3 49.4 51.1 49.4 47.9 49.8 50.9 49.3 51.9 ...      in which `nhtemp` is the `Yearly Average Temperatures in New Haven CT`.


title: Normal distribution probability calculation

text: I am stuck with a stats question..... help please > A television manufacturer is studying television remote control usage. One > of the criteria they are measuring is the distance at which people attempt > to activate the television set with the remote control. They have discovered > that the activation distances are normally distributed with a mean distance > of 6 feet and a standard deviation of 2.75 feet. If a remote’s maximum range > is 10 feet, what percentage of the time will users attempts to operate the > remote outside of its operating limit?


title: Decompose LibSVM model into binary classifiers

text: I have a multi-class .model trained with LibSVM. Is there a way to decompose this model into different binary .model files? For example, if I have a .model trained for $n$ classes, I would like to have $(n-1)/2$ binary classifiers separately. As far as I know, LibSVM trains one vs one in order to obtain the multi-class classifier.


title: comparing two groups on ordinal data with covariates

text: Hi I am using SPSS to compare two groups on ordingal data (frequency of substance use): never, 1-5/week, 6-10 week etc. The substance use data is highly skewed and non-transformable so ANOVA is out. I tried a Mann-Whitney U-test but SPSS gives error that ordinal data can not be entered. Moreover I need to covary for age, gender, and a few scale variables. Any ideas? I thought of ordinal regression, but wasn't sure if I need to run this for each group separately and if so would I compare statistically the F-values. Thanks for any help with this.


title: Can lack of main effect and lack of interaction be caused by the same confound?

text: Can the lack of main effect have the same underlying cause as the lack of interaction in 2-way ANOVA? My results failed to reach significance for variables of gender and language. Is it possible that a confounding variable which affected language (and thus made it impossible for main effect of language to reach significance) has influenced the interaction in the same way? Or am I being lazy and not analysing my design deeply enough?


title: Choosing which data-point to label (active learning)

text: For an online unsupervised learning algorithm, data-points are learned sequentially. The performance may improve if in addition to the unlabelled data we have some labelled data-points (i.e. semi-supervised learning with a small amount of labelled data). In this situation, it may be attractive to let the algorithm decide which data-points to label, that is, when the algorithm get a new data-point, it may actively decide to request the label of this data-point to the user, because it judges that it is an "important" example. As far as I know, this is called "active learning". My question is: how or in which situations can a learning algorithm decide that the current example (data-point) is important, and thus requesting its label ? that is, which measures or criteria can allow us to know if we should or not request the label of a given data-point ?


title: Small n: non-parametric of parametric tests?

text: I am in a situation where I have sample with $n=19$, and the the data are not normally distributed (according to the results of a Kolmogorov-Smirnov test) for all the variables. The data were collected using Likert scales in a questionnaire. So now I am unsure whether or not I can use a parametric test or whether I need to use non-parametric tests. I have found books stating that if you have a small $n$, you should always use non-parametric tests. However I have also found citations stating that the choice between parametric and non-parametric tests depends on the level of your data (Likert can be seen as nominal), so I should use parametric tests.


title: Multivariate meta analysis and metaregression with multi arm randomized controlled studies

text: Please, i need to do a network meta-analysis and metaregression or a multivariate meta analysis and metaregression. I have multi arm randomized controlled studies. I tried the metafor package # Dummy database with controlled multiarm treatment studies               Study = c(1,2,3,1:2,4:10,1,2,7)     Group1 = as.factor(c(rep(1,3),rep(2,9),rep(3,3)))     numberInt= c(rep(30,5),rep(28,7),rep(40,3))     InterventionMean1  = c(rnorm(3, mean = 350, sd = 2), rnorm(9, mean = 540, sd = 2),rnorm(3, mean = 860, sd = 2))     InterventionSD  = rnorm(15, mean = 80, sd = 15)     numberCont= c(rep(29,4), 28,rep(30,7),rep(42,3))     ControlMean  = rnorm(15, mean = 230, sd = 2)     ControlSD  = rnorm(15, mean = 55, sd = 9)     Dose = c(40,50,40,100,100,100,100,100,100,100,100,100,200,200,200)     x=as.factor(c(1,2,1,rep(1,5),rep(2,3),rep(1,4)))     data = data.frame(cbind(Study,Group, numberInt,InterventionMean1,InterventionSD,                  numberCont, ControlMean, ControlSD,Dose,x))      # x and Dose are covariates # Loading libraries               library(metafor)     library(Matrix)      # Effect size : mean difference dataEscalc = escalc(measure = "MD", m1i = InterventionMean1, sd1i = InterventionSD, n1i = numberInt, m2i = ControlMean, sd2i = ControlSD, n2i = numberCont, data=data) # Number of patients per studies               dataEscalc$Ni <- unlist(lapply(split(data, data$Study), function(x) rep(sum(x$numberInt) + x$numberCont[1], each=nrow(x))))      # Variance-covariance matrix                   calc.v <- function(x) {         v <- matrix(vi/x$numberCont[1] + outer(x$yi, x$yi, "*")/(2*x$Ni[1]), nrow=nrow(x), ncol=nrow(x))           diag(v) <- x$vi               v               }         V <- lapply(split(dataEscalc, dataEscalc$Study), calc.v)     V <- as.matrix(bdiag(V))     V      # Conducting multivariate meta anlaysis names(dataEscalc)               res <- rma.mv(yi, V, mods = ~ factor(Group1)  - 1,                   random = ~ factor(Group1) | Study,                   data=dataEscalc, method="REML")      > I get this error message > > Error in rma.mv(yi, V, mods = ~factor(Group1) - 1, random = ~factor(Group1) > | : Error during optimization. # Conducting multivariate metaregression               res <- rma.mv(yi, V, mods = ~ factor(Group1) + Group1:I(Dose) + Group1:I(x) - 1,                   random = ~ factor(Group1) | Study,                   data=dataEscalc, method="REML")      > I get this error message > > Error in rma.mv(yi, V, mods = ~factor(Group1) + Group1:I(Dose) + Group1:I(x) > - : Model matrix not of full rank. Cannot fit model. I followed the instructions from http://www.metafor- project.org/doku.php/analyses:gleser2009 http://www.metafor- project.org/doku.php/analyses:vanhouwelingen2002 Thanks in advance


title: A textbook error w.r.t structure and pattern loadings

text: I have this picture in Lattin representing structure and pattern loadings in factor analysis. If $Z$ (an observed variable) $=w_1 F_1+w_2 F_2$ (according to factor model), then the _pattern_ loadings of $Z$ on $F_1$ and $F_2$ must be $w_1$ and $w_2$ (see the diagram). Using some geometry, the _structure_ "loadings" (the correlations between $Z$ and the factors) of $Z$ on $F_1$ and $F_2$ must be $x_1=w_1+w_2 cos\psi$ and $x_2=w_2+w_1 cos\psi$ respectively. ![Lattin](http://i.stack.imgur.com/70wwQ.png) However the formula in the bottom is the opposite of whatever the diagram claims. Which is correct: the diagram or the formula?


title: Greek Gamma in a moment generating function

text: I'm puzzled: what does a Greek $\Gamma$ mean in statistics? For example, here http://en.wikipedia.org/wiki/Weibull_distribution in the definition of the mean.