title: What if your randomly formed groups are clearly not similar?

text: What if, before you begin the data collection for an experiment, you randomly divide your subject pool into two (or more) groups. Before implementing the experimental manipulation you notice the groups are clearly different on one or more variables of potential import. For example, the two (or more) groups have different proportions of subjects by gender or age or educational level, or job experience, etc. What is a reasonable course of action in such a situation? What are the dangers of discarding the original random division of the subject pool and dividing the pool again? For example, are the inferential statistics that you might calculate based on the second set of groups in any way inappropriate due to the discarded first set of groups? For example, if we subscribe to discarding the first division of the subject pool into groups, are we changing the sampling distribution that our statistical test is based on? If so, are we making it easier or harder to find statistical significance? Are the possible dangers involved in repeating the division of subjects greater than the obvious danger of confounding due to group differences in educational level, say? To make this question more concrete, assume for the sake of this discussion that the topic of the research is teaching method (and we have two teaching methods) and the difference noted between the two groups of subjects is level of formal education, with one group containing proportionally more people with highest educational attainment of high school level or less and the other group containing more people with some college or a college degree. Assume that we are training military recruits in a job that does not exist in the civilian world, so everyone entering that specialty has to learn the job from scratch. Assume, further, that the between group imbalance in previous educational attainment is statistically significant. Parenthetically note that this question is similar to What if your random sample is clearly not representative?. In a comment there, @stask perceptively noticed that I am a researcher not a surveyor and commented that I might have gotten more relevant answers had I tagged my question differently, including "experiment design" rather than "sampling." (It seems the sampling tag attracts people working with surveys rather than experiments). So the above is basically the same question, in an experimental context.


title: Filtering virtual economy transactions

text: I am currently writing a pet-project for a browser based game which consists of two parts, first part is a data miner which observes the transactions in the market. And second part is as you might guessed the analyzer. But between these two steps I want to filter out some transactions since people sometimes use the market to transfer gold between each other and that corrupts the data when someone pays 10x average price of an item. I haven't taken probability and statistics course yet so I am not sure what to use when filtering out those. First thing I thought of is to filter out 2x average prices but that doesn't feel right, what would you propose? On second thought I came up with filtering out average + 1 std dev. but I don't know what that border means in words, I mean I don't have the insight on it. At the end I am willing to come up with graphs and maybe functions of things like how does the average price and the median of prices change throughout the day(when do they approach eachother and when do they drift apart), or how differently the market behaves in weekdays or weekend. This all sounds exciting to me but first I have to filter my data set.


title: OLS: $E[\epsilon_{it}^T\epsilon_{it}] \not= 0$ in 1st equation biases standard errors in 2nd equation?

text: Suppose ${X_{it}},{Y_{it}}$ are time series with $X_{it}\sim N(0.1,1)$, ($\sigma^2(Y_{it}) = 1$ and $mean(Y_{it})$ is similar to that for $X_{it}$, but changes when the dummy = 1). and $t \in \\{1,2,...,200\\}$, $i \in \\{1,2,...,N\\}$. In a real world setting this will be periodic stock market returns over $N$ firms (but you can ignore this). There is a dummy, $D_t$ that's equal to unity over $t \in \\{150,151,...,200\\}$ and equal to zero otherwise. The time series model to be estimated with OLS $\forall i$ is: $(1) Y_{it} = \alpha_i + \beta_i X_{it} + \gamma_i D_{t} + \epsilon_{it}$ This model generally adheres to the Gauss-Markov assumptions for each $i$. However, we have $E[\epsilon_{it}^T \epsilon_{jt}] \not= 0$ for all $i$ and $j$. The next step is to construct a vector of gammas using the $N$ estimates of model $(1)$. Call this vector $\bf{\hat{\gamma}}$. We then use this in the cross-sectional model: $(2) \hat{\gamma}_i = a + b Z_i + u_i$ where $Z_i$ is some cross-sectional variable that doesn't cause any violations in OLS assumptions and is relevant for explaining $\hat{\gamma}_i$. * * * The claim in the applied econometrics literature is that $E[\epsilon_{it}^T \epsilon_{jt}] \not= 0$ in model $(1)$ leads to (i) No problem for the OLS coefficient estimates in $(2)$, but (ii) Biased standard errors in $(2)$.   * Can someone please post ideas about why this is the case?    * I don't understand what $\epsilon_{it}^T$ is in the expression $E[\epsilon_{it}^T \epsilon_{jt}] \not= 0$. Of course $\epsilon_{it}$ is a scalar and you can't transpose a scalar. This is seen HERE, where they apply this methodology.


title: Is optimizing α for partitions of a scale a sensible method for factor analysis? (was: relation between consistency and unidimensionality)

text: > **Cluelessness Disclaimer** > > I'm a statistics noob so if at all possible, please don't stone me. Also > write slowly and in simple terms. I'm wondering what the relation between the internal consistency of a scale as usually measured by Cronbach's $\alpha$ and unidimensionality is. The concrete reason I'm asking this fundamental question is that I tried to split a scale into three factors which separated quite nicely using a PCA and a direct oblimin rotation. Almost no intercorrelation between factors and the worst loading in any scale is $.59$. Now while the first new scale that fell out of this process has a nice $\alpha$ of $.81$ (the total scale had $.5$), the second scale has a whoppingly bad $\alpha$ of less than minus frickin' six point five! I've seen my fair share of negative alphas in my short statistical life but this is unprecedented. Now my naïve understanding is that internal consistency measures the existence of an underlying construct and unidimensionality of a scale means that this underlying construct has only one dimension. This would mean that internal consistency is a prerequisite for unidimensionality. This SSPS FAQ item at least seems to imply this and this section in Wikipedia that I found while Googling might imply anything because it's really poorly written but it refers to this paper that unfortunately goes over my head. Because this post got so lengthy, this is probably a good time for a short break for your eyes and brain. ![lolcat](http://i.stack.imgur.com/msr29.jpg) Sorry, this looked like a funny idea on paper. Anyway, there is also this question Assessing reliability of a questionnaire: dimensionality, problematic items, and whether to use alpha, lambda6 or some other index? I already found that I feel might contain the answer I'm seeking but it doesn't fully reveal itself to me. Based on this understanding however, I started with $\alpha$ and tried to work my way from there. What I did is somewhat … unorthodox (which probably translates to "wrong" but I hope we'll hear more on that in an answer to this post). I did a global optimization of $\alpha$ which means I went through all possible partitions of the items of my scale and computed the average $\alpha$ for each partition. Funny thing is, the best solution for three factors (three and four factors had an identical best average $\alpha$ of $.71$) actually makes sense in terms of interpretability (the whole PCA shebang didn't). I'm pretty sure this is a coincidence though, right? * * * **Edit:** After StasK's comment, I decided to include my code for the $\alpha$ optimization in case others might find it useful. Just put your scale into a file named `data.csv` with columns as variables and rows as cases. This script's not pretty and it's not fast; it did the job I personally wanted it to do. The optimization of my 12 items took almost 12 minutes on my Core i3. There can be optimizations done in the code but they probably won't fundamentally change the fact that optimizing a 30-item scale with this method would take not 30 minutes as one might naïvely assume but roughly 300 times the age of the universe. Blindly partitioning your 100-item questionnaire is probably not a good idea with this method. The number of partitions generated by `partitions()` can be checked here http://www.wolframalpha.com/input/?i=bell+number+of+12. Just replace 12 by your number of items. If you only want to generate partitions with a certain number of factors, replace `for p in partitions(range(len(variances))):` by `for p in n_partitions(range(len(variances)), 3):` for 3 factors say. This won't make it faster though (here a considerable speedup would probably be possible if you restrict the search depth in the first place). In the long run, only a cleverer way than to test all possible partitions will do for larger scales.               #!/usr/bin/env python          from numpy import *     import sys          def variance(data):         n = len(data)         variance = 0.0          mean = float(sum(data))/n         for x in data:             variance += (x - mean)**2         return variance/n          # from http://stackoverflow.com/q/2037327/1050373     def partitions(set_):         if not set_:             yield []             return         for i in xrange(2**len(set_)/2):             parts = [set(), set()]             for item in set_:                 parts[i&1].add(item)                 i >>= 1             for b in partitions(parts[1]):                 yield [parts[0]]+b          def n_partitions(set_, n):         for partition in partitions(set_):             if len(partition) == n:                 yield partition          def alpha(data, variances, cols):         n = len(cols)         cols = array(cols)         data_cols = data.transpose()[cols]         variances = variances[cols]         data_rows = data_cols.transpose()         row_sums = array([sum(row) for row in data_rows])         return n/float(n-1)*(1-sum(variances)/variance(row_sums))          data = genfromtxt('data.csv', delimiter=',')          print 'input data:'     print data          # pre-compute column variances     variances = array([variance(col) for col in data.transpose()])          data = data[0:len(data)]          print 'overall alpha:'     print alpha(data, variances, range(len(variances)))          best_alpha_average = float('-inf')          for p in partitions(range(len(variances))):         lengths = [len(subset) for subset in p]         if min(lengths) < 2:             continue         alphas = [alpha(data,variances,list(subset)) for subset in p]         alpha_average = sum(alphas)/len(p)         if alpha_average >= best_alpha_average:             print 'new best:',alpha_average,alphas             print 'partition of size '+ str(len(p)) + ': ' + str([list(array(list(subset))+1) for subset in p])             best_alpha_average = alpha_average             sys.stdout.flush()


title: Exponential Distribution - Rate - Bayesian Prior?

text: I have gone through WinBugs documentation (for example, http://www.mrc- bsu.cam.ac.uk/bugs/thebugsbook/examples/html/Chapter-11-Specialised/Example-11_7_2-leukaemia.html). And also through this book (http://www.amazon.ca/Bayesian-Survival-Analysis- Joseph-Ibrahim/dp/0387952772). Both use a gamma distribution prior for the rate (lambda) for the exponential distribution. They alternate between Gamma(0.01,0.01) and Gamma(0.001,0.001). I would like to use a non-informative prior.... but I don't think this is it? Could someone explain whether this is a non-informative prior? If not, can anyone suggest one?


title: Error on weighted mean

text: I have some measured datapoints. Each of the points has attached an error, depending on how precise the measurement was done. Now I'm using the squared weighted mean with the inverse of the errors as weights so that I emphasize more on the preciser values and less on the values where the error is large. What I'm wondering is what can I say about the precision of my squared weighted mean. Is there a way to calculate the expected error range for it?


title: reconstructing matched data set from R's matching package

text: I am using Mebane and Sekhon's Matching package for R, and my goal is to make casual inferences about a proportion and a mean. The proportion is straight forward, but the mean has missing values, which `Match` does not accept. My approach is trying to get access to the matched records through the row indexes, but the results are not quite right. Example code:               require(Matching)     demo(DehejiaWahba)     # mean before matching     by(lalonde$u75, lalonde$treat, mean)     # mean after matching     matched <- rbind(lalonde[rr$index.treated,],        lalonde[rr$index.control,])     by(matched$u75, matched$treat, mean)      The output of the demo shows the matched means for u75 are 0.6 and 0.62072:               ***** (V14) u75 *****                        Before Matching           After Matching     mean treatment........        0.6                       0.6      mean control..........    0.68462                   0.62072       However, the reconstructed data set has different means:                   matched$treat: 0         [1] 0.7861272         ---------------------------------------------------------------------------------          matched$treat: 1         [1] 0.7543353


title: Cost function NN with weight derivation

text: In a single layered neural network with a sigmoid function (to make it easy to understand) the cost function is $E_j = \frac{1}{2} \sum_{k=1}^{K}(\text{target}_{jk} - \text{observed}_{jk})^2 + a\sum_{i=1}^{I}w_{ij}^2.$ So this is a two part problem and I know how to derive the first half at both the hidden and output layers. But, for the second half, I am not sure. Any advice or help would be great. Thank you.


title: R statistics: Adjusted Chi-squared test for clustered binary / categorical data

text: I'm looking for some assistance in statistical analysis with R, but also some general stats advice. I am analysing cardiac phenotype data by comparing 2 groups. The 2 groups are unmatched individuals, but within each group, they are clustered in family subgroups (of between 1 and ~6). I want to report the difference in prevalence of a specific ECG appearance (binary - i.e. either present or absent in each individual) between the 2 groups. For example: > Group 1 consists of 157 individuals comprised of 41 family clusters. Group 2 > consists of 463 individuals comprised of 163 family clusters. Prevalence of > x in Group 1 = 22.9% Prevalence of x in Group 2 = 24.6%. Group 1 are cases > and Group 2 controls (i.e. not randomized and defined by phenotype in an > observational study). What test is most appropriate in this circumstance, and which package in R provides the easiest way to account for the clustering of relatives within families? Having looked around, I have found:   * Ratio estimate chi-square test   * Generalized estimating equation But I have no experience of either of these techniques, and can't find any examples of their use in R. Any advice on how best to proceed? EDIT: See comment below for update. I believe the Donner (1989) chi-square correction may be the most appropriate (provided by R function donner). Second opinions and correct use of R command appreciated. Thanks.


title: What test will tell me a normalized percentage of the data?

text: I'm working with baseball statistics and I have a list of the players' batting percentages. I want a formula that will tell me what percentage of the data is represented from a player's score compared against the others as a percentage. For example, I have `PlayerA`, he has the highest batting average in the sample, so he is represented as 100%. `PlayerB` has the second highest batting average, he is represented as 98%. The percentage is NOT dependent on the size of the population (i.e., if the total population were 50 and I wanted to measure the value from the population directly then it would go by twos down to zero). The measurement is a percentage of the range of values. So if five players all had close batting averages their percentages would all be nearly the same. Another way to phrase this would be, e.g., I am 6'2" and I know the height of each member of the population. After I plug in this formula I am able to discern that I am taller than 72% of the population. It seems like it would be related to a Z-score but I'm not positive.


title: How to add a random perturbation term to the response variable, with transformation in R?

text: I'm trying to add a random perturbation term to my response variable (to estimate the impact of multicolinnearity in the beta estimates), but I am not sure on how to accomplish it when the response variable is transformed using a 1/4 power, in R. Without any transformation in the response variable, according to Faraway's "Linear Models with R," the way to write it is:               g <- lm(hipcenter+10*rnorm(38) ~ ., seatpos).       So if I have a 1/4 power transformation for the response variable, then the code should be:               g <- lm((hipcenter^(1/4))+10*rnorm(38) ~ ., seatpos).       Is this right?


title: Combining sources of uncertainty/variation in a multi-layered linear model

text: I’m having trouble understanding how to combine parameter uncertainty and interannual variability from various levels in a bootstrapped linear model. Specifically, this model is designed to generate estimates of some phenomenon $Y$ for 0.5º grid cells in a spatial region (e.g, Central America). It looks like this: $$ \boldsymbol{Y_{m,y}} = \boldsymbol{\beta_{m,y}X_{m,y}} $$ where $\boldsymbol{\beta}$ is the vector of parameters being estimated (size 1x3), $\boldsymbol{X}$ is a matrix of observations of land cover fractions (size 3xN), $\boldsymbol{Y}$ is the vector of predicted values (size 1xN), and the subscripts indicate that the model is being parameterized for each month $m$ and year $y$ (e.g., May 2006, June 2006, and May 2007 can all have different parameter estimates). (No intercept, although I don’t think that matters for my question.) If we think of the model as $Y_{m,y}=\sum\limits_{i=1}^3{\beta_{i,m,y}X_{i,m,y}}$, then each term $\beta_iX_i$ can be considered a subset of $Y$ — i.e., $Y_i$, the amount of $Y$ that happens on land cover type $i$ — in a real physical sense. Thus, there will likely be some correlation among the various parameter estimates. This is embodied/exacerbated by the fact that for each grid cell $g$, $\sum\limits_{i=1}^3{\boldsymbol{X_{g,i,m,y}}} = 1$ — that is, the land cover type fractions add up to 1. To create uncertainty bounds for the parameters, I have done 10,000 bootstrapping runs. Let’s say there are 500 grid cells in the region. For the first bootstrapping run, I randomly choose 500 grid cells from that region with replacement. I then estimated the parameters for each month and year using that sample set. I then repeated this sampling-fitting procedure 9,999 more times. My goal is to compare the _measured/observed_ mean annual amount of phenomenon $Y$ with the _estimated/modeled_ amount. I would like the uncertainty bars around the measured amount (which will take into account just interannual variability in the measured amount) to overlap with the uncertainty bars around the estimated amount (which will take into account both uncertainty in the model parameters as well as interannual variability in the annual estimated amount). For the observations, it’s easy enough — sum up observed $Y$ for all the months in each year, then find the standard deviation across all the years. I get confused when thinking about how to do this for the estimates, though. The problem becomes even stickier if I want to compare observed vs. estimated annual GLOBAL $Y$ (i.e., across all regions). I know that adding standard deviations in quadrature is going to be key here, but the fact that there are so many levels is confusing me. I'm working in Matlab now and am also fluent in R, but any help you can provide would be much appreciated.


title: Identifying a mixture model

text: I am trying to fit a mixture model to a dataset that consists of counts (so every record is a count of something, like the number of attempts by an IP address to connect to a website). I know, a priori, that each point belongs to one of two groups (for example malicious IPs and legitimate IPs), but I don't know which one without manually inspecting the record and determining the grouping. These are the things that I know though:   * As the count value increases, it is more and more likely that the data point belongs to group 2.   * 80 percent of values are 1 and 85% are less than 2.    * Most values for group 1 are 1, but there are values as high as 10,000 that are from group 1.  Ideally, I'd like to be able to draw a line somewhere and say if the count is higher than some $n$, then I am $p$% sure that this record is from group 2. I think a mixture model is a good place to look at to model this data, but looking at the histogram and kernel density plot, I don't see any distinct clumps in data. Based on what I mentioned above, my guess is that a mixture of a log-logistic (for group 2) and pareto (for group 1) might be an ok approximation to the data. Is there an R package or python module that would estimate the parameters of such a mixture?


title: How to validate and compare models predicting a binary variable?

text: I have a question about determining which models are "better" and how to assess that info. Let's say I have three models, each which predicts our bid on won ping. Our bid is a continuous variable and won ping is either 0 (no) or 1 (yes). From the models (logit) and the system that I'm modeling, I constructed the following table with the number of bids, our win rate, predicted win rate, and average bid. Essentially, we run the model our bids and whether we win and use that to predict our win rate in the future bids.                           #_of_bids    win_rate    pred.win_rate    avg.bid               Mod_1       1792          46%       62%         $1.54             Mod_2        851          2%     2%         $0.98         Mod_3       17037         6%         0%         $2.83       According to the above table, there is a "significant" discrepancy between our win rate and the predicted win rate. As a result, I ran through the model and came up with the following "diagnostics."               Mod_1         - Won ping rate = 45.7%         - Bid distribution issues = politively skewed, but no single bid repeats significantly         - Nagelkerke's R squared = 0.23 ->  The predictor only exlained 23% of the variance found in the            response variable (won ping)  ->  poor         - Cross-Validation (10-fold) = 0.693  ->  how accurate were the model predictions (0 to 1)           - ROC Curve = 0.746 ->  fair(C)  ->  how well the model fits the data (0 to 1)      However, I was left wondering if the discrepancy between the win rate and predicted win rate, plus the validation statistics were enough to conclude that only that model was off (aka performed poorly or didn't predict the actual win rate properly). Therefore, I ran through the same process for the other two models and came up with the following stats.               Mod_2            - Won ping rate = 2.3%         - Bid distribution issues = 53.5% of our bids are $0.26 or $0.26         - Nagelkerke's R squared = 0.40 ->  The predictor only exlained 23% of the variance found in the            response variable (won ping)  ->  average         - Cross-Validation (10-fold) = 0.974  ->  how accurate were the model predictions (0 to 1)           - ROC Curve = 0.914 ->  excellent(A)  ->  how well the model fits the data (0 to 1)          Mod_3         - Won ping rate = 5.7%         - Bid distribution issues = 12.9% of our bids were either $0.25 or $0.26         - Nagelkerke's R squared = 0.03 ->  The predictor only exlained 3% of the variance found in the            response variable (won ping)  ->  horrible         - Cross-Validation (10-fold) = 0.942  ->  how accurate were the model predictions (0 to 1)           - ROC Curve = 0.339  ->  horrible  -> how well the model fits the data (0 to 1)      So in trying to figure out if it was just model 1 that was off, I looked at these various statistics. However, I'm not sure about what the proper steps are. ### Questions   1. Should I assume that because the win rate and predicted win rates were so different in the first table that the model performed poorly?   2. While the validation stats for model 1 are poor, they are not much better for model 2 or 3. What explains that result?   3. How should one go about comparing two or three different models to one another in terms of overall performance (aka fits the data/predictive quality)?


title: Shrinkage Estimator for Newey-West Covariance Matrix

text: This is a cross post. I would like to apply the Newey-West covariance estimator for portfolio optmization. Up to lag one it is given by $$ \Sigma = \Sigma(0) + \frac12 \left (\Sigma(1) + \Sigma(1)^T \right), $$ where $\Sigma(i)$ is the lag $i$ covariance matrix for $i=0,1$. Furthermore I like to use shrinkage estimators as implemented in the corpcor package for R. The identity matrix as shrinkage prior for $\Sigma(0)$ is plausible. What would you use as prior for $\Sigma(1)$ - the zero-matrix? Do you know an R implementation that allows to estimate lag-covariance matrices using shrinkage? There must be some basic difference as a lag-covariance matrix is not necessarily positive-definite (e.g. the zero-matrix). If I apply shrinkage to $\Sigma(0)$ and use the standard sample-estimator for $\Sigma(1)$ then it is not assured that $\Sigma$ is positive-definite. The above definition is taken from: Whitney K. Newey and Keneth D. West. A simple, positive semi-denite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3):703-708, 1987. It can also be found here in formula (1.9) on page 6.


title: Location, scale, and shape parameters of the lognormal (with notation ambiguity)

text: **Main question** :   If I read in a paper that a particular dataset's best fit is lognormal with $\mu=7.7$ and $\sigma=1.9$, what is the location, scale, and shape parameter of the lognormal? **Side questions** :   Is the scale parameter $e^\mu$ or $\sigma$? I see different statements in different reference literature. Is $\mu$ given in the paper equal to the mean of the transformed normal or to some new transformed $\mu$?


title: Why is Entropy maximised when the probability distribution is uniform?

text: I know that entropy is the measure of randomness of a process/variable and it can be defined as follows. for a random variable $X \in$ set $A$ :- $H(X)= \sum_{x_i \in A} -p(x_i) \log (p(x_i)) $. In the book on Entropy and Information Theory by MacKay, he provides this statement in Ch2 Entropy is maximized if p is uniform. Intuitively, I am able to understand it, like if all datapoints in set $A$ are picked with equal probability $1/m$ ($m$ being cardinality of set $A$), then the randomness or the entropy increases. But if we know that some points in set $A$ are going to occur with more probability than others (say in the case of normal distribution, where the maximum concentration of data points is around the mean and small standard deviation area around it, then the randomness or entropy should decrease. But is there any mathematical proof for this ? Like the equation for $H(X)$ I differentiate it with respect to $p(x)$ and set it to 0 or something like that. On a side note, is there any connnection between the entropy that occurs information theory and the entropy calculations in chemistry (thermodynamics) ?


title: Interpretation in log linear regressions with coefficients bigger than 1

text: Just to clarify something regarding coefficients bigger than 1 in log-linear regressions. If we have this regression, how would we go about to interpret the 1.12?   D1 is dummy variable for having a car D2 is dummy variable for being us resident lninv is the ln of the investment               Linear regression                                  Number of obs =     580                                                        F(  4,   575) =   45.58                                                        Prob > F      =  0.0000                                                        R-squared     =  0.1983                                                        Root MSE      =  1.2828                           Robust     lninv   |    Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     --------+---------------------------------------------------------------     age     |  .1376738   .0221547     6.21   0.000     .0941598    .1811878     1.car   |  .5398736   .1179852     4.58   0.000     .3081391     .771608     1.USA   |  1.124221   .1959553     5.74   0.000     .7393456    1.509097     car#USA |        1 1  |  -.0202977  .2457097    -0.08   0.934    -.5028958    .4623003     _cons   |  -.0585373  .1228857    -0.48   0.634    -.2998969    .1828223      so we have lninv = age + d1(car) + d2(usa) + d1.d2 The trouble is the interpretation of the USA. (1) Is it correct to say that all things equal being a US resident increases the investment by 112%?   Also the interaction dummy is not significant. (2) Is it ok to check for interaction between two binary variables?


title: Choice Modeling and MaxDiff

text: I am new to the market research industry and have been exposed to a procedure called MaxDiff. I have seen this also called Best-Worst scaling. I am looking for good literature overview for this method and how it fits with respect to other survey design question choices. I am especially looking for a technical analysis of the method itself rather than an application. Thanks!


title: Generate covariance matrix with fixed values in certain cells

text: I want to be able to generate a covariance matrix of dimensions $D$ x $D$, such that certain specified cells of this matrix contain a fixed predetermined values (at least approximately). For e.g. For matrix, $S$ = $$ \begin{matrix} a_{11} & a_{12} & \ldots & a_{1D}\\\ a_{21} & a_{22} & \ldots & a_{2D}\\\ \vdots & \vdots & \ddots & \vdots\\\ a_{D1} & a_{D2} &\ldots & a_{DD} \end{matrix} $$ I want to make sure that certain $a_{i,j}$ 's have a predetermined value. (If I were generating a covariance matrix without this constraint, I would just use a Wishart Random Generator. In Matlab, it would be something like - `W = wishrnd([1 0.5; 0.5 3],30)/30`) One way I can think of this problem is that different cells of the covariance matrix have different degrees of freedom. So that the cells with fixed values can be assumed to have infinite degrees of freedom and the rest as some finite number.


title: Nearest neighbor methods

text: Given a training set $\\{(x_1, y_1), \dots, (x_n, y_n)\\}$, consider the following: $$\hat{Y}(x_i) = \frac{1}{k} \sum_{x_i \in N_{k}(x)} y_i$$ Suppose $k=3$. Is this formula saying that we pick the $3$ closest neighbors of $x$ and average over their responses? So if $x_i = x_1$ then we would choose the closest $3$ values to $x_1$ (say $x_2,x_3,x_4$) and compute the following: $$ \hat{Y}(x_1) = \frac{y_2+y_3+y_4}{3}$$ Is the idea to see whether $\hat{Y}(x_1)$ agrees with $y_1$? If not, are different classification rules needed?


title: Representing categorical data in a hybrid dataset

text: I'm using Matlab to discover some clusters in a dataset that has both numerical (e.g., $) and categorical (e.g.,zip code) variables. I know some tools can handle categorical variables as factors but I'm not sure if this is the case for Matlab. I'm looking for a good way of representing my categorical variables and I'm not sure that just using dummy variables would be the best way to go. One idea is to represent each digit from my zip code, for example, as a 10 digits long binary variable. By doing so, I would have 50 new binary variables. Questions:   * Is this an acceptable way of approaching this problem?   * What other and/or more efficient ways of representing categorical variables in a hybrid dataset are there?   * Do I actually need to numerically represent my categorical variables or can Matlab actually handle the mixed dataset?


title: Test of Difference between 2 Sets of Pearson Correlations

text: Suppose I have two sets of Pearson correlation coefficients -- call them set A & set B, and they are of the same size. How do I systematically compare the correlations in A against B? E.g., I want to test the hypothesis that A is less than B -- something similar to a two-sample t-test, but the problem here is the two samples are two samples of correlation coefficients. I've done some research and found out that there is Fisher's z-transform for this purpose. But it only tests the difference between one correlation $\rho_1$ against another $\rho_2$. I couldn't find a way to systematically test the difference between to sets of correlations. Is there such a way? Or I can only compare each pair of correlations using Fisher's method and somehow derive the difference, if any?


title: Comparing two sets of time-series (physiological measures)

text: I need help in finding which statistical approach I need to use. I have two groups: G1 (disease) and control. For both groups we measured their blood level every 10 minutes: 0 to 650 min. We are trying to find the differences between G1 (disease) and control and the significance of the difference. The variables in factors are Age, Diagnosis (disease or control), and time. The dependent variable is the blood level. I do not know which approach I need to use to find this significance (difference). I tried using Multiple Regression Analysis: P = constant + B1(age) + B2(Diagnosis) + B3(Time) Diagnosis 0:constant 1:Patient If I also want to include the Interaction (Diagnosis x Time) how would I calculate that to put into Excel for analysis?


title: Similarity / dissimilarity of two large bimodal datasets

text: I am interested in assessing the divergence, or similarity or dissimilarity of 2 datasets that are the results of 2 different lidar instrument measurements. Each dataset has over 90,000 values and they do measure a continuous variable, let's say altitude. Each dataset has a bimodal distribution with a long tail on the second mode, and are highly correlated at 0.8 (Spearman). I was looking into both Kolmogorov-Smirnov test as well as Kullback-Leibler. I am very concerned that my distributions are bimodal and very large. I know that tests that will give a p-value are very sensitive to distribution types and number of measurements. For this reason I think that the KS test is not appropriate. So I remain only with the KL test. But ….. I wonder how this test is influenced by the number of data, the distribution, and how to interpret the result. For example I gather that if I have 2 KL tests, for each the minimum values in one case is 0.5 and the other case is 0.9, then in the first case the 2 datasets are less divergent than in the second case. But how small the KL should be until I can say that "one dataset is a good approximation of the other"? Granted, even that is debatable, because I can look at divergence itself, or how different one dataset is from the other in a statistical sense (this is one exercise) or look at the similarity / dissimilarity of data interpretation, or what the data "tell" about the variable measured. So …. It is enough that actually I look only at correlations in this case? Maybe I should classify my data in a number of convenient classes and compare one classification to the other through kappa-statistics? Maybe this will tell me at least if the interpretation of one dataset is similar enough with the interpretation of the second dataset? Do you have any other ideas? References? I would appreciate any thoughts on the matter. ## Edit Thanks for your answers. Now i realized i should have included few more details. Both raw lidar data go through a standard processing phase that results in a regular grid xyz of pre - determined resolution (in my case 5 m by 5 m). Both datasets have same xy coordinates and only z varies. I know that one sensor has usually an error of +/- 15 cm vertical, and probably the other one is very similar as well (i can find out actually). Both sensors were collecting data over the same area. We don't have any ground truth for altitude there. We want to compare these processed data and not the actual xyz lidar cloud (in other words the DEM result .....). I hope this is now a little bit more clear. Thanks.


title: How should I characterize the "typical" variance of a large collection of variables?

text: For an experiment I'm analyzing, I took many measurements of many different values. (To make it concrete, let's pretend I measured the height of people from many different age groups -- I have the height of ten 20-year-olds, ten 21-year-olds, ten 22-year-olds, etc.) I'd like to demonstrate that the variance of all my readings, per variable, is pretty low. (In my example, I'd like to show that people of any given age are all about the same height. That is, say 20-year-olds are all close to 5'4" with low variance and 21-year-olds are all close to 5'6". The particular heights don't matter; just that people of any given age are similar in height to other people of that age. That's obviously not true, but let's pretend it is.) What's the best metric to report? I could, for example, take the arithmetic mean of the variances -- or maybe the harmonic mean or geometric mean would be better? Or perhaps I should be averaging the standard deviations or the standard errors?


title: Length normalization in a naive Bayes classifier for documents

text: I'm trying to implement a naive Bayes classifier to classify documents that are essentially sets (as opposed to bags) of features, i.e. each feature contains a set of unique features, each of which can appear at most once in the document. For example, you can think of the features as unique keywords for documents. I've closely followed the Rennie, et. al. paper at http://www.aaai.org/Papers/ICML/2003/ICML03-081.pdf, but I am running into a problem that doesn't seem to be addressed. Namely, classifying short documents are resulting in much higher posterior probabilities due to the documents having a smaller number of features; vice versa for long documents. This is because the posterior probabilities are defined as (ignoring the denominator): $$ P(class|document) = P(class) * P(document|class) $$ which expands to $$ P(class|document) = P(class) * P(feature_1|class) * ... * P(feature_k|class) $$ From that, it's clear that short documents with fewer features will have higher posterior probabilities simply because there are fewer terms to multiply together. For example, suppose the features "foo", "bar", and "baz" all show up in positive training observations. Then, a document with single feature "foo" will have a higher posterior probability of being classified in the positive class than a document with features {"foo", "bar", "baz"}. This seems counter- intuitive, but I'm not quite sure how to solve this. Is there some sort of length normalization that can be done? One idea is to add the size of the document as a feature, but that doesn't seem quite right since results would then be skewed by the size of documents in the training data.


title: Exploiting features in a multiarmed bandit scenario

text: I am facing a challenging problem: Say I have shirts of three different colors (same price). And say I am running a strange kind of store in into which people come in one by one, and I can show them only one shirt, and they decide weather to buy or not to buy before they leave. I wish to optimize my sells. Up to this point, this looks to me like a classic multiarmed bandit problem, which would be classified as reinforcement-learning problem. Now say that I have a number of features about each person that comes into my store: age/gender and so on (a fairly large mix of categorical and non- categorical measures). How can I use this information to optimize my sells? I though about taking the following approach: Assuming that the expected rate of purchase for each of the colors is low, and though NOT equal, it is about the same (say around 1%), I can look at this problem as a classification problem (supervised learning). I will forget about the people who did not buy a shirt, and ask, from those who did buy a shirt, to which of the three groups of buyers does this new person resembles most? I`ll use a classification algorithm and try to show him the colored shirt that my trained classifier predicted. Even though I do not know weather the buyers would have bought other colored shirt had they been presented, my hope is that with large number of buyers, and hopefully relevant features, a classifier can find similarities between the groups of buyers. I know that this is far from a classic classification problem, but I think this can do better than just ignoring the feature I know of the people, and just use a multiarmed bandit like approach. My questions: Is my approach reasonable? Any different ideas to approach this? I would appreciate any though. Thanks!


title: Analyze and generate "clumpy" distributions?

text: Are there standard ways of analyzing and generating "clumpy" distributions?   * analyze: how clumpy is a given point cloud (in 1d, 2d, nd), what are its clumpy coefficients?   * generate or synthesize a pseudo-random cloud with coefficients C (These are the basics for any family of distributions, e.g. normal.) There are many kinds of clumpiness in nature (traffic jams, clumpiness of climate changes), so it's a wide term with I imagine various attempts at description and various links to "classical" statistics.   I'm looking for an overview; pictures would be nice. Added Friday 22 Oct: I had hoped to find methods for _both_ analysis and synthesis,   abstract <-> real both ways; surely a wide-ranging abstraction must do both. Still looking ... (Experts please add tags).


title: How do I get people to take better care of data?

text: My workplace has employees from a very wide range of disciplines, so we generate data in lots of different forms. Consequently, each team has developed its own system for storing data. Some use Access or SQL databases; some teams (to my horror) are reliant almost entirely on Excel spreadsheets. Often, data formats change from project to project. In some cases calling it a 'system' is too kind. The problems this entails are that I have to write new code to clean the data for every project, which is expensive; people manually editing spreadsheets make reproducibility and auditing of data near impossible; and even worse, there is a chance that data gets lost or made incorrect. I've been given the opportunity to discuss these problems with a board member of the company and I need to work out what to tell him. I think I've already persuaded him that we have a problem and that getting this right will enable better science and saving money. The question is: what should we be aiming for, and how do we get there? More specifically: How should we store data, in way that lets us track it from creation to publishing in a paper? (Databases stored on a central server?) How do you go about standardising database formats? Are there any good resources for educating people on how to care for data? (As a general rule, occupational hygienists and explosives engineers aren't data nerds; so non-technical content preferred.)


title: Identification of significant variables where dependent variable is categorical

text: I have a set of some 50,000 data points. There is one dependent variable which is categorical in nature and there are some 100 possible explanatory variables. Out of these 100, I have to select some 10 variables which are most significant in terms of determination/explanation of dependent variable. Can you please suggest which is the best method to do that?


title: Least stupid way to forecast a short multivariate time series

text: I need to forecast the following 4 variables for the 29th unit of time. I have roughly 2 years worth of historical data, where 1 and 14 and 27 are all the same period (or time of year). In the end, I am doing a Oaxaca-Blinder style decomposition on $W$, $wd$, $wc$, and $p$.               time    W               wd              wc               p     1       4.920725        4.684342        4.065288        .5962985     2       4.956172        4.73998         4.092179        .6151785     3       4.85532         4.725982        4.002519        .6028712     4       4.754887        4.674568        3.988028        .5943888     5       4.862039        4.758899        4.045568        .5925704     6       5.039032        4.791101        4.071131        .590314     7       4.612594        4.656253        4.136271        .529247     8       4.722339        4.631588        3.994956        .5801989     9       4.679251        4.647347        3.954906        .5832723     10      4.736177        4.679152        3.974465        .5843731     11      4.738954        4.759482        4.037036        .5868722     12      4.571325        4.707446        4.110281        .556147     13      4.883891        4.750031        4.168203        .602057     14      4.652408        4.703114        4.042872        .6059471     15      4.677363        4.744875        4.232081        .5672519     16      4.695732        4.614248        3.998735        .5838578     17      4.633575        4.6025          3.943488        .5914644     18      4.61025         4.67733         4.066427        .548952     19      4.678374        4.741046        4.060458        .5416393     20      4.48309         4.609238        4.000201        .5372143     21      4.477549        4.583907        3.94821         .5515663     22      4.555191        4.627404        3.93675         .5542806     23      4.508585        4.595927        3.881685        .5572687     24      4.467037        4.619762        3.909551        .5645944     25      4.326283        4.544351        3.877583        .5738906     26      4.672741        4.599463        3.953772        .5769604     27      4.53551         4.506167        3.808779        .5831352     28      4.528004        4.622972        3.90481         .5968299      I believe that $W$ can be approximated by $p\cdot wd + (1 - p)\cdot wc$ plus measurement error, but you can see that $W$ always considerably exceeds that quantity because of waste, approximation error, or theft. Here are my 2 questions.   1. My first thought was to try vector autoregression on these variables with 1 lag and an exogenous time and period variable, but that seems like a bad idea given how little data I have. Are there any time-series methods that (1) perform better in the face of "micro-numerosity" and (2) would be able to exploit the link between the variables?    2. On the other hand, the moduli of the eigenvalues for the VAR are all less than 1, so I don't think I need to worry about non-stationarity (though the Dickey-Fuller test suggest otherwise). The predictions seem mostly in line with projections from a flexible univariate model with a time trend, except for $W$ and $p$, which are lower. The coefficients on the lags seem mostly reasonable, though they are insignificant for the most part. The linear trend coefficient is significant, as are some of the period dummies. Still, are there any theoretical reasons to prefer this simpler approach over the VAR model? Full disclosure: I asked a similar question on Statalist with no response.


title: R package for visualizing and exploring large datasets

text: I recently read about a package, which I believe had a name like "tabstat" or "tablestats" which produced a really useful plot of univariate distributions of multiple variables for large datasets. I believe this package was developed by google, but am not sure. I can't seem to find it anywhere (tabstat on google returns a lot of results). If I'm completely crazy, and this package doesn't exist, what are some good packages for explore really huge datasets? I know plyr has some methods for binning 2D scatterplots to reduce overplotting, but I don't know of any other good packages for visually summarize huge datasets.


title: How to increase size of label fonts in barplot

text: I have the following _R_ code for creating _barplot_ :               # ... here I read cvs file          mx <- rbind(results$"AVG.P10.")        colnames(mx) <- results$"RUN"         rownames(mx) <- "AVG P"         postscript(file="avg_p_result.ps")         barplot(mx, beside=T,  col=c("grey"), names.arg= results$"RUN",              cex.axis = 1.5, cex.lab=1.5)      I have tried `cex.lab=1.5` but it does not work at all.


title: Learning a mapping from one time series to another with a Kalman Filter

text: I am interested in finding the relation between two (possibly multi dimensional) time series $x_{1:T}$ and $y_{1:T}$. I wonder how I can do that with a linear dynamical system/Kalman filter. My approach would be sth like estimating the parameters $\theta$ with EM for the joint probability $p(x_{1:T}, y_{1:T}, h_{1:T})$ where $h_{1:T}$ is my hidden state. Then say that I have observed a sequence $x'_{1:T}$, I want to determine the most likely sequence $y'_{1:T}$ by finding $p(y'_{1:T}|x'_{1:T})$. Does anyone know whether this is the standard approach for this problem? If not, what is the standard approach? Does anyone know a source where this is described?


title: What are alternatives to uniform distribution when trying to fit observed data distribution?

text: I have a dataset of 500 observations and I have to find the best fitting distribution. If I look at the histogram of the data and the QQ plot my data seems to follow an uniform distribution. But a chi-square test rejects the corresponding null hypothesis, and I can't use Kolmogorov-Smirnov test (which doesn't reject null hypothesis) because parameters are estimated from the data. What other distributions could I try to fit?


title: Estimating a distribution from a dataset with multiple parameters

text: How would you go about solving the following problem? You're an insurance company who writes workers compensation policies. You want to build a probability distribution for the number of annual claims of an account based on certain characteristics. You determine that the three most important predictor variables are the number of employees, the assigned risk level ("high", "medium", "low"), and the primary state where the business is located ("AL", "AK", ...). You think a poisson distribution would be a good model, but if you try to build a separate model using only the data in each unique combination of the variable classes, you don't have enough data to be credible. For instance, if you try to fit a distribution to AL, high risk accounts with 10 employees, you only have 1 or 2 observations to use. How do people typically solve this sort of problem. I'm sure its been done a thousand times before, but I don't have much experience with it. Perhaps some sort of glm? What makes it more confusing to me is how you might combine the three very different types of predictors - a numeric predictor, and ordinal factor, and a non ordered factor.


title: Trick in survey design for "problematic" questions

text: Years ago I read that there's a certain trick researches resort to when they want to obtain survey data regarding questions where one wouldn't expect people to answer honestly, such as "Do you beat your children?". The trick had something to do with randomization and mixing the contentious questions with innocent questions and then doing some sort of analysis. Unfortunately I can't remember what the method was called and how exactly it works. Does anyone have a pointer for me?


title: College Football Dataset

text: Several years ago I read about a contest to generate an algorithm to predict the outcome of a college football game. The contest provided a comprehensive data set of games from previous years. I didn't get the data set at the time and as a college football fan and amateur statistician I wish I had. Is anybody aware of a data set or source like this?


title: How many possibilities is there to choose 10 boxes of given size?

text: We have a collection of 40 boxes, each is one of four sizes (evenly distributed 10 of each): `sizeA`, `sizeB`, `sizeC`, `sizeD`. Our person is presented with 1 box of each size, they must select one box from the different sizes, then the remaining three will be taken away and four new boxes of each size is presented. I know this gives me $4^{10}$ possibilities. But consider this now. Of the 10 boxes the person must end with these sizes:               1 sizeA,      2 sizeB,      3 sizeC,      1 sizeD,      3 of any size.      I'm trying to figure out based on this how many possibilities I would have now? I'm just a web developer trying a side project thing and I'm very fresh with this.


title: What makes values from two features equal for the purposes of constructing a maximum entropy model?

text: I'm looking to construct a maximum entropy model using general iterative scaling. In the process, I've come across the SharpEntropy project for .NET and I'm using that as a reference for the implementation I'm building. That said, the simple example uses the following observations (all nominal values):               Day       Temp.       Moisture  Time of Day   Take Umbrella (outcome)     ---       -----       --------  -----------   -----------------------       1       Warm        Dry                     No_Umbrella       2       Cold        Dry                     No_Umbrella       3       Cold        Rainy                   Umbrella       4       Cold        Dry                     Umbrella       5       Warm        Dry                     No_Umbrella       6       Cold        Dry       Early         Umbrella       7       Cold        Rainy     Early         Umbrella       8       Cold        Dry       Late          No_Umbrella       9       Warm        Rainy     Late          No_Umbrella      10       Warm        Dry       Late          No_Umbrella      Looking through the code, it seems that in order to generate the binary functions (predicates) for the maximum entropy model, it's taking the labels across _all the attributes_ , no matter where the value is found. Meaning, if `Warm` was in the `Moisture` column, then that would be tallied when figuring out the entropy associated with the binary function `observation of warm`. This seems _very_ incorrect to me. Granted, for the data above, it _happens_ to work because all of the labels are unique across all of the attributes. Let's add another feature to the set above, the temperature of the prior day:               Day   Prior Temp.  Temp.  Moisture  Time of Day  Take Umbrella (outcome)     ---   -----------  -----  --------  -----------  -----------------------       1   Warm         Warm   Dry                    No_Umbrella       2   Cold         Cold   Dry                    No_Umbrella       3   Cold         Cold   Rainy                  Umbrella       4   Warm         Cold   Dry                    Umbrella       5   Warm         Warm   Dry                    No_Umbrella       6   Cold         Cold   Dry       Early        Umbrella       7   Cold         Cold   Rainy     Early        Umbrella       8   Warm         Cold   Dry       Late         No_Umbrella       9   Cold         Warm   Rainy     Late         No_Umbrella      10   Warm         Warm   Dry       Late         No_Umbrella      Now what SharpEntropy seems to do is generate the following binary function:               // I wanted to do this in LaTeX, but my LaTeX-fu failed me.     // This would be, f0, for example.     return (priorTemp == "Warm" || priorTemp == "Warm");      This seems wrong to me, and that you really want _two_ binary functions, one for the prior day's temperature:               // f0     return priorTemp == "Warm";      And then one for the current day:               // f1     return temp == "Warm";      That said, is this the proper way to generate the binary functions, that each attribute on an observation is has it's own universe of values, even though they might semantically mean the same thing when considering a larger universe of values (say, temperature in general versus temperature on a given day)?


title: ARIMA identification process

text: I have a 2000 SKUs, I want to identify the ARIMA process for them, what is the best way to do that?, which tools should I use? thank you very much. Roji


title: Specifying the Form of Prior, Likelihood and Posterior Distributions for Bayesian Analysis

text: I have recently begun to look into Bayesian Analysis, and, although I'm beginning to get to grips with the general framework (i.e. $\text{posterior} \propto \text{likelihood} \times \text{prior}$), I'm still having some difficulties in playing around with distributions. I hope to gain a more concrete understanding of these processes through worked examples; what follows is a template-like rendering of the kinds of problems I've seen in various books and online resources - I hope it makes sense. Assume I want to carry out a Bayesian analysis on some process. Data ($x$) gathered on this process over a certain time period suggests that the process is normally distributed with variance of $\sigma$ (known); however, the goal of the research is to discover the value of the mean ($\mu$) (unknown). Another researcher is asked his/her opinion on the same process and suggests a mean value of $\mu_{0}$ with standard deviation of $\sigma_{0}$.   * Now, assume I wanted to specify a _gamma_ prior distribution from the fellow researcher's opinion. What would the pdf look like in this instance? (i.e. How do I embody the mean and variance into the gamma distribution?)   * Furthermore, for completeness, what would the likelihood and posterior pdfs of $\mu$ look like? As you can see, it's mostly the manipulation of the distributions that are causing my difficulties.


title: Covariance matrix of uniform spherical distribution

text: I need to figure out the covariance matrix of a uniform spherical distribution. But there I can't even find a closed form of the distribution. This link says it is $\frac{1}{n}\mathbf{I}$, where $\mathbf{I}$ is the $n \times n$ identity matrix. But can someone please help me derive it? Or even point to literature on how to? Thanks! p.s. By uniform spherical distribution, yes, I meant the unit sphere on $\mathbb{R}^n$, that is $S^{n-1} = \\{ s \in \mathbb{R}^n : |s| = 1 \\}$.


title: Propagation of uncertainty/standard deviations in biological experiments

text: I am trying to work out how to propagate standard deviations in a biological experiment, but I have some difficulties. I have the following (semi)fictitious data originating from an image analysis resulting in optical density (OD) measurements:               ds <- data.frame(         Sample = c(rep(seq(1,6),2)),         Control.OD = c(7000, 7100, 5600, 5200, 5900, 7000,                         8100, 7700, 6100, 5500, 6600, 7500),          Target.OD  = c(1000,  330,   35, 9300, 5570, 8700,                         1300,  400,   62, 9100, 5817, 9000),         Group = c(1,1,1,2,2,2,1,1,1,2,2,2)     )          print(ds)             Sample Control.OD Target.OD Group     1       1       7000      1000     1     2       2       7100       330     1     3       3       5600        35     1     4       4       5200      9300     2     5       5       5900      5570     2     6       6       7000      8700     2     7       1       8100      1300     1     8       2       7700       400     1     9       3       6100        62     1     10      4       5500      9100     2     11      5       6600      5817     2     12      6       7500      9000     2      The data represents duplicate measurements from one experiment with 6 biological samples from two different groups. I do the following:               # Aggregate technical duplicates     ds.Agg <- aggregate(cbind(Control.OD, Target.OD) ~ Sample + Group,                          data = ds, mean)          # Rename columns for easier reading     names(ds.Agg)[names(ds.Agg)=="Control.OD"] <- "Control.mean"     names(ds.Agg)[names(ds.Agg)=="Target.OD"] <- "Target.mean"          # Calculate initial ratio between the target and the control     # This is the value that represents each biological sample     ds.Agg$Ratio <- ds.Agg$Target.mean / ds.Agg$Control.mean           # Calculate the group means and normal sd (no propagation)     ds.group <- aggregate(Ratio ~ Group, data = ds.Agg, mean)     ds.group$sd <- aggregate(Ratio ~ Group, data = ds.Agg, sd)$Ratio     ds.group$sd.rel <- ds.group$sd / ds.group$Ratio * 100          # Print results for classical way of doing this (no propagation)     format(ds.group, digits=3)            Group Ratio     sd sd.rel     1     1  0.07 0.0742  106.0     2     2  1.28 0.4080   31.8          A final read-out would be:               library(ggplot2) # version 0.9.0             png(type = "windows", file = "c:/temp/propagation1.png",width=140,height=200)     ggplot(ds.group, mapping = aes(x = factor(Group), y = Ratio)) +       geom_bar(width = 0.8, position="dodge") +       geom_errorbar(mapping = aes(ymin = Ratio - sd, ymax = Ratio + sd),                      width = 0.40, position=position_dodge(width = 0.8)) +       scale_x_discrete("Group") +       scale_y_continuous(limits = c(-0.1, 1.7)) +       theme_bw()     message(dev.off())      ![No propagation](http://i.stack.imgur.com/RLQzm.png) Then an attempt on propagation:               # Calculate sd for technical duplicates which will be propagated     sdTmp <- aggregate(cbind(Control.OD, Target.OD) ~ Sample + Group, data=ds, sd)      ds.Agg$Control.sd <- sdTmp$Control.OD     ds.Agg$Target.sd <- sdTmp$Target.OD          # Calculate the relative propagated sd associated with the calculated Ratio     fltRelativeSd <- with(ds.Agg,          sqrt( (Control.sd / Control.mean)^2 + (Target.sd / Target.mean)^2 )      )     ds.Agg$Ratio.sd <- ds.Agg$Ratio * fltRelativeSd     ds.Agg$Ratio.sd.rel <- ds.Agg$Ratio.sd / ds.Agg$Ratio * 100      **The follow code is my biggest concern when it comes to calculation of the propagated standard deviation:**               # Calculate propagated sd for the groups     ds.group$sd.prop <- aggregate(Ratio.sd ~ Group, data = ds.Agg,          function(x) {              sqrt( sum((x)^2) ) / length(x);         })$Ratio      **to here.**               ds.group$sd.prop.rel <- ds.group$sd.prop / ds.group$Ratio * 100      Re-print the results with the propagated standard deviations:               format(ds.group, digits=3)            Group Ratio     sd sd.rel sd.prop sd.prop.rel     1     1  0.07 0.0742  106.0  0.0111       15.79     2     2  1.28 0.4080   31.8  0.0418        3.26      And a graph:               png(type = "windows", file = "c:/temp/propagation2.png",width=140,height=200)     ggplot(ds.group, mapping = aes(x = factor(Group), y = Ratio)) +       geom_bar(width = 0.8, position="dodge") +       geom_errorbar(mapping = aes(ymin = Ratio - sd.prop, ymax = Ratio + sd.prop),                      width = 0.40, position=position_dodge(width = 0.8)) +       scale_x_discrete("Group") +       scale_y_continuous(limits = c(-0.1, 1.7)) +       theme_bw()     message(dev.off())      So a read-out for the propagation: ![With propagation](http://i.stack.imgur.com/1vFrk.png) I would have expected standard deviations in the same neighborhood as seen on the first graph, however they seem a factor 10 off. So my thoughts/questions at this point:   1. Am I doing the calculations right?   2. Does it make sense to do this calculation (it seems that the biological variance is ignored and only the initial duplicate sd is actually propagated and presented)?  Any comments are most welcome. Additional info: The experiment uses a technique called Western Blotting and the image analysis was done in ImageJ in a similar fashion as described here.


title: Inequality for bivariate normal distribution

text: Let $X_1$ and $X_2$ be bivariate normal with mean $\mu=(0,\mu_2)$, for any $\mu_2$, and correlation $\rho$. Consider the following inequality: \begin{align*} Pr\left\\{|X_1| \ge \Phi^{-1}(1-\alpha/2), |X_2| \ge \Phi^{-1}(1-\alpha/2)\right\\} \le 2Pr\left\\{|X_1| \ge \Phi^{-1}(1-\alpha/4), |X_2| \ge \Phi^{-1}(1-\alpha/2)\right\\} \end{align*} Is the inequality true for correlation $|\rho| \le 1$? Are there relevant results on general bivariate normal distributions?


title: Test for significance with no comparator

text: I have roughly 15,000 non-independent features for **each** of which I have a sample size of 100. For each of these feature-sample data points, I have a value of TRUE or FALSE corresponding to some status of the feature, such that for each feature, I have 100 TRUE or FALSE values. How can I test the "trueness" of each feature and thereby determine the percentage of features with a significant TRUE value? A suitable null hypothesis may be that the feature does not have a TRUE status. My problem though is that an $\alpha$ of 0.05 corresponding to a rejection of the null hypothesis if there are >95 TRUES for a feature, seems too high. I have no reference dataset by which to make a comparison. ![enter image description here](http://i.stack.imgur.com/Fa2qj.jpg) * * * UPDATE: Re: Use of Fisher's exact test. I am not sure what parameters to use for the 2 x 2 contingency table using Fisher's exact test or chi-squared test. What would class 1, class 2, sample 1 and sample 2 be? 'TRUE count', etc.? Here are the first few rows and columns of the data:                         Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 Sample7 Sample8 Sample9 Sample10     Feature1    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature2    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature3    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature4    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature5    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature6    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature7    TRUE    FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature8    FALSE   FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature9    FALSE   FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature10   FALSE   FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   FALSE     Feature11   FALSE   FALSE   TRUE    FALSE   TRUE    TRUE    FALSE   TRUE    FALSE   TRUE     Feature12   FALSE   FALSE   TRUE    FALSE   TRUE    TRUE    TRUE    TRUE    FALSE   TRUE     Feature13   FALSE   TRUE    TRUE    FALSE   FALSE   TRUE    TRUE    TRUE    FALSE   TRUE     Feature14   FALSE   TRUE    TRUE    FALSE   FALSE   TRUE    TRUE    TRUE    FALSE   TRUE     Feature15   FALSE   TRUE    FALSE   FALSE   FALSE   TRUE    TRUE    TRUE    FALSE   TRUE


title: Appropriate residual degrees of freedom after dropping terms from a model

text: I am reflecting on the discussion around this question and particularly Frank Harrell's comment that the estimate for variance in a reduced model (ie one from which a number of explanatory variables have been tested and rejected) should use Ye's Generalized Degrees of Freedom. Professor Harrell points out this will be much closer to the residual degrees of freedom of the original "full" model (with all the variables in) than that from a final model (from which a number of variables have been rejected). Question 1. If I want to use an appropriate approach to all the standard summaries and statistics from a reduced model (but short of a full implementation of Generalized Degrees of Freedom), would a reasonable approach be to just use the residual degrees of freedom from the full model in my estimates of residual variance, etc? Question 2. If the above is true and I want to do it in `R`, might it be as simple as setting               finalModel$df.residual <- fullModel$df.residual      at some point in the model fitting exercise, where finalModel and fullModel were created with lm() or a similar function. After which functions such as summary() and confint() seem to work with the desired df.residual, albeit returning an error message that someone has clearly mucked around with the finalModel object. Thanks.


title: Machine learning self-learning book?

text: > **Possible Duplicate:**   >  Machine learning cookbook / reference card / cheatsheet? I wonder if there is a good self-learning textbook for machine learning? I am particularly looking for those in application-level, not so much theoretical. I wouldn't mind a math-heavy book, as long as you guys think it is suitable for self-learning. Thanks!


title: Appropriate problems for different type of ML classifiers

text: I have learned several classifiers in Machine learning - Decision tree, Neural network, SVM, Bayesian classifier, K-NN, Markov process...etc. Can anyone please help to understand when I should prefer one of the classifier over other - for example - in which situation(nature of data sets, etc) I should prefer decision tree over neural net OR which situation SVM might work better than Bayesian OR what types of problems are appropriate to apply decision tree, neural net, k-nn, markov processs or SVM or bayesian ?? Thanks.


title: General rules for choosing machine-learning algorithm?

text: > **Possible Duplicate:**   >  Machine learning cookbook / reference card / cheatsheet? There are numerous machine-learning approaches out there. Also there are numerous ways how to optimize their parameters and numerous ways to preprocess your data beforehand. It seems to be a common practice to initially implement some off-the-shelf approach (SVM, random forest) and then trying to beat the accuracy with more sophisticated approaches. I would like to know if there are any guidelines, recommendations or rules of thumb which approach to use in which context? Are recommendations data-specific or domain-specific (e.g. genetics, econometrics)? Also I would like to know how to account this multitude when writing results up. Trying out multiple approaches and only reporting the best one in your publication cannot be valid, right?


title: Overview of (common) classifiers and their characteristics

text: > **Possible Duplicate:**   >  Machine learning cookbook / reference card / cheatsheet? Each classifier has it's own advantages and disadvantages. E.g. train/test speed, classification/regression (and how many classes can be handled), how many degrees of freedom, suitable/not suitable for abstract kinds of datasets (of possible to say), how interpretable is the resulting model (like for max-ent: importance of single features), and so on and so on Do you know a good overview? If not, would you be in to create one?


title: Learning Algorithms/Neural networks books?

text: > **Possible Duplicate:**   >  Machine learning cookbook / reference card / cheatsheet?   >  Machine learning self-learning book? What are some good theoretical and practical machine learning algorithms books? I am particularly looking for a book that explains the most important and prominent machine learning algorithms being used today. Also is there a good book for learning basic neural networks? I am looking for both theory and application so suggestions on either is good. I am only a beginner in the field, so the easier to understand the book is the better!


title: Machine learning cookbook / reference card / cheatsheet?

text: I find resources like the Probability and Statistics Cookbook and The R Reference Card for Data Mining incredibly useful. They obviously serve well as references but also help me to organize my thoughts on a subject and get the lay of the land. Q: Does anything like these resources exist for machine learning methods? I'm imagining a reference card which for each ML method would include:   * General properties   * When the method works well   * When the method does poorly   * From which or to which other methods the method generalizes. Has it been mostly superseded?   * Seminal papers on the method   * Open problems associated with the method   * Computational intensity All these things can be found with some minimal digging through textbooks I'm sure. It would just be really convenient to have them on a few pages.


title: Estimating misclassification rate from summary classification values

text: I'm trying to work out the misclassification rate in a published study, where the results do not include this value. However, I have the `sensitivity`, `specificity`, `positive predictive value` and `negative predictive value` results. The results are:                Sensitivity     Specificity     PPV      NPV      69.9            99.3            95.6     93.7      The sample size is 860. I figure this can be solved using a set of 4 simultaneous equations, but my algebra isn't up to the task. If we use the following letters for replacement in the equations:   * a = true positives   * b = false positives   * c = false negatives   * d = true negatives Then the set of simultaneous equations appears to be:                a/(a+c)=.699 (i.e. sensitivity)      b/(b+d)=.993 (i.e. specificity)      a/(a+b)=.956 (i.e. PPV)      d/(c+d)=.937 (i.e. NPV)      The fractional terms in the equations have been driving me nuts. Surely someone else must have had the same problem in the past and solved it, but I have been unable to find a ready reckoner on how to solve this problem. Once I know what the four values are, I can easily calculate the misclassification rate. I would appreciate any assistance on this back calculation problem. I've posted here as I'm assuming an epidemiologist or biostatistician is the most likely type of person to answer and they may not read the mathematics questions.


title: Calculating $\tan 2\theta$ in terms of correlation coefficients and variances

text: I know I am asking a few questions and I apologize for that. However I shall try to read the solutions and know the loopholes in my theory. Here it is: If $u=x\cos \theta+y\sin\theta$ and $v=y \cos \theta-x\sin\theta$. The variables $u$ and $v$ are uncorrelated, then how may I prove that $$\tan 2\theta=\frac{2r_{xy}s_xs_y}{s_x^2-s_y^2}$$ where $r_{xy}$ is the correlation co-efficient between $x$ and $y$ and $s_i$ is the standard deviation of a variable $i$. Sorry for asking so many questions without much effort but it is late at night now and I need to identify my faults soon.


title: Latent variables, overparameterization and MCMC convergence in bayesian models

text: Sometimes I have a large number of latent variables in a Bayesian hierarchical model to which, but I am only interested in estimating projected transformations of those latent variables (for example, I will parameterize a binomial parameter as an inverse logit of a set of possibly-non-identifiable covariates, even though the result I'm interested in is the binomial parameter estimate). The projected transformations will often converge very quickly (based on convergence diagnostics such as the Gelman/Rubin or by eyeballing the posterior samples) even if the latent variables have not yet converged. Intuitively this makes sense, the model may be an overparameterization where the latent parameters are not identifiable - the derived quantities are constrained to be in a constrained a narrow high-likelihood region of the transformed variables' parameter space which maps to a much larger largely flat likelihood (but bounded) region of the latent variable parameter space. So is the intuition correct that I shouldn't be concerned that the overparameterized latent variables are not identifiable and aren't fully converged when I take my posterior samples? Are there some good references which discuss the use of non-identified latent variables in this way? I've heard some discussion on overparameterizing to speed up mcmc convergence, but I'm not entirely clear on how to think about this, as the approaches and attitudes towards overparameterization and non-identifiability in bayesian methods seems to be a bit different than in other areas of modeling.


title: How important is $R^2$ when prediction is not the research aim?

text: I have the problem that in my binomial-glmm model with a nested random effects structure (30,000 rows located 30 groups nested in 25 groups) the diagnostic plots look okay (as far as I can tell), residuals do not show heavy patterns and some predictors are even significant. However, the $R^2$ is 0.055. My goal is not to predict new values but only to understand the relationships in the present dataset. Is $R^2$ still important? Or is it possible that my structure is nonsense and thus leads to the low $R^2$? Additional info: there are 25 Wind energy stations (WEA) beeing build. During building these, noise levels are believed to have an impact on porpoises (dolphin like species). There 30 acoustic stations measuring porpoise noises (they can differentiate between ships and fish), temperature value and noise levels (SEL). the WEA sometimes use a noise-muting-device. The aim is to find out whether this device leads to more porpoise clickx during a certain time period. Measuremnts of the stations were done over a year, every minute. however, of interest is only the 100 hours interval after the WEA-build was finished. This is called the HRW-phase (hour related work). There is always one WEA-build per time, not all together. So WEA can be thought of as an event. Dependent variable:   * ppm (counted porpoise clicks per minute during an hour) so hourly values      * zero-inflated count data Independent variable:   * SEL a measure of noise, averaged per day (continous)   * average daily temperature (continous)   * total building time (continuos = minutes)   * locations x/y of the stations (still in geographic system, i.e. 54,333 9,3333)   * hrw (time phase after work, from 0,1,2,...,100. stations do not always measure the full lentgh ,sometimes only 17 hours)   * noise muting device (binary factor)   * month as a factor with 7 level (cat.)   * WEA ID (cat, random effect)   * acoustic station (cat, random effect)   * podid at one station, sometimes pod-acoutstic measure systems per station were exchanged (cat)   * potentially also 4 distance classes  I build the random effects with (hrw|WEA/station) although I am still not sure if this is the right design. the thing is that WEA is an event (there 25) at during an event 30 stations measure ppm. is this then crossed or nested? i am not sure. although I expected temporal and spatial autocorrelation I did not find any hints on that. I hope this clarifies the story a bit more.


title: Outlier treatment in Vector Autoregression (VAR) Model

text: Data: Multivariate Time Series, Series 1) Demand of a product 2) Rainfall data both available at monthly level from 2010-2013. Approach: I am trying to estimate the effect of rainfall on demand of the product using VAR( Vector Autoregression) model. Demand data has some outliers, like a month of sudden high demand and followed by zero values. **Question:** How to treat these outliers (I am working in R), since I already have few data and deleting them is not an option for me.


title: What is discriminative training?

text: I was reading M. Collins' paper introducing structured perceptron (here is a link). I do not understand what is meant by discriminative training of the HMM. How can you train a generative model in a discriminative manner? (and what exactly does that mean)


title: How to visualize price distribution in marketplaces?

text: I am working in a hobby project, and I'd like some help with the following problem: I'm pulling pricing information from an internet marketplace, so if I provide a keyword (say, iphone cases) I'm going to get the first 500 items that show up in that marketplace, along with their prices. My question is, how can I properly visualize price distribution, so that I can get meaningful information? I've tried scatter plots, but I don't think the graph tell me much about the prices in that marketplace: ![my attempt to visualize price distribution](http://i.stack.imgur.com/Y3W10.png) I've also tried to plot the mean, and the standard deviation in the scatter plot, but I'm not sure if I'm mixing pears and apples, when I plot those values in a scatter chart. I'm not a data miner, nor a stats person. I'm doing this just as a hobby and using this small project as a personal learning experience. Any help would be greatly appreciated.


title: Hypothesis testing: am I interpreting the results from the OLS regression correctly?

text: Consider the following (partial) regression results, in which the explainatory value of 5 independent variables are measured for 3 dependent variables. In this post, I'm referring to the models `(1d)`, `(2d)` and `(3d)`. For each model, the first column is the coefficient, the second the p-value: ![enter image description here](http://i.stack.imgur.com/AP1O1.png) In my thesis, I defined several hypotheses, a.o.   1. if `hvVol0LN` increases, `trading volume` increases   2. if `hbVol0LN` increases, `returns` increase   3. if `hbAgreeQ0` decreases, `returns` increase   4. if `svi0` increases, `trading volume` increases It is now time to check whether these assumptions were right. I just want to double check whether I'm interpreting these results right:   1. hypothesis accepted: there is a relatively strong (`r=.146`) and statistically significant (`p=.000`) relationship   2. hypothesis rejected; the relationship is very small to begin with, but what is more it is statistically insignificant (`p=.142`). It would only be significant at the 14,2% significance level   3. hypothesis accepted at 5% level: the relationship between `hbAgreeQ0` and `returns` is indeed negative and relatively strong (`r=-.100`), and statistically significant at the 5% level (`p=.005`). If `p` would have been larger than `.005`, the relationship would not be statistically significant, since in practive significance levels of 1% to 5% are accepted (depending on the type of research, of course)   4. hypothesis accepted: the relationship is significant (`p=.000`) but not very strong ('r=.052') Is this a correct way of describing the outcome of the hpyotheses? In other words:   * can you derive the acceptance / rejection directly from the regression output, and   * am I interpreting the `r` and `p`-values correctly? Your help in finally finishing my thesis is greatly appreciated :-) **edit: some SPSS regression output to be more clear of what I mean by`r`** ![enter image description here](http://i.stack.imgur.com/LzhBI.png)


title: How to test association in contingency tables with very small numbers and proportions?

text: I have a sample of 197 responses. 8.6% (17)are from Group A, the remainder from Group B. (The groups are mutually exclusive, and not independent). (If it helps visualize the issue, Group A is under 10's, Group B respondents aged 11+.) The responses are then sorted in to groups by type. Type 1 is "head injuries", and there are 3 members of this group. 2 are from Group B, 1 is from Group A. I would expect the results for head injuries to be under 10 (Group A) & head injury: 0.3 Under 10 & no head injury: 16.7 Over 10 (Group B) & head injury: 2.7 Over 10 & no head injury: 179.3 I have asked a few friends how to test if the expected value being so much lower than the observed value is significant. So far responses seem to favour a Z test, or a chi squared test with a correction for _very small numbers_. I'm innumerate in the extreme, but can manage either well enough -- just not sure what is the appropriate test here?


title: Flat maximum effect?

text: I've seen this referred to in relation to model selection but haven't come across a definitive description of it. What the flat maximum effect?


title: Recovering raw coefficients and variances from orthogonal polynomial regression

text: It seems that if I have a regression model such as $y_i \sim \beta_0 + \beta_1 x_i+\beta_2 x_i^2 +\beta_3 x_i^3$ I can either fit a raw polynomial and get unreliable results or fit an orthogonal polynomial and get coefficients that don't have a direct physical interpretation (e.g. I cannot use them to find the locations of the extrema on the original scale). Seems like I should be able to have the best of both worlds and be able to transform the fitted orthogonal coefficients and their variances back to the raw scale. I've taken a graduate course in applied linear regression (using Kutner, 5ed) and I looked through the polynomial regression chapter in Draper (3ed, referred to by Kutner) but found no discussion of how to do this. The help text for the `poly()` function in R does not. Nor have I found anything in my web searching, including here. Is reconstructing raw coefficients (and obtaining their variances) from coefficients fitted to an orthogonal polynomial...   1. impossible to do and I'm wasting my time.   2. maybe possible but not known how in the general case.   3. possible but not discussed because "who would want to?"   4. possible but not discussed because "it's obvious". If the answer is 3 or 4, I would be very grateful if someone would have the patience to explain how to do this or point to a source that does so. If it's 1 or 2, I'd still be curious to know what the obstacle is. Thank you very much for reading this, and I apologize in advance if I'm overlooking something obvious.


title: What is the difference between normal distribution and standard normal distribution?

text: Two terms normal distribution and standard normal distribution are used in statistics. Does standard term contribute to the normal distribution anything? Please give a simple-however a substantive reasoning in the back of these terms.


title: How can I estimate the likelihood of being in a traffic accident on the US101 by day of week and time of day?

text: My plan is this:   * Find source of data about when and where accidents occurred on the US 101   * Find a source of data about traffic volume on the same road   * Subset the data to include only accidents that occurred on the US101 between San Francisco and Palo Alto.   * Divide accidents by traffic volume for as small of time intervals as I can get traffic volume data. For example, if I can get traffic volume per hour, that would be great, because then I can divide average # of accidents in that hour on a given day by the volume of traffic in that window, and then assume, for lack of a better idea, that every car has an equal chance of being involved. Maybe I can get some data about different risk levels by driver age or type of car, but I imagine the insurance companies hold that data and aren't likely to share. Suggestions for sources of data are greatly appreciated. Even if I have to do something lousy like multiplying national risk per capita by traffic volume, that would be good enough for now, my main problem is getting the data. Any clever suggestions appreciated. (FYI: this is just for personal interest)


title: Dividing a uniform by a normal random variable: What's the distribution?

text: I have two random variables $s\sim \mathcal{N}(\nu,\sigma)$ and $a\sim \mathcal{U}(0,A)$, $0<A<1$, calculate a third r.v. $t=(1-a)/s$ and want to find its distribution $p(t|\nu,\sigma,A)$. My reasoning is, that for any value of $a$, there is exactly one $s=(1-a)/t$ such that the $(a,s)$ pair will produce $t$ and I will therefore have to integrate over the probabilities for these values: $$p(t|\nu,\sigma,A) = \int_0^A p_a(x)p_s((1-x)/t)\,\mathrm dx\\\$$ which can be expressed as the sum of two error functions. But already when checking this step numerically, I get a discrepancy between estimated (by sampling) and calculated distribution: ![](http://i45.tinypic.com/2da0ew5.png) _Why do the histogram and the calculated distribution do not match?_ This is the code I used:               n=10000000          A=0.7     nu=.7     sigma=.1          # sampling from the target distribution     s=rnorm( n, mean=nu, sd=sigma )     a=runif( n, min=0, max=A )     t=(1-a)/s     hist(t,200, freq=F, xlim=c(0,5), ylim=c(0,2.5))               # plot the analytical result     analytic <- function(t, A, nu, sigma ){       tmp <- function(x, A, t, nu, sigma){          return (1/A*dnorm( (1-x)/t, mean=nu, sd=sigma) )       }       return( integrate( tmp, 0, A, A=A, t=t, nu=nu, sigma=sigma)$value)     }          x=seq(0,5,by=.01)     y=rep(0,length(x))     for( i in seq(1,length(x)) ){       y[i]=analytic(x[i], A, nu, sigma)     }     lines( x, y, col="red", type="l")


title: Expectation of Sums of Squares

text: $\mu$=mean, $\sigma$=standard dev show that if $H_0:\mu_1=...=\mu_k$ in the normal one way classification then the Total Sum of Squares has expectation $E(T)=(N-1)\sigma^2$ the Within Groups Mean Square has expectation $E(W/(N-k))=\sigma^2$ and hence that the Between Groups Mean Square has expectation $E(B/(k-1))=\sigma^2$ I understand the definitions of each of them but i am completely unsure about how to go about getting the expectation. I assume then i will just be algebraic manipulation.


title: Mathematical formulation SVM Model

text: Good morning, For a homework I used a support vector machines (classification) with a RBF kernel, k-fold cross validation=10, cost of constraints violation=100, and gamma=0.001. It works great and my predictions are very good. However, we didn't see how to formulate the model in mathematical terms (well we are not going to cover svm in the course at all), and I believe it's important. I reviewed Hastie's and other books, but I still don't know how to formulate this model. Can somebody point me in the right direction? Thanks, Diego


title: Comparing multiple time series in R

text: If I have many time series that I'd like to compare to see if there are relationships between the variables, (I have several dependent variables and many more independent variables) how might I go about doing this (I'm working in R, just fyi)? I haven't really found too many examples seeking to compare and explore the relationships between many variables. In particular, I'd like to see if the variation in my independent time series drives variation in the dependent time series, and I'm so new to statistics (and R) that I'm really not sure how to approach this problem. Here is some sample data (I know there are 2 missing values, and I may choose to treat Y2 as an independent variable since I think Y2 and Y1 may be correlated. The column labeled "covariate" is years since last flood, because I think this might be important too):               year Y1 Y2  X1  X2  X3      X4  X5  covariate     1    40  92 0   0   20.6    91  503 3     2    54  65 0   0   21.7    33  175 4     3    59  75 1   1   22.2    34  94  5     4    68  53 8   9   22.2    24  86  6     5           5   20  20.6    5   185 7     6    76  65 8   13  22.2    32  119 8     7    76  55 16  18  23.3    0   153 9     8    82  58 18  2   24.4    19  0   1     9    60  57 28  24  23.33   0   223 2     10   58  46 18  3   22.78   0   184 3     11   49  48 2   1   23.33   0   110 4     12   28  76 0   3   22.78   0   213 5     13   56  61 0   1   22.78   12  123 6     14   105 53 56  24  23.33   0   122 7     15   99  43 28  13  24.44   0   154 8     16   119 47 46  35  23.33   0   182 9      Any guidance would be much appreciated as I'm not really sure how to proceed here.


title: To find significant differences between means using repeated measures GLM

text: I am new in using SPSS and stuck with General Linear Model [which is just linear model] procedure. I have 8 cognition variables like latency, thigmotaxis, spatial, quadrant, nt quadrant, proximity, proximity nt and active speed. All these were analysed on day 1,2,3,4 and 5 And were tested for four different diet groups and only one gender (female). My boss wants to find the significant differences in means between the days and among diets and represent it as a table using GLM. I tried with just one variable using the GLM repeated measures but I don't know which p value must be reported. I know so far that we have to go Analyse> general linear model> repeated measures> within-subject (entered as time)> levels (entered as 5)> define> within-subject variables (chosen for all days)> between-subjects (entered diet)> contrasts(repeated)> options(emmeans)> ok. **I still don't understand how p values were obtained as shown in the image below**. I didn't see a table that shows like this or something similar to this. I don't know what I am missing here. My data is in Means +/- SEM and N. ![enter image description here](http://i.stack.imgur.com/XIBEq.png)


title: Performance comparison between reduction techniques

text: If I want to compare the performance of PCA and PLS with nonlinear data, what do I have to measure ? I know that I can use the MSE (but if I do not use a linear profile, I can't estimate the parameters of the linear regression, so how do I calculate it? For the PLS, the data should follow a normal ditribution or it is an unsupervised method. Which other methods would you propose to me for dimension reduction techniques?


title: How to calculate Eigenvalue without using eign() function with R?

text: I understand how to calculate EIGN by hand but when I try to write code without function EIGN(), I did not have a clue. To calculate Eigenvalue is to count all the possible c in `det(A-c*I)=0`. A is the matrix, and c is a number.


title: PCA vs ridge regression for multicollinearity?

text: Orthogonalization via PCA and ridge regression are two common methods to account for multicollinearity for linear regression models. When would you use one over the other?


title: Transformation of standard normal distribution $\Phi(x)$

text: Can I write $\Phi(ax)$ as $b\Phi(x)$ for some constant $b$, please? FYI, $\Phi(\cdot)$ is the cdf for standard normal random variable. Is there such relationship?


title: Bayes Decision Boundary and classifier

text: Is it correct to say that the purpose of classifier (e.g. K-NN, Logistic Regression, LDA) is to approximate the Bayes Decision boundary?


I have the following data vector(in fact its `2392` data points long):               > (dput(head(data, 35)))     c(0.009917216, 0.158975732, -0.008065753, -0.018648987, 0, 0.283519647,      -0.108234184, 0, -0.151124958, -0.049062385, -0.12970766, 0,      -0.281741056, -0.04834985, -0.025923954, 0.167527783, 0.080834965,      0.457401398, -0.007328027, -0.026647639, 0, 0.146048389, 0.059651004,      -0.115781325, -0.022670935, -0.025330951, -0.011990763, 0, -0.051658383,      -0.054335988, -0.190313988, -0.031937749, -0.000498226, -0.064351546,      0.108053675)      Now I want to do a non parametric wilcoxon signed rank test for this data frame:               > (wilcox.test(data))              Wilcoxon signed rank test with continuity correction          data:  data     V = 1100747, p-value = 8.992e-12     alternative hypothesis: true location is not equal to 0      As you can see V is extremely high and I think its therefore probably wrong. Any recommendations what I am doing wrong? What I personally think is, that my R input could be wrong. Any recommendations how to specify it "much more" for a `Non parametric Wilcoxon Signed Rank test`. I really appreciate your answer!


title: Sparse principal component analysis

text: Does anyone know where I can find an algorithm (as well as an `R` implementation of it) to carry out sparse principal component analysis?


title: Why do Bayesian methods achieve sparsity via normal likelihood and Gamma hyperprior?

text: In the view of $\ell_1$ regularization, a natural way to incorporate sparsity into a Bayesian model is to add a Laplacian prior in my opinion. But instead of this, a usual fashion to handle sparsity in Bayesian model is the method proposed in the paper Sparse Bayesian Learning and the Relevance Vector Machine. What is the difference between these two in practice? Are there comparative studies of them? Thank you.


title: Soccer contest rule puzzled me

text: I have a little knowledge about soccer contest rule, there is a question I cannot understand. The dotplot below shows the numbers of goals scored by 20 teams playing in a city's high school soccer games on a particular day. Please use the following code to create the dotplot.               A <- c(0,1,1,1,1,2,2,2,2,2,2,3,3,3,4,4,5,5,7,10)       stripchart(A, method = "stack", offset = .5, at = .15, pch = 19,xaxt="n")       axis(1, at = 0:9, labels = 0:9)          1. One superstar scored six goals, but his team still lost, what are all possible final scores for that game? Explain.    2. Is it possible that all the teams scoring exactly two goals won their games? Explain.  Answers:   1. If the player scored six goals, his/her team must have scored either 7 or 10, but they lost, so they scored 7, and the only possible final score is that they lost by a score of 10 to 7.    2. No, there were six teams that scored exactly two-goals, but there were only five teams that scored less than two goals, so not all the two-goal teams could have won.  I read the answer over and over, but I still cannot understand.


title: How to test for differences in 3 related measures assessed pre- & post- treatment on the same sample?

text: I am designing a study to compare whether one specific treatment improves `Quality of life`. I am using a convenience sample ($N=20$). I hope to administer three different questionnaires to each subject pre- and post- intervention (all validated and standardised) to compare pre and post `QOL` measures. One of questionnaires will look at the severity and frequency of presenting symptoms pre and post treatment. I hope to analyse results to show whether symptoms are improved and what impact if any this has on each subjects Quality of life. Should I use paired t-tests?