title: Clustering quality measure

text: I have a clustering algorithm (not k-means) with input parameter $k$ (number of clusters). After performing clustering I'd like to get some quantitative measure of quality of this clustering. The clustering algorithm has one important property. For $k=2$ if I feed $N$ data points without any significant distinction among them to this algorithm as a result I will get one cluster containing $N-1$ data points and one cluster with $1$ data point. Obviously this is not what I want. So I want to calculate this quality measure to estimate reasonability of this clustering. Ideally I will be able to compare this measures for different $k$. So I will run clustering in the range of $k$ and choose the one with the best quality. How do I calculate such quality measure? UPDATE: Here's an example when $(N-1, 1)$ is a bad clustering. Let's say there are 3 points on a plane forming equilateral triangle. Splitting these points into 2 clusters is obviously worse than splitting them into 1 or 3 clusters.


title: Balanced sampling for network training?

text: Hi I have a very imbalanced sample set, e.g. 99% true and 1% false. Is it reasonable to select a balanced subset with a 50/50-distribution for neural network trainining? the reason for this is, that i guess training on the original dataset may induce a bias on the true-samples. Can you suggest my some literature that covers this topic especially for neural netowrks? Thank you very much, Julian


title: Inconsistent results using Bayesian Model Averaging for model selection?

text: I have attempted to use Bayesian Model Averaging (BMA) as method to give guidance in the decision to include/omit factors in the linear factor stock return model: $$r^i = X \beta^{i} + \epsilon^{i} \;\;\mathrm{for} \; i = 1,2,...,m $$ where $r^i$ is a (T x 1) vector, $X$ is an (N x k) matrix, $\beta^{i}$ is a (k x 1) vector, and $\epsilon^{i}$ is $N(0, \sigma_{i}^{2})$ I assumed a g-prior for $\beta^{i}$ and a diffuse prior for $\sigma_{i}^{2}$. After deriving the posterior probability of each model, I get the highest- probability model (i.e. the set of best predictors). I then move on to testing these selected predictors using the **Cross- Sectional Approach** , also known as Fama-MacBeth procedure. This test reveals whether the selected predictors have a significant premium or not. The odd thing is that only one of the BMA selected predictors has a statistically significant premium. I am confused as to how to reconcile the results of BMA and Fama-MacBeth. Any help would be greatly appreciated. Thanks!


title: Sample balancing with batting averages

text: I'm trying to compare regional batting averages year over year (using comparable time periods i.e. August of each year) but I'm running into what I think is a sample imbalance. This is an extreme example, but say hypothetically a player goes 10-for-10 one day (double-header) and on the same relative day the next year, he goes 1-for-1. In both cases, his batting average for that day is 1.000. However, he obviously did much better the first day, but from an average perspective, there is no difference. Because of the regional data and time period assumptions I need to use, I will sometimes have 25 observations one day and 1 observation the same relative day the next year. I ran a Wilcoxin-Mann-Whitney test on the monthly samples comparing the two years (since the daily data seems to be too small) and they generally show no difference in the medians, but I'm looking for some opinions on if/how to balance the data by at-bats, rather then going strictly off of the batting average. Because of the collection method, I only have samples and at this point I am unable to gather additional data.


title: Activation Functions

text: In the linear regression models, the model prediction $y(x, w)$ is given by a linear function of the parameters $w$. In the simplest case, the model is also linear in the input variables and therefore takes the form $y(x) = w^T x + w_0$, so that $y$ is a real number. For classification problems, however, we wish to predict discrete class labels, or more generally posterior probabilities that lie in the range $(0, 1)$. To achieve this, we consider a generalization of this model in which we transform the linear function of w using a nonlinear function $f(\cdot)$ so that $$y(x) = f(w^T x + w_0).$$ Question: I am taking my first course in machine learning and am unable to wrap my head around the reason for the need to transform the linear function of $w$ into a non-linear function. And also could someone post a simple example of such a transformation.


title: Simple example that shows the advantages of Bayesian Model Averaging (BMA)

text: I'm incorporating a Bayesian Model Averaging (BMA) approach in my research and will soon give a presentation about my work to my colleagues. However, BMA isn't really that well-known in my field, so after presenting them with all the theory and before actually applying it to my problem, I want to present a simple, yet instructive example on why BMA works. I was thinking about a simple example with two models one can choose from, but the true data generating model (DGM) is somewhere in between and the evidence doesn't really favor any one of them. So if you choose one and continue from them on, you would ignore model uncertainty and make an error, but BMA, although the true model is not part of the model set, at least gives correct posterior density of the parameter of interest. For instance, there are two weather forecasts each day (A and B) and one wants to predict the weather best, so in classical statistics you would first try to find the best forecaster between the two, but what if the truth is somewhere in between (that is, sometimes A is right, sometimes B). But I couldn't formalize it. Something like that but I'm very open to ideas. I hope this question is specific enough! In the literature, I haven't found any nice examples from what I have read so far:   * Kruschke (2011), while a great introduction to Bayesian statistics, doesn't really focus on BMA and the coin toss example he has in chapter 4 is great for introducing Bayesian statistics, but doesn't really convince a fellow researcher to use BMA. ("Why again do I have three models, one saying the coin is fair and two saying it's biased in either direction?")   * All the other stuff I read (Koop 2003, Koop/Poirier/Tobias (2007), Hoeting et al. (1999) and tons of others) are great references, but I haven't found a simple toy example in them. But maybe I just missed a good source here. So does anyone have a good example he or she uses to introduce BMA? Maybe by even showing the likelihoods and posteriors because I think that would be quite instructive.


title: Chi square test for correlated parameters

text: I am doing a chi squared test with, in principle, 2 degrees of freedom. But it turns out that because my 2 parameters are correlated, my critical chi squared values are smaller than the canonical ones. Is there any test I can do to know _a priori_ the _effective_ number of degrees of freedom I should expect? Thanks.


title: Interpreting PCA scores

text: Just wondering if anyone can help me interpret PCA scores. So my data comes from a questionnaire on atititudes toward bears. According to the loadings, I have interpreted one of my principal components as "fear of bears". Then would the scores of that principal component be related to how each respondent measures up to that principal component (whether he/she scores positively/negatively on it)? Any help would be much appreciated! Thank you. A


title: Comparing sample profiles

text: I have some samples that have been treated with various chemical agents. Users score the treatment effect using categorical values/bins (i.e $< 1, 2.5, 8, ... > 100$). Within a range of say 1-100 they can get some fairly precise metric but outside of that range they input greater or less than some value. So the question is a user wants to know how similar one sample is to another... So I could do some sort of hierarchical clustering but I'm wondering with this mixed data type:   * What could be an alternative approach?    * What distance metric would be suitable?   * Are there alternative methods to clustering?


title: Why does the bootstrapped correlation revolve around zero while the original correlation $\approx 0.52$?

text: I have two data series containing 132 log-returns. One is for EURUSD, the other is for NZDUSD. The `head()` function shows you how some of the data looks. The correlation coefficient between the two, as calculated by `cor()` is $0.5178912$. To get a better sense of the correlation coefficient I bootstrap it by running `cor()` 1000 times on different 132 long samples. I run this in a loop and update `euro.nzd.corr` on every iteration. This is the R code I'm using:               head(euro)     [1] -0.001257862 -0.011637970  0.002428757  0.003602590 -0.003457319 -0.002012728     head(nzd)     [1]  0.008773255 -0.007744927  0.005498693  0.005642524 -0.000896363  0.003449576     cor(euro,nzd)     [1] 0.5178912     euro.nzd.corr <- numeric(1000)     for(i in 1:1000){     euro.nzd.corr[i] = cor(euro[sample(132,132,replace=TRUE)],nzd[sample(132,132,replace=TRUE)])     }     plot(density(euro.nzd.corr), lwd=3, col="steelblue")      Once I have the data, I plot the density chart, and get this: ![density](http://i.imgur.com/7IRIENi.png) Bootstrapped data has mean $\approx 0$ and mostly spreads between $-0.3$ and $0.3$. Where has the initial `cor()` result of $0.5178912$ gone? What am I to make of this? That it is better to conclude the two variables are uncorrelated versus correlated with a coefficient of $\approx 0.52$? Have I made any coding mistakes, or is the applied methodology simply flawed?


title: Panel regression: what to do when Hausman test fails and want to keep time invariant regressors?

text: I am running a panel data regression. First, I did a pooled OLS regression. Then I did a random effects (re) one. I carried out the Hausman test, and it refuted the null hypothesis (ie. I am discouraged to use random effects over fixed effects). So, I did the following: (1) I carried out a Hausman-Taylor regression (in Stata, `xthtaylor`). Using `xtoverid`, I got that this model is OK when compared with the fixed effects one. However, this model displays no $\text{R}^2$ value. Furthermore, I was told that Hausman-Taylor is not a good model to use when you goal is to use the model to estimate outcomes. I am trying to find out about the minimum distance estimator. It should be a procedure that allows me to combine fixed effects with time invariant regressors. Is there a Stata command for that? Do you know any reference on this topic? Any suggestion on how to deal with my current problem is very welcome!


title: illustrate interaction without ANOVA, other methods?

text: I am wondering whether there are acceptable ways to illustrate the presence of an interaction besides using an ANOVA. I can't use an ANOVA because I have unequal sample sizes, and due to the nature of the topic being investigated, there was no experimental design we could have used to avoid this.   Is using a scatterplot and adding R2 lines for both subgroups an acceptable method for demonstrating an interaction? If not, what are other methods? We replicated the effect twice (unequal sample sizes in both experiments) so we are confident the interaction exists. The only problem is demonstrating this in an acceptable-for-publication manner. Thank you in advance, all help is tremendously appreciated.


title: Creating a recommender system for shoe sizing

text: (Sorry if I have mis-categorized or titled this question.) Suppose you have shoe size data for a bunch of people across a variety of shoe brands. So, one person might have Nike = 11, Reebok = 10.5, Adidas = 11, Converse = 11.5, and so on. You might define a matrix $M_{ij}$ where $i$ represents the $i$-th individual and $j$ their size in the $j$-th brand. Then suppose you have a person who knows their Nike size is, say, 9, and wants to know their most likely Adidas size (or maybe they know their size in multiple brands, but not Adidas). What is the best way to predict their Adidas size? A simple method might be to look at the size difference of each person between Nike and Adidas ($M_{i,j=Nike} - M_{i,j=Adidas}$ for each $i$), and then to take the most common difference and add it to the person's Nike size to get their Adidas size. For example, if most people size up one half size when going from Nike to Adidas, the most likely Adidas size for this person would be 9.5. But is there a more sophisticated method? For example, maybe there are people who are outliers and have really weird feet, and the relationship of their sizes among all brands does not fit well with most other people. Is there some natural way to weight the data of these people less than that of others, and would you want to do this? I have a feeling that this is a common problem, for example anything online these days that makes a recommendation based on consumer data.


title: The product distribution: how fast does dissimilarity increase as a function of number of samples?

text: If $\mathcal{D}$ is a distribution, let $\mathcal{D}^n$ denote the $n$-fold Cartesian product of $\mathcal{D}$. In other words, $\mathcal{D}^n$ is the distribution of $n$-tuples $(x_1,\dots,x_n)$ where each $x_i$ is drawn iid from $\mathcal{D}$. If $\mathcal{D},\mathcal{D}'$ are two distributions, let $d(\mathcal{D},\mathcal{D}')$ denote a dissimilarity measure; it measures how different the two distributions are. For instance, we might use the $\ell_1$ distance measure (i.e., $||\mathcal{D}-\mathcal{D}'||_1$, aka total variation distance), the $\ell_2$ distance (i.e., $||\mathcal{D}-\mathcal{D}'||_2$), the $\ell_\infty$ distance (i.e., $||\mathcal{D}-\mathcal{D}'||_\infty$), the KL divergence ($D_\textrm{KL}(\mathcal{D} || \mathcal{D}')$), or something else entirely. Often we might expect these to be normalized so that the dissimilarity is between 0 and 1, but I suppose that's not necessary. Let $\mathcal{D}_0,\mathcal{D}_1$ be two distributions, and $d_1,d_2$ any two dissimilarity measures (maybe the same one, maybe different). Are there any bounds, estimates, approximations, or rules of thumb that relate $d_2(\mathcal{D}_0^n,\mathcal{D}_1^n)$ to $d_1(\mathcal{D}_0,\mathcal{D}_1)$?


title: training approaches for highly-imbalanced data set

text: I have a highly-imbalanced test data set. The positive set consists of 100 cases while the negative set consists of 1500 cases. On the training side, I have a larger candidate pool: the positive training set has 1200 cases and the negative training set has 12000 cases. For this kind of scenario, I have several choices: 1) Using weighted SVM for the whole training set (P: 1200, N: 12000) 2) Using SVM based on the sampled training set (P:1200, N :1200), the 1200 negative cases are sampled from 12000 cases. Is there any theoretical guidance on deciding which approach is better? Since the test data set is highly imbalanced, should I use the imbalanced training set as well?


title: probability by permutation

text: I tried to solve this question. My solution is below. Please could you tell me if it is correct or not? If 3 books are picked at random from a shelf containing 5 novels, 3 books of poems, and a dictionary, what is the probability that: (a) the dictionary is selected?   (b) 2 novels and 1 book of poems are selected? Answers: (a) $$\frac{\binom{1}{1}}{\binom{9}{3}} +\frac{\binom{8}{2}}{\binom{9}{3}} $$ (b) $$\frac{\binom{5}{2}}{\binom{9}{3}} + \frac{\binom{3}{1}}{\binom{9}{3}}$$


title: R lrm model with no predictors

text: I am creating multiple logistic regression models using lrm from Harrell's Design package in R. One model I would like to make is the model with no predictors. For example, I want to predict a constant c such that:               logit(Y) ~ c      I know I how to compute c (divide the number of "1"s by the total), what I would like is to use `lrm` so I can manipulate it as a model in a consistent way with the other models I am making. Is this possible, and if so how? I have tried so far:               library(Design)     data(mtcars)     lrm(am ~ 1, data=mtcars)      which gives the error:               Error in dimnames(stats) <- list(names(cof), c("Coef", "S.E.", "Wald Z",  :         length of 'dimnames' [1] not equal to array extent      and I have tried:               lrm(am ~ ., data=mtcars)      But this uses all the predictors, rather then none of the predictors.


title: What is AIC? Looking for a formal but intuitive answer

text: I've heard that AIC can be used to choose among several models (which regressor to use). But i would like to understand formally what it is in a kind of "advanced undergraduated" level,which i think would be something formal but with intuition arising from the formula. And is it possible to implement AIC in stata with complex survey data? Thank you very much in advance!! Regards!


title: Prediction using multiple training sets

text: I have multiple different training sets($TS_1$,$TS_2$,..,$TS_n$) and one test set $TS$.   I have calculated the prediction measures precision, recall, and F-measure for each pair ($TS_n$,$TS$). Is there any statistical method that can combine the prediction measures calculated for the different training sets?


title: Alternatives to MAD to find a yardstick to assess data

text: In a paper by Rousseeuw and Croux from 1993 ("Alternatives to the Median Absolute deviation", page 1274, link to pdf), I came across an indicator I'm considering using. The formula is:   $$ Sn = C\, {\rm median}_i \\{{\rm median}_j |x_i-x_j|\\} $$ The way I read the formula is that for every $X_i$ I take the differences from each of the other values in the sample. Then, again for every $X_i$, I take the median of the values I get. Finally I take the median of the median values I obtained in the previous step and multiply it by the parameter $C$ which according to the authors is $$ C= 1.1926. $$ I'm not a statistician by any means and I don't understand what the authors mean when they say that the value $C$ is derived through an asymptotic assumption. What I know is that my data is not normally distributed and resists normalization if I try to standardize it. This is why I'm trying to find an indicator where I can rank the data around a central point that does not demand symmetry in the distribution (hence, not using the MAD). **My question is: do I have to maintain the $C$ factor in the equation?** What is the underlying assumption? I seem to understand that it involves some kind of Gaussian distribution (what I know as a normal distribution) but since my data is not normally distributed, I am unsure what value for $C$ I should apply and if it's OK to just leave it out. **My end goal is the following:   1. I have data not distributed in any meaningful usable way. Yes, there may be ways to transform or fit it in some kind of distribution but that would really complicate things, beyond my ability to explain it to my colleagues.   2. Each data point is a score, ranging from 0 to several thousand, which is a function of the answers given in a survey.   3. I need to assess each of the data points in relation to one another, to determine which stand at the bottom of the risk scale (lowest scores) and which at the top (highest scores). To do so I need to define some cutoffs. My plan was to use the MAD or better, the Sn indicator as of above, to have an indication of what the score is for a common risk profile and then use percentiles to associate a lower risk profile to lower scores and higher risk profile to higher scores. I hope this is not confusing but I'm venturing into really uncharted territory for me here.


title: p-value for correlation/co-occurrence in sets of intervals

text: I'd like to have a probability measure for the **co-occurrence** of intervals in two datasets of intervals.               dataset1: [125,500],[900,1300],[2220,2500], ...     dataset2: [600,800],[1200,1400],[3020,3500], ...      (Two exemplary datasets with an overlap in the second intervals) The intervals within a dataset are represented by the coordinates of their start and end (just two integer values). Within a dataset intervals are non- overlapping. The number, average lengths and distribution of intervals may be very different between two datasets, but the length of the dataset that is the maximal integer a coordinate can have is identical and known for both. I am looking for a kind of p-value for the probability that there is a "correlation" between two datasets, that is that intervals overlap more frequently (or on the whole more widely in terms of intersection length) than can be expected by chance. I would be grateful for suggestions on how to address the problem in practical as well as statistical terms (perhaps using R). Perhaps there is even a solution for the "correlation" between multiple datasets ... Thanks a lot :-)


title: Should I check Multicollinearity When There is An Inverse?

text: At Machine Learning algorithms there are usually inversion process about matrices and sometimes Matlab throws error when Multicollinearity occurs. **Should I check Multicollinearity(and how) everytime before I inverse a Matrix (and use pseudo inverse - pinv)**


title: Circular statistics and bidirectional data

text: I am curious if tests such as Rayleigh's test, Kuiper's test, and Watson's test are valid for bidirectional data (i.e. 180 degrees) as well as unidirectional data.   If not, what are the appropriate corresponding tests?


title: Missing value replacement when 3 out of four scale items are missing

text: I have a subscale with four items and 316 participants. Five participants missed all four items, two participants missed three of the four items and five participants missed only one item (three participants missed the same item, the other two missed different items). The missingness on the subscale is MCAR. What would be the best approach for handling missingness (in SPSS)? Thanks very much, Tim


title: Model for Minimizing or Preventing Jet Lag

text: I've been relying on a relatively simple model for preventing jet lag, and was looking for better input/knowledge for more accurate modeling. So here goes mine so far: Our goal is to find the ideal number of hours to sleep on the plane, $S$, to minimize jet lag. To start, I assume, very simplistically, that the human sleep and wake cycle is like a charge cycle, with linear charge and discharge of energy. This means that, for an average of 8 hours slept in one 24 hour cycle, 1 hour of sleep $S$ = 2 hours of awake "energy" $E$. Now let's suppose we want to fly from San Francisco (SFO) to Shanghai (PVG). One possible nonstop flight is from 13:00 SF time to about 19:00 Shanghai time, with a flight duration of 14 hours. Further, let's suppose we have a regular sleeping cycle in SF of, say, 0:00 - 8:00, and we want to immediately start a normal sleeping cycle of 0:00 - 8:00 in China as well. From this point, it's a simple algebraic calculation. When we wake up at 8am in SF, we will have "energy" $E=16hr$. We want to end with $E=0hr$ at Shanghai midnight after the flight. The number of hours we stay awake before leaving SF and after landing in Shanghai but before sleeping is $13:00-8:00+24:00-19:00=10hr$. This leaves $E=6hr$ outside of the flight. On the plane, the number of hours slept $S$ plus the number of hours awake $A$ equals the total duration of the flight, 14 hours. It then follows that: $\begin{equation} \begin{cases} 6+2S-A=0\\\ S+A=14 \end{cases} \end{equation} $ Solve to get: $S=8/3hr\approx 2.7hr$ In general, the equations to solve for $S$ are: $\begin{equation} \begin{cases} E_{init}-(t_{start}-t_{wake}+t_{sleep}-t_{end})+2S-A=0\\\ S+A=t_{flight} \end{cases} \end{equation} $ where $E_{init}$ is the initial waking "energy", $t_{start}$ is the local flight start time, $t_{end}$ is the destination flight end time, $t_{wake}$ is the local waking time, $t_{sleep}$ is the destination sleep time, and $t_{flight}$ is the duration of the flight. (Using this model on another example from JFK to LAX, a 6.5 hour flight, I get $S=1hr$) I am looking to improve this model's accuracy in reflecting human sleep. I have considered perhaps modeling the sleep "charge" cycle as a sigmoid function, but that doesn't seem to be really plausible either. Perhaps considering how REM sleep cycles work, this model can be improved. Also, the solution does not take into account _when_ to sleep on the plane and if to sleep in intervals or a single block, so the answer could also use refining. From the equations, it should also be noted that $S$ is essentially independent of the flight start/end times as long as flight duration is the same. Interestingly enough, this model seems to work well for me on such long- distance flights, but only when going in the westward direction. If I use this model for eastward travel, it fails horribly. Any possible reason for this?


title: Network Visualisation for huge dataset composed of many unconnected clusters

text: I have a dataset which consists of 190,455 nodes and 1,241,638 edges. These can be broken down into a set of 2300 subgraphs which are not connected to each other. I am having trouble creating visualisations for this using tools such as Gephi (Yifan-Yu seems to work best) as the layout takes a very long time and doesn't seem to separate the subgraphs well. Is there a tool which performs a layout for such a dataset well? It would seem best to plot each subgraph individually, and then perform some merging at the end to avoid overlap while the algorithms seem to be trying to plot and arrange everything at once. Thanks!


title: Determining conserved features using a Bayesian approach

text: I would like to perform some sort of binary classification, and my data set consists of 100 examples (for each class), which are vectors with 2500 elements. Ideally, I would like to determine which the most important features are (i.e. best discriminating between those two classes) and reduce the number of features to a more reasonable dimension, e.g. 50. I decided to start off with a naïve approach, using the Standard score (by only determining the indices of the most prominent features, but without modifying the actual values), which obviously does not take into account the distribution of values per each class. What I basically do is I calculate the score and order in descending order of prominence, and store the indices. Then I just keep the first $n$ (say 50, as suggested above) of them and extract only those features. Just recently, it was recommended to me to take a look at the following paper: Bayesian automatic relevance determination algorithms for classifying gene expression data, and try to estimate the probability of common feature occurring together over a given number of iterations (i.e. in order to estimate an optimal $n$, as opposed to picking one arbitrarily). Ignoring the presented two algorithms in the paper, I focused on the following: $\mathbf{P}(N_c=n_c|n,n_1,n_2) = \frac{\binom{n_a}{n_c}\binom{n - n_a}{n_b - n_c}}{\binom{n}{n_b}}$, where $n_a = \max(n_1, n_2)$, and $n_b = \min(n_1, n_2)$ So basically, my approach is to perform the following over a number of iterations: >   1. Put together the training sets of the two classes (and keep the indices > of the examples of each one of them) >   2. Split the resulting training set in half and keep working separately on > each half (subset) >   3. Determine the m (which is not necessarily equal to n) most prominent > features using the Standard score, and keep track of the frequency of > occurrence of each > At the end, I calculate $N_c$ for $c$ from 1 to my feature vector length (i.e. 2500) and thus I get an optimal number for $n$ (the reduced number of features, which if picked arbitrarily could lead to missing out important features, or including not so prominent ones). So here comes my confusion - is my approach so far correct at all? And if yes - what do I do with the $N_c$ which has the highest probability. Do I put together the $N_c$ most prominent features (from each subset), based on the accumulated values of the frequency vectors, or I should do something entirely different (e.g. like taking the best $N_c$ features from a single scoring - using the Standard score - over the entire training set)?


title: gaussian mixture model - approximate a matrix

text: I have a similarity matrix M - the value M(i,j) indicates the similarity between two elements i and j. I want to approximate that matrix using a Gaussian Mixture model or I want to cluster that matrix into set of K clusters. How can I do it? I am not sure about the input of the clustering algorithm.


title: Interpreting the results of post-hoc ANOVA

text: After concluding from a one-way ANOVA, that a variable (country) is indeed statistically significant, I conducted Tukey's multiple comparison post hoc test. Here is the plot for mean Cost (DV) vs Country with the bars representing Confidence Intervals: ![enter image description here](http://i.stack.imgur.com/PLKAP.png) Now, since there are multiple countries being tested here (even more countries in the actual data), it becomes very complex trying to manually analyse which ones differ. What I could make out of it is that certain countries could be considered as 'clustered' together (eg. diagram below). Then it naturally follows that tukey's test be performed for _within_ a cluster, and obviously be performed again _across_ clusters. Is that the way to go about this, creating clusters ? How to account for effect sizes ? Too many groups actually blur any conclusions. Here is another plot with more countries: ![enter image description here](http://i.stack.imgur.com/ZYafj.png) Many Thanks!


title: Trying to understand unbiased estimator

text: After reading this post, I still don't thoroughly understand what an estimator is. Suppose the samples $D_i={(x_1,y_1),...(x_n,y_n)}$ are drawn randomly from function $$f(x)=sin(2\pi x),$$ so my ultimate goal is to come up with an estimated function $h(x)$, which should approximate $f$ as close as possible, right? Since I don't know what the true $f$ looks like, so I might choose different models to do the approximation, here I consider the linear regression model, which is $$h(x)=\beta^Tx,$$ then I try to estimate the model parameters $\beta$ via OLS, and finally given training sample $D_i$, I have the estimates $\hat\beta_{|D_i}$, and absolutely dffierent $\hat\beta$ for different training set $D_i$. Here is my problem,   1. What is the estimator here? According to the post I read, I think that the process of _choosing the linear regression model and estimating $\hat\beta$ via OLS_ are components of the estimator, right?   2. How to check the estimator is biased or not? Since I chose the linear regression model $h(x)$, and now have the estimates $\hat\beta$, I think to check the estimator is biased or not, is to find out $$B(\hat\beta)=\beta_{true}-E[\hat\beta],$$ right? Now this is the thing confuses me badly, since we don't have this $\beta_{true}$, because the true $f$ is actually $sin(\theta x)$ with $\theta=2\pi$, we only have $\theta_{true}$, no such $\beta_{true}$, right? Hope you can help me to get things right. **UPDATE** To address @whuber's comment, my question originates from this lecture note (from page 11) 1. _how do you "randomly draw" values from a function?_ Let's say the target function $f(x)=sin(2\pi x)$ is defined on $x\in [-1,1]$, I could randomly (via a uniform distribution over $[-1,1]$) draw $N$ points $\\{x_1,...,x_N\\}$, where $N$ is the size of the sample, after that I could compute the output value $y_i=f(x_i)$, no noise added. Now the sample is $\\{(x_1,y_1),...,(x_N,y_N)\\}$. 2. _are either (or both) of the $xi$ or $yi$ measured? what is the expected statistical behavior of the measurement errors?_ Here I only measure the training/test error of output value $y_i$. I mean, with different training sample, I have different $\hat\beta$, then for a fixed point $x_0$, the predicted output $\hat y_0=h(x_0)=\beta^Tx_0$ will be different, but the true value is $y_0=f(x_0)$. 3. _why are you choosing a linear estimator β′x to approximate a function sin(2πx) that is a fortiori strongly nonlinear?_ As what is done in the lecture note, what I have is just the sample, which looks nonlinear, maybe it's not wise to choose a linear estimator, and maybe I should choose this one $$h(x)=\beta_0+\beta_1x+\beta_2x^2$$, a polynomial estimator, but it's still far from the target $sin(2\pi x)$.


title: How to achieve strictly positive forecasts?

text: I am working on a time series whose values are **strictly positive**. Working with various models including AR, MA, ARMA, etc, I couldn't find an easy way to achieve strictly positive forecasts. I'm using **R** for doing my forecasts, and all that I could find was forecast.hts {hts} that has a **positive** parameter described here: **Forecast a hierarchical or grouped time series, package hts**               ## S3 method for class 'gts':     forecast((object, h,       method = c("comb", "bu", "mo", "tdgsf", "tdgsa", "tdfp", "all"),       fmethod = c("ets", "rw", "arima"), level, positive = FALSE,         xreg = NULL, newxreg = NULL, ...))          positive         If TRUE, forecasts are forced to be strictly positive      http://www.inside-r.org/packages/cran/hts/docs/forecast.gts Any suggestions for non-hierarchical time series? What about generalization on using other constraints like minimum, maximum, etc? Even if not implemented in R, suggestions on articles, models or helpful general variable transformations would be appreciated.


title: Computing a multi-sample (i.e., pooled) Akaike Information Criterion

text: I have a set of multivariate time series observations that I am trying to model using VAR processes, using AIC to choose the best model. However, instead of determining the best model order for each individual sample (itself multivariate), I would like to figure out the best model order for fitting _all_ the samples. However, each sample will still have unique parameters. In other words, how can multiple AICs be combined? More information: > All criteria [AIC, HQC, SWC] add a penalty to the one-step ahead MSE which > depends on the sample size $T$, the number of variables $m$ and the number > of lags $q$. (F. Canova, _Methods for Applied Macroeconomic Research_ ) I have $N$ samples, all with the same number of variables $m$, all cropped to identical length $T$. I would like to compute the pool information criterion for some lag $q$. I can easily sum or average the one-step ahead MSEs for each sample. But what is the best way to pool the penalty terms? The naive approaches would be to either sum (/average) the penalty terms or just apply a single penalty term. The latter seems to capture the situation where you are fitting the same VAR model to each sample. Does the former capture the situation where you are allowing the parameters to vary for each sample but keeping the number of lags $q$ constant?


title: How to calculate power of different normality tests such as Shapiro-Wilk, Ryan test etc

text: Normality test for small samples with power comparison


title: Inconsistent delta.chisq.scaled using Lavaan and semTools?

text: I've recently been using Lavaan and semTools to test for measurement invariance in confirmatory factor analysis (CFA) models and I noticed an (apparent) inconsistency in a calculation that confused me. When you use the "measurementInvariance" command from semTools to test nested models (using MLR), the chi-square difference value (delta.chisq.scaled) is not the same as when you manually extract the fit statistic from each model and subtract them yourself. It also does not match the difference in the non-scaled chi-square values. Am I missing something about how this is calculated? I've given an example below using demo data from Lavaan. Thanks in-advance!               require(lavaan)     require(semTools)     HW.model <-       'visual =~ x1 + x2 + x3      textual =~ x4 + x5 + x6      speed =~ x7 + x8 + x9      '     out<-measurementInvariance(HW.model, data=HolzingerSwineford1939, estimator = "MLR",group="school")     modelDiff<-compareFit(out)     summary(modelDiff, fit.measures="all")      The resulting values are:   1. measurementinvariance command chi-square (scaled) diff = 6.567   2. manual subtraction of extracted chi-square (scaled) values = 4.25   3. manual subtraction of extracted non-scaled chi-square values = 7.68 Does anyone know why the values from list item 1 & 2 above differ?


title: Simple, reliable, open, and interoperable plain text format for storing data

text: In a previous question I asked about tools for editing CSV files. Gavin linked to a comment on R Help by Duncan Murdoch suggesting that Data Interchange Format is a more reliable way to store data than CSV. For some applications a dedicated database management system is what is needed. However, for small scale data analysis projects something more light weight seems more suitable. Consider the following criteria for evaluating a file format:   * **reliabile** : the data entered should stay true to what has been entered; data should open consistently in different software;    * **simple** : it would be nice if the file format is easy to understand and ideally be readable with a simple text editor; it should be easy to write a simple program to read and write the format.   * **open** : the format should be open   * **interoperable** : the file format should be supported by many systems I find tab and comma separated value formats fail on the reliability criterion. Although I suppose I could blame the importing and exporting programs rather than the file format. I often find myself having to make little adjustments to the options in `read.table` in order to prevent some strange character from breaking the loading of the data frame. Questions   * Which file format best meets these needs?   * Is Data Interchange Format a better alternative? or does it have its own problems?   * Is there some other format that is preferable?   * Am I unfairly evaluating TSV and CSV? Is there a simple set of tips for working with such files that make the file format more reliable?


title: Export/format output from pairwise.t.test to LaTeX

text: I am trying to format the output from pairwise.t.test into LaTeX, but have not found a way of doing this. Has anyone got any suggestions? **EDIT:** As this is a one-time only report where I do need to customize the variable names, and row-/column headings, I was hoping to avoid using Sweave.


title: How to calculate the probability for pattern that violates "control chart rules" to occur?

text: I am going to host a training session to teach healthcare staff how to use control chart (c-chart to be specific), and I need to tell that why these rules (Western Electric Rules) are so called rules that when pattern is matched, outbreak can be considered as a probability. I have read the article "Control Charts 101: A Guide to Health Care Applications" a bit and it teaches me a bit, saying the chance for one data point to exceed UCL is less than 0.5%, I can calculate that using the following R code without problem:               1-pnorm(3, mean=0, sd=1)     [1] 0.001349898      I have written the following R code to verify the same thing in c-chart (i.e. poisson rather than normal distribution), for lambda ranged from 1.0 to 100.0               UCL <- integer(0)          for (loop.UCL in (1:1000))         {lambda <- loop.UCL/10          result <- 1-ppois(lambda+ceiling(sqrt(lambda)*3),lambda)          UCL <- c(UCL, result)         }          summary(UCL)          Min.   1st Qu.    Median      Mean   3rd Qu.      0.0009679 0.0016160 0.0018260 0.0019960 0.0021220           Max.      0.0134600       I think I can handle the first Western Electric Rules, but I am puzzled for the second and third. Can anyone give me some insight or hint? Thanks!


title: Learning to predict maximum of parameterized function class

text: I am interested in a multi-task regression problem: I have a parametrized function $f_x : \mathcal{R}^n -> \mathcal{R}$ where $x \in \mathcal{R}$ is a real-valued parameter. For some values of $x$, I have a large set of (noise- less) training examples. I would like to learn to predict for some new $x'$ where the maximum of $f_{x'}$ lies. $f_x$ and $f_{x'}$ should be similar if $x'$ is close to $x$. The $f_x$ are not restricted to any parametric function class but are upper-bounded by $0$. As an add-on, it would also be very nice if not only the estimate of the maximum but also some uncertainty-measure would be available. Think of the estimated maximum of $f_{x'}$ as a starting point for a black-box optimization of $f_{x'}$. My first idea was to treat the entire problem as a single regression problem $f: \mathcal{R} \times \mathcal{R}^n -> \mathcal{R}$. However, this is not very amenable for estimating the maximum of $f_{x'}$. An other option would be to learn a mapping $g : \mathcal{R} -> \mathcal{R}^n$, which maps task parameter to the estimated maximum of $f_{x'}$, i.e., $g(x') = \arg\max\limits_{y \in \mathcal{R}^n} f_{x'}(y)$. This could be trained based on the pairs $(x, \arg\max\limits_{y \in \mathcal{R}^n} f_{x}(y))$ for the $x$ for which training data exists. However, it would loose all information contained in the non-maximal datapoints. Any suggestions?


title: Bayesian classifier and discovery of new classifications

text: I've written Naive Bayesian classifiers before, they work wonderfully. But I'd like a classifier which will learn like a Bayesian classifier and identify new classifications when a new cluster emerges. Lets presume I have sensor data for some sort of machinery. I have training data for different conditions - too hot, to cold, misfires, etc. Now lets assume some new condition occurs - say a misalignment, or say over revving. What techniques would allow for that new cluster of data to build a new classification? Naive Bayes+K-Means?


title: Trend in data - determining according to angle between regression line and vertical line?

text: This question is a little bit referring to this question How to determine trend strength from linear regression slope? but I found a another solution so I am creating new question to confirm my idea. In my previous question I need to determine how is a trend in dataset steep. I were trying to do it with slope parameter from linear regression equation. I have a another idea and I would like to consult it. According to linear regression equation I can determine two points that lie on this regression line. Next I can create a line that crosses one of these two points and is collinear with x-axis. Now I can compute degree between this horizontal line and regression line. Maximum value of this angle is 90°. So I can compute my angle (e.g. 48°) and divide it with 90°.               Trend strength = 48 / 90 = 0.53      what means 53% of its maximum value and with this I can work in my application. **Question:** Is my idea correct from statistical perspective?


title: Designing a test for a psychic who says he can influence dice rolls

text: Let's say that I have a friend (let's call him "George") who says that he can control the roll of dice using his mind (ie, make the dice more likely to fall on a specific number he is thinking of). How can I design a _scientifically rigorous_ test to determine if he can actually do this? (I don't really think he can, of course, but I want him to agree to the details of a test, Amazing Randi-style, before the test begins.) I want to reduce the (very likely) post test-excuses that he'll come up with. Here is what I have so far:   1. Determine the physical dice-rolling technique (which dice, shaker cup, landing surface, etc)   2. Define a "test session", consisting of **X** rolls of the dice. This needs to be small enough to do in one sitting, but large enough to determine (after analysis) within 95%-99% confidence whether the dice fell fair, or favored one side   3. Run **Y** sessions on the chosen dice (with no influence from George), as a "control" to make sure that the dice show "fair" results on their own   4. Run **Z** sessions with George. Before each, roll a separate die to determine what number George will "concentrate on" during that entire session.   5. Compile and analyze the results.   6. George makes up some excuses for his dismal performance. So my questions for you:   * Any flaws or problems with my overall methodology? Anything George would likely object to?   * Should I use a D6? Or a D20? Does it matter? Would a die with more faces require more rolls to produce similarly confident results? Or the opposite? I'd prefer fewer rolls than more, due to practical considerations :)   * What are reasonable values for **X** , **Y** , and **Z**? They are not entirely unrelated; if my chosen value of **X** only allows 95% confidence for a single session, then 1 out of every 20 sessions could "fail", even without George's influence   * How do I define "success" or "failure" for an individual session? (I did find this question which goes through the details of a chi-squared test, so I think that's my evaluation methodology, but what are reasonable confidence thresholds?)   * How do I define "success" or "failure" for the overall test? George might "win" a single session by sheer chance, but how many of the **Z** sessions would he have to pass to pass the whole test? I'll probably be analyzing these results in a MS Excel spreadsheet, if that makes any difference.


title: Extracting multiple columns from a matrix in R

text: If I have a matrix `M` of 15 columns, what is R syntax to extract a matrix `M1` consisting of 1,7,9,11,13 and 15 columns?


title: Help to analyze results of a multinomial logit

text: After I have followed the advice I have been given on that website, I have been able to get some results, unfortunately I have some trouble understanding what I got: original post : Here edit : Direct link to image of results Here So I need some clues on how to see if what I have is relevant. Also, to what does the 1, 2, 3 (e.g: month: **1** , month: **2** , month: **3** ) represent? Hope someone will be able to help me Thanks in advance! Edit: So I have followed a link from the original post and used the VGAM library instead of the mlogit since my data was not formated for it. I ran then the vglm function and ran a summary. I got, in the end, those results. I have red the presentation available on this website : http://www.stats24x7.com/3dMultinomialLogisticRegression.pdf which says (I quote): > Educ is not significant (indicated by small magnitudes of t‐values) – run > model without educ. So I was wondering what means the "small magnitudes", shall I get rid of Month and Age (i.e: does the differences between the number counts or the fact that it is near 0 ?) Also, how do I get the p-value ?


title: Robust regression - a better understanding

text: I looked at robust regression for the first time today and I am a bit confused, comparing it to something like ordinary least squares and I am not sure if I am on the right track. I read a few articles and they say that with robust regression you don't need to worry too much about outliers and heteroscedasticy and the normality restrictions on the residuals is not as important as with OLS. When I do robust regression in R however, I also don't get any significant indicators regarding coefficients (I am using the `rlm()` function). So when robust regression like M or MM-estimation is applied in a regression model, is the significance of coefficients still important or is the idea of robust regression just to fit the best possible plane through the data points and find the coefficients?


title: Why is there a difference in signs of regression coefficients of same variable in a simple regression and multiple regression on spss

text: I have checked the relationships of dimensions of independent variable with the dependent variable in a simple linear regression on SPSS, but when I performed multiple linear regression on SPSS of the same dimension collectively then some dimensions have changed signs of regression coefficients against simple linear regression. what is the reason for this change in relationship? e.g., in a simple linear regression analysis debt to assets ratio showed positive relationship with ROA (Return on Assets) but when I add more variables then debt to asset ratio shows negative relationship with ROA.


title: Kolmogorov-Smirnov CDF

text: I'm writing some C++ code that uses a two-sample K-S test. I've computed the largest difference between the two cumulative density functions, $D$, and then compute the test statistic $x = \sqrt{\frac{n_1 n_2}{n_1 + n_2}} D$. Does anyone know of a simple form for the cumulative density function under the null hypothesis? I found the asymptotic form when $n \rightarrow \infty$ on Wikipedia: ![](http://i.stack.imgur.com/l9kEk.png) But I haven't been able to find any approximation where $n$ is not large. Can you give me the formula or provide a reference? Thanks a lot.


title: Covariance when exponentiating a normal random variable

text: Say $Z \sim N(\mu, \sigma^2)$. I am trying to figure out what is the variance of $Y = exp(Z)$. The first thing that came to mind was approximate $exp(Z)$ with a second order taylor series, then it boils down to finding the variance of a polynomial. But then what is the covariance term of $(Z, Z^2)$?


title: How to get F-test p-value

text: Suppose we can choose from two different catalysers. 10 observations are taken from the first one and 12 from the other one. If $s_1 = 14$ and $s_2 = 28$, can we reject at $\alpha = 5\%$ the hypothesis that the variances are equal? Here is what the teacher did : The ratio is : $s_1/s_2 = 0.5.$ Then $$P(F_{n=9, m=11} \le 0.5) = 0.1538$$ Then he says : the p-value is $2 \times \min(0.1539; 0.8461) = 0.3074$ and he rejects $H_0$. **How do I get the 0.1538 ?** I think I can check a F-table for n=9, m=11 , but what do I do then to get the probability that this value is $\le 0.5$ ?


title: Independent variables in ordinal logistic regression

text: One of my IV's for my ordinal logistic regression is a nominal categorical variable with 4 categories. Most examples I see for this type of logistic regression have only binary categorical variables. I have tried with and without dummies for my 4-category independent... but I am not confident of the outcome. When I renumber the categories (switch them around as 1,2,3,4 because they are not ordered in any way, I get different parameters and p-values each time. So something is wrong.


title: Problem with Poisson distribution

text: Let X and Y be two independent real random variables following two Poisson distributions that are denoted respectively by $P(\lambda_1)$ and $P(\lambda_2)$. Let $Z = X + Y$. Now:   1. Determine the generating function of $Z$. Deduce the distribution of $Z$.   2. Find the conditional distribution of $X$ given $Z = n$. I'm concerned about the answer of part 2.


title: Constrain decision boundary to fall on grid lines in multiple class logistic regression

text: I would like to use multiple class logistic regression to learn the decision boundaries separating the different classes (denoted by color) in the image below. Kernel logistic regression with a RBF kernel seems like a good choice, but I would like the decision boundary, when projected back to the 2-d space, to fall along the white grid lines. One way of proceeding would be to introduce a penalty term in the objective function, but I'm not sure how to proceed. Does anyone have any references for using a grid as a contraint in this kind of problem? A perhaps different way of asking this is how can I combine, on the one hand, a prior belief of the form $\Pr(\int^\mathbf{x}\pi=.5 \mid \mathbf{x})$ about the distribution over $\mathbf{x}$ of when the cumulative probability distribution of $\pi$ takes on a particular value with, on the other hand, a likelihood $\Pr(\mathbf{x} \mid \pi)$? I ![Different classes of points on a grid](http://i.stack.imgur.com/BnFVP.png)


title: Selecting a model for predicting a multinomial response from multinomial data

text: I am trying to create a predictive model from a set of categorical observations:               task   user     ----   ----     S      Alice     M      Bob     M      Alice     S      Charlie     L      Bob     M      Charlie     S      Alice     ...    ...      I'm not interested in predicting what type of task will be next for a given user. Rather, I'm interested in predicting a user's most likely distribution of tasks. For example, if a user has a set of 20 tasks, how many of them will be S, M, and L? The actual dataset has several thousand distinct users that can be grouped hierarchically. Each user belongs to a workgroup, and each workgroup belongs to a department. It's pretty easy to come up with a contingency table and empirical distribution, and I can see that the distribution varies among departments, and among workgroups within each department. I believe that variations among individuals within the same workgroup are insignificant, but this is a statement from intuition, and I don't know how to justify it statistically. Currently, my predictive model just calculates the distribution of tasks across the entire dataset and uses the same prior distribution for all users: if S/M/L are 0.2/0.5/0.3, then the prediction for every individual is based on this distribution. But in practice, this model doesn't work as well as I would like for several workgroups (and one entire department) because their task distribution appears to be different from that of the whole company. I'm thinking that in these situations, my model needs to use a workgroup-specific or department-specific distribution for the individuals within it. But I would like to avoid over-fitting the model. So, a few questions:   1. How can I determine if the distribution in a given department or workgroup differs from the total distribution sufficiently to justify using a different set of parameters for that department or workgroup?   2. Let's say I generate a model in which each department has its own distribution. How do I determine whether this model is "better" than the model that uses the same distribution for all departments, or whether the department-specific model is overfitted?   3. Let's say I build a model in which the predictive distribution for an individual is a weighted mixture of the total, the department and the workgroup distributions. How can I figure out the appropriate weights for this mixture? My dataset is huge (>1M rows), and I have access to Python+SciPy and R. My programming knowledge is much better than my statistical knowledge, so if you describe the statistical technique I can figure out how to implement it in code. Thanks.


title: Modelling receipts, volume and price

text: I am trying to model receipts from a transaction tax, where the the tax rate varies by band according to the price. The basic equation is something like this: receipts = volume * price. I am using average prices and total volumes as the covariates. I am not sure if there are possible simultaneity issues here, because volume and price are correlated. Can anyone give me some advice on how to test/address this? Do I need instruments for price and volume?


title: ANOVA or other test for summed data

text: I would appreciate any suggestions of how to analyse the following dataset. It was suggested that I use ANOVA, but I wanted to check with the community. I have a protein where I have calculated the entropy of each reading frame (blue, red and yellow). In the attached figure, the bars are a sum of the values calculated. It would probably be better to represent the values as means rather than total sums. Anyway, eye-balling it, the yellow section has a higher entropy than the blue and red regions, whereas red and blue look similar. However, I need to add some statistical weight to this and I would only compare reading frames within one genotype. I'm unsure whether to use ANOVA (repeated as it's 3 groups ?), a pairwise.t.test in R, or even Friedman's test, it's been a while. Perhaps this question is similar to How to test hypothesis of no group differences? ? I didn't find questions and answers that I identified with, so apologies if I missed one. ![enter image description here](http://i.stack.imgur.com/ICIbD.jpg)


title: Bayesian estimators

text: I have a model $f(x|\theta)$ ($\theta$ is a vector) for which I want to specify a prior $\pi(\theta)$. I only know that $\theta$ is in some interval. There are ways to specify an ignorance prior (Jeffery's, transformation group, reference prior etc.). Not knowing how to choose, I have simulated $f(x|\theta=\theta_0)$ many times and tried to estimate with each prior candidate, so that I have a distribution $D(\hat{\theta_0}|prior)$. Based on these distributions, I could choose the "optimal prior" using some criterion (e.g. MSE). Is it reasonable at all to proceed this way?


title: Correlation significancy measure for multiple predictions

text: I'm not sure what statistical test to use in the following situation: I have three different algorithms producing regressions for an observation. I now want to know, which of the algorithm yields the best prediction for the true values, i.e. which algorithms output correlates best to the observation. Specifically: is one of the algorithms significantly better than the others?


title: Asymmetric Dynamic Conditional Correlation in Matlab

text: I need to interpret the results of the estimation of the A-DCC model by Cappiello, Engle and Sheppard, 2006, Journal of Financial Econometrics. I used the Matlab routines of the MFE Toolboox http://www.kevinsheppard.com/MFE_Toolbox in this form: parameters = dcc(data,[],1,1,1,[],[],[],[],'2-stage') where data is a 100x2 matrix, hence the estimation is performed on two variables observed 100 times. The output "parameters" is a vector of 10 components. Unfortunately I am not able to derive a precise expression of the model and hence to interpret the parameters. The help of the routine says: _OUTPUTS: PARAMETERS - Estimated parameters. Output depends on METHOD. 3-stage: [VOL(1) ... VOL(K) corr_vech(R)' vech(N)' alpha gamma beta] 2-stage: [VOL(1) ... VOL(K) corr_vech(R)' alpha gamma beta] where VOL(j) is a (1+P(i)+O(i)+Q(i)) vector containing the parameters from volatility model i._ Which is the volatility model i? Unfortunately the documentation on the A-DCC model in the MFE toolbox page is missing. Note that by defult P(i)=1, O(i)=0 and Q(i)=1. Something more is said in the end of the help:                 2-stage:       Q(t) = R.*scale + a(1)*e(t-1)'*e(t-1) + ... + a(m)*e(t-m)'*e(t-m)       + g(1)*v(t-1)'*v(t-1) + ... + g(l)*v(t-l)*v(t-l) + b(1)*Q(t-1) + ... + b(n)*Q(t-1)      where v(t,:) = e(t,;).*(e(t,:)<0) and s = sqrt((1-sum(a)-sum(b)-gScale*sum(g))) and scale = s*s' Note that, having specified 1,1,1 in the call of the command, only one auto- regressive component is present for e(t), v(t) and Q(t). Q(t) should be a 2x2 matrix and hence a(1), b(1) and g(1) should be 2x2 matrices. My problem is to find how the element of the three matrices are linked to the entries in the vector PARAMETERS.


title: Testing for confounding

text: Suppose we want to test whether $Z$ is a confounding variable for the effect of $X$ on $Y$. Is it enough to just check unadjusted and adjusted estimates of the coefficient of $X$ and see if they differ? Also, if the difference between the two is very small (e.g $1/1000$) can we say that there is no confounding? It seems in survival analysis it would be complicated to check for the association of $Z$ vs. $X$ and $Z$ vs. $Y$.


title: How to build a prediction model for exam score based on previous scores

text: I am trying to construct a formula, which will take student's previous exam results (for ex: SAT) taken at particular dates and predict his future test result. One X is previous test result 1; another X is date of previous test 1 (can be converted to number of days between this test and the last test for simplicity); and other Xs are similar variables for additional previous tests, of which there are 3-5 per person. My Y is the result of the last test. Normally I would use simple linear regression to model this relation, but the problem is that this relation is not linear, because improving one's score from 100 to 200 is easier than from 300 to 400, for example. And also because of the upper limit of the test score (700 for example). Is there a way to create a more or less meaningful model for such prediction given 3-5 previous test results? Thank you!


title: Getting NaN value for X-Squared and N/A for p-value in Chi-Square Test

text: I have to perform a chi-square test on this given data set: ![http://img27.imageshack.us/img27/2600/dataset.jpg](http://i.stack.imgur.com/ilyCw.jpg) Upon doing so with the following code:   ![http://img199.imageshack.us/img199/4837/chisqcode.jpg](http://i.stack.imgur.com/Wr5iA.jpg) I get the following result:                       Pearson's Chi-squared test          data:  data.table      X-squared = NaN, df = 44, p-value = NA          Warning message:     In chisq.test(data.table) : Chi-squared approximation may be incorrect      Why is the X-squared value "NaN" and the p-value "N/A"? What am I doing wrong here?


title: VAR and Granger causality test

text: Is it necessary to calculate VAR before Granger causality test so that we can have the lag length to be used in Granger causality test


title: Robust estimation in SPSS generalized mixed models

text: I'm using mixed models in SPSS 19 to analyse dietary data. The mixed procedure is used because we have more than one measurement from many of the participants. My problem is that many of my dependent variables have a high number of zeroes. This does not mean missing data, but that the person hasn't consumed the food item in question. I've seen that there is an option called "Use robust estimation to handle violations of model assumptions" under the Generalized Mixed Models in SPSS. Does anyone know if this would be helpful for me? I have a sample size of about 1800 participants. Using the robust estimates does not seem to influence my results much. Thanks! :)


title: Absolute risk and prevalence?

text: From Wikipedia > If the absolute risk in the control group is available, conversion between > the two is calculated by: $$ RR \approx \frac{OR}{1 - R_C + (R_C \times > OR)}$$ where: > > RR = relative risk > > OR = odds ratio > > RC = absolute risk in the unexposed group, given as a fraction (for example: > fill in 10% risk as 0.1) What is the difference between absolute risk and (point) prevalence (used in rare disease assumption)? Thanks!


title: Why PCA performed on two similar data sets result in different number of components?

text: I have two data sets (2048 dimensions) collected under slightly different circumstances. I am using PCA to reduce the dimension of the data before passing it further for classification. Both data sets are spectra of the same alloy samples collected with slightly different lasers. My problem is that of the two sets of _visually_ similar data one requires 60 components and the other 275 given the same amount of variance retention (99.9%). I would like to ask is there some technique that would allow me to identify the source (features or combinations of thereof) of this behavior? Is there some standard method other then simply comparing the total variance?


title: Spearman Rho or ANOVA when examining effect of ordinal IV on continuous DV?

text: I have an ordinal independent variable (<1hr, 1hr, 1-3hrs, 3-5hrs, >5hrs) and an interval dependent variable (test scores). Apparently I am expected to examine the relationship between the two by conducting Spearman Rho correlations. However, would it not be more useful to perform a one-way ANOVA and post hoc test? Especially since the mean scores increase across ranks of the IV until the highest rank, at which point they drop. Or will reporting results from both tests be useful in data analysis?


title: Measure for the uniformity of a distribution

text: I can't seem to find a well established and simple statistical measure of uniformity in occurrence datasets in the presence of zero-valued categories. I've looked at Shannon's entropy which seems to be the closest but the problem is that it can't tell me anything about variables which have no occurrences in them. I always have a set number of variables e.g. 5 with each consisting of 0-10 occurrences. Although entropy can tell me about the distribution across all variables if they all have more than zero occurrences, in the presence of zero it becomes somewhat meaningless. Had a look at How does one measure the non-uniformity of a distribution? but I believe this is a different case. In my case I have a minimum and maximum number of occurrences that can occur in each category, (0-10) and a fixed number of these categories (5). Can anyone point me in the right direction? Example;               10 10 10 10 10      H(X) = 2.32193 and metric entropy = 0.04644      Ideally a metric would give me an extreme value of 1 for uniformity               10  0  0  0  0      H(X) = 0 and metric entropy = 0      Ideally a metric would give me an extreme value (or near it) of 0 for non- uniformity               10 5  0  0  0     H(X) = 0.9183 and metric entropy = 0.06122      Many thanks and sorry for my lack of statistical knowledge


title: Which test can I use to compare two sets of categorical data?

text: I am really stuck with the statistics section for my dissertation. I collected data as to whether items were arriving on a hospital ward; it was YES/NO, I made a change to the system and re-measured (different sample) using the same YES/NO answer, as there is no other categorisation of the data I am unsure which tests I can use, I understand this is really basic compared to what people are talking about on this site but I am sooooooo stuck (head in hands). The data is in SPSS.


title: In R, why are residuals of lm() and summary(lm()) different when using weights for the regression?

text: I am using a decay weighting scheme for time series regression in financial markets (idea is that more recent data is more relevant). However `lm()` residuals are different from `summary(lm())` residuals. Is there a reason for this?               rr <- rnorm(500) * 50     xx <- 1:500     yy <- xx + rr     dd <- decay(500, 150)     sd(rr)     [1] 50.31787     sd(lm(yy ~ xx)$residuals)         [1] 50.29099         sd(summary(lm(yy ~ xx))$residuals)     [1] 50.29099     sd(lm(yy ~ xx, weights = dd)$residuals)         [1] 50.29154         sd(summary(lm(yy ~ xx, weights = dd))$residuals)     [1] 2.278805     dd2 <- c(rep(1, 250), rep(2, 250))     sd(lm(yy ~ xx, weights = dd2)$residuals)         [1] 50.29102         sd(summary(lm(yy ~ xx, weights = dd2))$residuals)     [1] 62.66624      here is my decay function, by the way:               decay <- function(len, halflife, sumone = TRUE) {     #function generates an exponentially decaying series         t <- len:1 # generate a series of numbers reverse order so biggest weights last         lambda <- log(2) / halflife #figure out the lambda for the halflife         w <- exp(-lambda * t) #create the weights series           if(sumone) w <- w / sum(w) #normalise sum to 1 if necessary         return(w)      }


title: What is the manifold assumption in semi-supervised learning?

text: I am trying to figure out what the manifold assumption means in semi- supervised learning. Can anyone explain in a simple way? I cannot get the intuition behind it. It says that your data lie on a low-dimensional manifold embedded in a higher- dimensional space. I didn't get what that means.


title: What are typically encountered condition numbers in social science?

text: As part of my thesis, I'm proving (or attempting to prove...) a few asymptotic results. Because these results depend on the condition number, I'd like to have some idea about the typical sizes of a condition numbers that crop up in social science research. That way, I can give some guidance about how large the sample size has to be before we reach the happy land of asymptopia. I'd be happy for any guidance. **My very specific** setup is as follows. For the standard Generalized Least Squares (GLS) model $$Y = X\beta + e \quad \quad \quad e \sim N(0, V\sigma^2) $$ where $V$ is assumed to be known and positive definite, we define $$ X^- = (X^\top X)^{-1} X^\top \quad \quad \quad U = (I-XX^-)V$$ and the condition number $\kappa$ $$ \kappa = \frac{ \lambda_{\text{max}} }{ \lambda_{\text{min}} } $$ where the $\lambda_\star$ values are the maximum and minimum eigenvalues of the matrix $U$. Does anyone have pointers to references for the sizes of condition numbers in social science research? I don't even know where to look. Any pointers for either   1. OLS estimators (used incorrectly in a GLS context as posed above)   2. GLS estimators (correctly analyzed)    3. REML/ML estimators where $V$ is estimated and then conditioned upon, or    4. OLS fixed effect only models where $V$ is the identity matrix would be most welcome!


title: How to compare pre-post outcome when each subject react differently?

text: I have 11 subjects with pre and post outcome (two repeated measures). I'd like to test the difference between both conditions. The problem is that each subject react very differently. Here's an example :                 Subjects Pre-score Post-score              S646     1.52      8.18         S657    -1.43     -3.69         S680     1.57      1.33         S689     1.46      0.38         S690     2.02      3.48         S691     4.26      8.79         S692     2.20      2.83         S695    -0.88     -0.28         S699     5.05      5.13         S701     2.90      1.04         S711     2.17      2.09      There is something happening, but the mean and a t-test reveal nothing. Is there a way to conceptualize or analyze that kind of "effect" between two conditions on the same subjects ? If I use the absolute value between the scores, |PostScore - PreScore|, then I use an independent t-test to evaluate if the absolute difference is different from zero, can I say that I examine the magnitude of the effect without consideration of the direction of the effect? Or maybe, I can use an effect size with that kind of data, with pre-post or repeated measures ? In fact, I'm interested by the magnitude of the effect, not by the direction.


title: How to perform training and testing using DPpackage of R for Non parametric Bayesian analysis?

text: I want to do some Non-parametric Bayesian predictive anlysis through R and I found out that in R there is `DPpackage` which has functions for Non- parametric Bayesian modeling. But I am having trouble in understanding how exactly to perform training and testing using the DPpackage as I am unable to find some examples which show such implementations. If anyone has some sort of example which shows how training and testing is performed using DPpackage in R then it will be very helpful for me to get the understanding of the whole process.


title: Likelihood ratio test to determine if average number of accidents has dropped?

text: I can't find any worked (non-trivial) practical example for a likelihood ratio test, believe me I have spent hours looking. Here is a question I've been trying to complete but I can't get any further. Can anyone show me how to complete this, and also tell if I am on the right track with the work I have done? **Question** The average time between accidents on a road is 4 days. This year a random sample of accidents showed that the time between accidents was 2, 3, 6, 1 days. Is there evidence that the number of accidents has dropped on average? Use a likelihood ratio test. **Attempt** Our hypotheses are $H_0: \lambda = 4$ $H_1: \lambda < 4$ Where $\lambda$ is the population mean time between accidents Our likelihood ratio test is $$\phi = \frac{\sup_{\lambda \in \Omega_{H_{0}}}L(\lambda | \underline x)}{\sup_{\lambda \in \Omega}L(\theta|\underline x)}$$ Now, $L(\lambda|\underline x) = \prod_{i=1}^4 f(x_i|\lambda)$ $= \lambda^n e^{-\lambda \sum_{i=1}^n x_i}$ $= \lambda^n e^{-\lambda n \overline{x}}$ So we have - $$\phi = \frac{\sup_{\lambda\in \Omega_{H_{0}}}\lambda^n e^{-\lambda n \overline{x}}}{\sup_{\lambda\in \Omega}\lambda^n e^{-\lambda n \overline{x}}}$$ We only have one choice for $\lambda$ in $H_0$ (although I'm not sure if I should be using $4$ or $\frac{1}{4}$ for $\lambda$) and $n=4$ giving us $$= \frac{4^4 e^{-4 * 4 * \overline{x}}}{\sup_{\lambda\in \Omega}\lambda^4 e^{-\lambda 4 \overline{x}}}$$ $$= \frac{256 e^{-16 \overline{x}}}{\sup_{\lambda\in \Omega}\lambda^4 e^{-\lambda 4 \overline{x}}}$$ I'm not sure what to do now. How do I deal with the sup in the denominator? Can I sub in $3$ for the $\overline x$'s due to the random observations of $2, 3, 6$ and $1$? Has my work up to this point been correct?


title: how calculate inflation observed and expected p-values from uniform distribution in QQ plot

text: I am trying to quantify degree of inflation (ie. how best the observed data points fit to expected). One way is too look at QQ plot. But I would like to calculate some numeric indicator for inflation - means that how well the observed fits the theoretical uniform distribution. ![enter image description here](http://i.stack.imgur.com/Y0H1d.jpg) Example data:               # random uniform distribution      pvalue <- runif(100, min=0, max=1)     # with inflation expected i.e. not uniform distribution       pvalue1 <- rnorm(100, mean = 0.5, sd=0.1)


title: Residual autocorrelation versus lagged dependent variable

text: When modeling time series one has the possibility to (1) model the correlational structure of the error terms as e.g. an AR(1) process (2) include the lagged dependent variable as an explanatory variable (on the right hand side) I understand that their are sometimes substantial reasons to go for (2). However, what are the **methodological reasons** to do either (1) or (2) or even both?


title: Measure of dispersion over unordered set

text: I'm looking for a measure of dispersion, such as standard deviation, that can be used when distributing to an _unordered set_. Specifically: A _bucket distribution_ assigns a non-negative value to each _bucket_ in a finite set. The sum of all assigned values is one. (So far, it's like a probability distribution). However: The buckets in the set can be distinguished but have no order. It's this which makes it different than a probability distribution. We can represent bucket distributions as lists of decreasing non-negative numbers that sum to 1. Eg [1] or [1/2,1/2] or [1/3,1/3,1/6,1/6] or [1/2,1/4,1/8,1/16...]. I'd like to be able to measure how _dispersed_ a particular _bucket distribution_ is. Intuitively, [1] has 0 dispersion, [9/10,1/10] has some, [9/10,1/20,1/20] has more, [1/3,1/3,1/6,1/6] has more, etc. But I haven't been able to quantify this. I've tried using standard deviation, variance, moment of inertia, etc. But I can't find a good way to do this. I can arbitrarily rank order the buckets, but this seems, well, _arbitrary_. * * * _(I'm the original poster, but can't seem to comment any more)_ Ray: Great work, thanks for sharing this original research. If you can provide more info on how you developed it, it would be fascinating. whuber: Do you still feel entropy is a better measure than $k'$? Why? Intuitively, $k'$ fits like a glove. _Just noticed:_ The Wikipedia entry http://en.wikipedia.org/wiki/Diversity_index gives both $k'$, entropy, and other measures of diversity. But, alas, no comparison (what are fundamental assumptions of each, when is one better suited, etc.)


title: Comparing observed Vs predicted window overlap/intersection

text: I have a large number (~1.5 million) of protein sequences, each of them of different lengths.There are 6 schematic examples in the attached image. Within each of these sequences, there are >= 0 domains (shown as colored boxes). Each domain is a match to one amongst a collection of approx 15000 profile HMMs (D1 ... D15000). Domains cannot overlap. The sum of length of all domains identified in a protein <= length of entire protein. Some proteins may not possess a domain D type region in them. Apart from these 15K 'D type domains', I have also detected matches to a different type of region in each protein, lets call it 'IDR type'. These IDRs are shown as red lines above the protein sequences. Some proteins may not possess an IDR type region in them. The IDRs may be of different lengths, but the sum of their lengths in one protein, cannot be longer than the length of the protein.The location of these IDRs can be anywhere between the start and end of the protein, the length of each individual IDR is also quite variable. Two IDRs will not overlap with each other. Now I want to know how to formulate and statistically test my hypothesis that 1. the IDR(s) in each protein (across the entire collection of proteins) do NOT overlap AT ALL with domain D-type regions (i.e. red lines and colored boxed never overlap). and / or the IDR(s) in each protein (across the entire collection of proteins) do NOT significantly intersect with domain D-type regions (i.e. red lines and colored boxed never overlap). In other words, is the frequency of intersect between D type regions and IDR type regions significantly lower or higher than expected just out of chance....and what is the statistical support for that assertion? ![pictorial representation of dataset and problem](http://i.stack.imgur.com/WIcHM.jpg) For each of the proteins in question, I know the start/stop locations of each of the D type matches, and the IDR type regions.But I have not processed this information to know, even casually, what the reality of these overlaps are. In the attached image that I made up with imaginary data, the position of the IDRs (red lines) is completely non-intersecting with the domains (colored boxes). Which is not necessarily the case with the real data I have. Please note that the position / location / coordinates of any "domain of D type" or "region of IDR type" is random. Could the answer be something like what is explained athttp://www.mathpages.com/home/kmath580/kmath580.htm. Or perhaps I am overthinking this? Also, I suppose because I am testing the hypothesis on multiple proteins, there needs to be some sort of multiple testing correction? Thank you in advance for your help.


title: Function of random variable

text: If I have a random variable $X$ which has mean $\mu$ and variance $\sigma^2$, what is the approximate expression of $log(X)$ and $\sqrt{X}$? Do I assume normal approximation or use Taylor expansion?


title: kernel density estimation of a Dirichlet distribution

text: I have a Dirichlet distribution over the parameters of a multinomial, and I want to estimate its posterior density given some set of evidence. Due to some pecularities of my problem (e.g. presence of "soft evidence", among others), I unfortunately cannot use a simple update counts on the parameters, but I have to resort to sampling techniques. Once I have collected enough samples for the posterior distribution, I then need to re-estimate its density. Note that I don't necessarily need to estimate a Dirichlet again, the only thing I need is a probability density function. I thought of two possible strategies, one parametric and one non- parametric:   1. Try to re-estimate the Dirichlet parameters (alpha counts) from the samples. It seems theoretically possible (see here for instance), but quite difficult to implement, since there is no closed-form solution.   2. Alternatively, I could give up trying to re-estimate the Dirichlet counts entirely, and resort to kernel density estimations based on my samples.  I'm currently trying to implement the second strategy, based on a kernel density estimation with Gaussian kernels (see e.g. here). I unfortunately don't seem to get the right result, the calculated densities being too high. As I understand it, the problem might be that the samples of my Dirichlet are not standard multivariate values: namely, they are constrained to be between 0 and 1 and sum up to 1. They are therefore very _strongly correlated_. For instance, for a Dirichlet with 2 dimensions (aka a Beta distribution), if the first variable has a value 0.6, the second variable is constrained to have a value 0.4. All the other values have a null density. I am therefore wondering if it is possible to perform a multivariate kernel density estimation for values drawn from a Dirichlet distribution. If yes, how do you think I should do it? Should I use another kernel? If no, do you have any alternative to KDE to suggest? Thanks!


title: Establishing relationship between discrete vector and categorical vector

text: I have a vector of positive integer counts taken from a time-series (365 observations/days). I believe that there is variance in these counts due to the day of week they are observed and as I am looking for abnormal observations I'd like to remove this variance before placing a threshold on what is abnormal. I have tried creating a vector of days (1-7) which goes from t1 to t365 and counting the activity on those days. I then divided each of the count observations by the activity for that day. Is this the right approach and how do I measure the strength of the correlation? (I have looked at point-biserial correlation but is that just for binary categories?) Ultimately I would like to see how much of the variance in the discrete vector is explained by the day of the week.


title: Can we use as predictor a variable that was used in the calculation of the dependent (a ratio)?

text: I wonder if someone could give me some advice on using ratios as a dependent variable in a Generalized Linear Model. I have a variable referring to the increase of "size at Time 1" to "size at Time 2" for $n=37$ individuals. I want my dependent variable to be percentage increase from Time 1 to Time 2. Because "size at Time 2" is equal to or greater than "size at Time 1" my dependent may range from 0 to +infinity. I'm interested in testing a number of potential explanatory variables for this percent increase while also taking into account a potential confounding factor (that will be represented as a continuous variable). I've read that I might log-transform my dependent variable and use a GLM to conduct the analysis. This would allow me to include my confounding factor as a covariate in the analysis in order to measure its effect. My question is can I also test the influence of "size at Time 1" as an explanatory variable? What are the implications of using as predictor a variable previously used for calculating the dependent?


title: Negative and positive AICc/BIC for two models with transformed data - how to compare?

text: I am using AICc for model selection for transformed data (continous variable). One model I used $\log_{10}(PWV)$ as response and the other $\log(PWV)$ as response, but I am not sure which one to use if I compare the AICc. Se underneath:               Response Log10(PWV)     AICc    BIC     -2745,2 -2699,71          Response Log(PWV)     AICc    BIC     907,8648    953,3541      Should I use the model with the value that is smallest absolute value no matter negative or positive, that is, the $\log(PWV)$ or the model with the smallest value, that is $\log_{10}(PWV)$. I also get values of 2 LogLikelihood, which I am not sure whether I can use instead for to comparison. **Ekstra information about my very exciting study:** Using SAS. JMP to my mixed model. My study is to validate an apparatus that measures PWV (pulse wave velocity), and to show reproducebility and intra/interobserver variation. We measured 8 minipigs, two observers, several times on each pig, and we repeated this 3 times on different days. So have an unbalanced data set, with Pigs as random effect, Observer (2 levels) and examinationday (3 levels) as fixed effect, and then a lot of variables. My data shown no normaldistribution by Shapiro-Wilks test, so I tried to transform by log(PWV), log10(PWV) and root(PWV) and then make the models with transformed data to se if residual error is normaldistributed. But don't really understand how to do this in JMP, so I'm looking for the answer to do this. How do I check residual error for normaldistribution in JMP? Further what do I use I my residual aren't normaldistributed? Should I use Wilkinson, and is that fine with a repeated measurement as mine? This is the test with PWV not transformed.               ***Response PWV***          Summary of Fit          RSquare 0,656081     RSquare Adj 0,651699     Root Mean Square Error  0,955595     Mean of Response    6,236935     Observations (or Sum Wgts)  796          AICc    BIC     2256,695    2317,064          Parameter Estimates          Term        Estimate    Std Error   DFDen   t Ratio Prob>|t|     Intercept       5,9679039   0,665618    2,001   8,97    0,0122     W[26]       -2,485353   1,603072    2,001   -1,55   0,2611     W[27,5]     0,7102739   1,602099    1,996   0,44    0,7009     W[28]       -0,069661   1,602321    1,997   -0,04   0,9693     W[28,5]     1,3252744   1,60304 2,001   0,83    0,4953     W[29]       0,9563412   1,226491    1,995   0,78    0,5173     E[1]        0,5646209   0,061558    783 9,17    <,0001     E[2]        -0,444435   0,055241    783,1   -8,05   <,0001     O[1]        -0,191662   0,039188    783 -4,89   <,0001     E[1]*O[1]       0,2489414   0,060013    783 4,15    <,0001     E[2]*O[1]       -0,25976    0,054522    783,1   -4,76   <,0001               REML Variance Component Estimates          Random Effect   Var Ratio   Var Component   Std Error   95% Lower   95% Upper   Pct of Total     P[W]    3,4790351   3,1769234   3,1893213   -3,074031   9,4278783   77,674     Residual        0,9131622   0,0461511   0,8290524   1,0108232   22,326     Total       4,0900856               100,000                 -2 LogLikelihood =      2230,229647          Fixed Effect Tests          Source  Nparm   DF  DFDen   F Ratio Prob > F          W   5   5   1,994   0,6658  0,6917       E   2   2   783,1   45,5276 <,0001       O   1   1   783 23,9207 <,0001       E*O 2   2   783,1   12,2427 <,0001        Here is another test with transformed PWV. _**Response Log(PWV)_**               Summary of Fit          RSquare 0,693231     RSquare Adj 0,689323     Root Mean Square Error  0,151507     Mean of Response    1,795197     Observations (or Sum Wgts)  796          AICc    BIC     -634,624    -574,255          Parameter Estimates          Term        Estimate    Std Error   DFDen   t Ratio Prob>|t|     Intercept       1,7446075   0,109662    1,999   15,91   0,0039     W[26]       -0,484959   0,26411 2   -1,84   0,2078     W[27,5]     0,1478977   0,263961    1,995   0,56    0,6318     W[28]       0,01762 0,263995    1,996   0,07    0,9529     W[28,5]     0,2290124   0,264105    2   0,87    0,4773     W[29]       0,1623861   0,202078    1,994   0,80    0,5062     E[1]        0,082577    0,00976 783 8,46    <,0001     E[2]        -0,060043   0,008758    783,1   -6,86   <,0001     O[1]        -0,03304    0,006213    783 -5,32   <,0001     E[1]*O[1]       0,0288572   0,009515    783 3,03    0,0025     E[2]*O[1]       -0,03895    0,008644    783,1   -4,51   <,0001               REML Variance Component Estimates          Random Effect   Var Ratio   Var Component   Std Error   95% Lower   95% Upper   Pct of Total     P[W]    3,7578306   0,0862586   0,086597    -0,083468   0,2559857   78,982     Residual        0,0229544   0,0011601   0,0208401   0,0254093   21,018     Total       0,109213                100,000            -2 LogLikelihood =      -661,0898412               Fixed Effect Tests          Source  Nparm   DF  DFDen   F Ratio Prob > F          W   5   5   1,996   0,8590  0,6156       E   2   2   783,1   37,1920 <,0001       O   1   1   783 28,2795 <,0001       E*O 2   2   783,1   10,1548 <,0001        Generally the Residual by Predicted Plot doesn't look really nice, but they aren't really bad. But I want to now how to see whether the residual is normaldistributed? Hope to get some help.


title: non-EM algorithm approach to mixture model?

text: I have a mixture model and the components are further parameterized by ~200 variables. Originally I use EM-algorithm to get a MLE estimation of the parameters. The algorithm works quite well and converges quickly. However, when I scale up the problem to ~400 variables, the problem become ill-conditioned. After some regularization, I have already overcome the numerical difficulty and I could get a meaningful answer, but the computation time of EM algorithm is too long for practical use. I am currently thinking whether I should switch to gradient-based method. But the gradient of the log- likelihood is very complicated (and may be numerical unstable too) if I do not employ EM-algorithm. May I have some advice on what I algorithm or method I should look for? P.S. I don't want to use MCMC / Fully Bayesian approach to the problem as I want a MLE estimation.


title: SEM or path analysis using AMOS, issues with sample size and model complexity

text: I am new to SEM / AMOS and as such I have done a lot of reading around my problem but find I am getting lost, any help would be appreciated and I am sorry if this sounds rather basic;   1. I have a model that looks at a number of factors (distal, indirect factors) which include 5 factors of personality and a couple of 'other factors' shown to be related to 2 distinct accident behaviours (outcome variables, DVs etc)   2. There are 2 mediating factors between the distal and outcome variables   3. the mediating variables and outcome variables are bi directional so I have an non recursive model.   4. I have relatively small sample of 200 for the purposes of SEM depending on the perspective you take on this. I understand the principles of SEM and here lies my problem, to include a measurement model with observed variables that underlie these factors is not possible (the 5 individual factors are measured by a 120 item questionnaire), this puts unrealistic demands on sample size and simply drawing a model like this seems impossible. I was looking at the use of Path Analysis which I understand is a form of SEM. Here I treat the factors as observed variables so lose the ability to account for measurement error. There are also issues with testing a non-recursive model here. I should say that I conducted a series of multiple regressions to establsh the best predictors of the 2 mediating and 2 outcome variables to form I suppose a 'conceptual' model that I would like to test and test the mediating variables. I hope this makes sense to someone..The more I read, the more I am getting a little overwhelmed and was hoping someone could offer some guidance or comments? Many thanks in advance


title: Should one compute multiple q-values when power of tests differs or one global q-value?

text: **Background:** If I have 40 hypotheses that represent the same statistical test based on equally large data sets, I can compute 40 q-values (Benjamini-Hochberg). **True Problem** What if I have 100 sets of these 40 hypotheses (for a total of 4000 hypotheses)? Each set of 40 hypotheses (of which there are 100) may incorporate a lot less data, meaning not all 100 sets of hypotheses have the same statistical power. However, within each set, all 40 hypotheses do have the same statistical power. Wouldn't computing a global q-value based on all 4000 hypotheses be overly conservative? Because the q-value is in reality a rate, wouldn't it make sense and be perfectly valid to simply compute each q-value based on the 40 hypotheses in its set rather than on the global set of 4000 hypotheses?


title: Change point analysis using R's nls()

text: I'm trying to implement a "change point" analysis, or a multiphase regression using `nls()` in R. Here's some fake data I've made. The formula I want to use to fit the data is: $y = \beta_0 + \beta_1x + \beta_2\max(0,x-\delta)$ What this is supposed to do is fit the data up to a certain point with a certain intercept and slope ($\beta_0$ and $\beta_1$), then, after a certain x value ($\delta$), augment the slope by $\beta_2$. That's what the whole max thing is about. Before the $\delta$ point, it'll equal 0, and $\beta_2$ will be zeroed out. So, here's my function to do this:               changePoint <- function(x, b0, slope1, slope2, delta){         fit <- b0 + (x*slope1) + (max(0, x-delta) * slope2)     }      And I try to fit the model this way               nls(y ~ changePoint(x, b0, slope1, slope2, delta),          data = data,          start = c(b0 = 50, slope1 = 0, slope2 = 2, delta = 48))      I chose those starting parameters, because I _know_ those are the starting parameters, because I made the data up. However, I get this error:               Error in nlsModel(formula, mf, start, wts) :        singular gradient matrix at initial parameter estimates      Have I just made unfortunate data? I tried fitting this on real data first, and was getting the same error, and I just figured that my initial starting parameters weren't good enough.


title: Determine over/under of regression prediction

text: I've created a regression analysis in Excel and using the model I have created a prediction point for Y for a given input of X. I'd like to add an "over/under" to the Y value to show the range that Y could reasonably be expected to fall in. Is using Y +/- the "standard error value for y" (Sey from the LINEST function) a reasonable approach or is there a better way?


title: R code for two sample t-test in case of equal variances

text: In R, how can I perform a two sample t-test in the case of equal variances? My actual problem is to compare different algorithms for length measurements in machine vision. I have studied _Applied Statistics and Probability for Engineers_ by Montgomery & Runger and I have tried to implement the test in R and compare the results to the R command `t.test`. This is my code:               # Douglas C. Montgomery, George C. Runger     # Applied Statistics and Probability for Engineers     # Third Edition     # 10-3.1 Hypotheses Tests for a Difference in Means, Variances Unknown, p.337     # Case 1: sigma_1^2 = sigma_2^2 = sigma^2     #     # See also: EXAMPLE 10-5, p.339     two_sample_t_test_equal_variance <- function(x1,x2,Delta_0,alpha)     {        n1 <- length(x1)        n2 <- length(x2)        dof <- n1+n2-2        S1_squared <- var(x1)        S2_squared <- var(x2)        #pooled estimator:        Sp_squared <- ((n1-1)*S1_squared + (n2-1)*S2_squared)/(n1+n2-2)         x1_bar <- mean(x1)        x2_bar <- mean(x2)        Sp <- sqrt(Sp_squared)        t0 <- (x1_bar - x2_bar - Delta_0)/(Sp*sqrt(1/n1+1/n2))        t_half_alpha_dof <- -qt(alpha/2,dof)        reject_H0 <- (t0 > t_half_alpha_dof || t0 < - t_half_alpha_dof)        if ( reject_H0 ) {           cat("Reject H0 (alpha =",alpha,").\nH0 is mu1 - mu2 =", Delta_0, "\n\n" )        } else {           cat("Cannot reject H0 (alpha =",alpha,").\nH0 is mu1 - mu2 =", Delta_0, "\n\n" )        }             test_result <- t.test(x1, x2, alternative="two.sided", mu=Delta_0,                                paired=FALSE, var.equal = TRUE )             accept_H0_by_R <- (test_result$p.value > alpha)        if ( reject_H0 != !accept_H0_by_R ) {           cat("WARNING: R t.test gives a different answer for accepting H0.\n" )        }             rel_error_warn <- 0.01        if ( abs(test_result$statistic - t0)/abs(test_result$statistic) > rel_error_warn ) {           cat("WARNING: t0 relative error is >", rel_error_warn, "\n"  )        }     }          # data from EXAMPLE 10-5, p.339     cat1<-c(91.50, 94.18, 92.18, 95.39, 91.79, 89.07, 94.72, 89.21)     cat2<-c(89.19, 90.95, 90.46, 93.21, 97.19, 97.04, 91.07, 92.75)     two_sample_t_test_equal_variance(cat1,cat2,0,0.05)      Is my invocation of t.test correct for the case at hand? Is there any error in the code? I tagged the question with "homework" even if I am not at school since more than 10 years and I am self-studying statistics for job-related reasons. Thank you. Alessandro


title: Examples of text mining with R (tm package)

text: I spent three days dabbling with `tm` after reading a draft paper by a friend where he explored a text corpus with UCINET, showing text clouds, two-mode network graphs and Single Value Decomposition (with graphics, using Stata). I ran under a large number of issues: on Mac OS X, there are issues with the Java behind libraries like Snowball (stemming) or Rgraphviz (graphs). Could someone point out not packages – I have looked at `tm`, `wordfish` and `wordscores`, and know about NLTK – but **research, if possible with code, on textual data, that successfully uses`tm` or something else to analyse data like parliamentary debates or legislative documents?** I cannot seem to find much on the issue, and even less code to learn from. My own project is a two-month parliamentary debate, with these variables informed in a CSV file: parliamentary session, speaker, parliamentary group, text of oral intervention. I am looking for divergence between speakers and especially between parliamentary groups in the use of rare and less rare terms, e.g. "security talk" against "civil liberties" talk.


title: Variability of a curve with 4 parameters

text: Say I have 10,000 data in 2-D and I want to fit a curve to them. There are many functional forms this curve could take -- polynomial, B-spline, trigonometric, and so on. I've decided that I only want to use 4 parameters. Is there a way to figure out what is the most accurate functional form? That is, considering all possible functions with 4 parameters, which one fits the best in, e.g., an $L_2$ sense? _edit: maybe I should ask about the most accurate function with the same V-C dimension as a polynomial of degree 4 rather than "4 parameters"?_


title: How to specify the null hypothesis in hypothesis testing

text: What is a good rule of thumb for how to choose the question for the null hypothesis. For instance, if I want to check if the hypothesis B is true, should I use B as the null, B as the alternative hypothesis, or NOT B as the null? I hope the question is clear. I know that it has something to do with the error I want to minimize (Type I?), but I keep forgetting how it goes, because I don't have a clear intuition built for it. Thanks.