title: Model comparison across multiple correlated responses

text: I have two multivariate linear regression models (multiple outcomes, i.e., the responses are a matrix), and I'm measuring their performance using $R^2$ in cross-validation, over these individual responses, and across all responses. I'd also like to quantify the difference in $R^2$ using a p-value, for example t-test (after Fisher transform of $R^2$ for better normality etc). Now, these responses are highly correlated with each other. Therefore, averaging the $R^2$ over the responses into a single estimate worries me as they are not independent, and a significance test that ignores this might be badly biased. Should I instead get a separate p-value for each outcome and then use some multiple-testing correction?


title: Experimental Design for Comparative Responses

text: Suppose I were looking to optimize the amount of certain spices in a chili spice recipe. The textbook experimental design would have me encode the amount of each spice in the design variable, choose some kind of design matrix, cook up a pot of chili for each row of the design matrix (one would have to take some shortcuts in the kitchen!), and assign a numeric dependent variable for each choice of the design variables. This latter operation is tricky for something like chili (or my real application, which is a secret); more likely a taste-tester could taste two different recipes and determine which is better, but assigning a numeric score is subject to all kinds of problems (drift over time, etc). So I have a two part question, and they are intertwined:   1. Given the results of a number of such comparisons, how should one pick the optimal vector of design variables? (As a first pass, I might look into "Noisy Sorting", but something less fancy, based on logistic regression, say, might be a better choice.)   2. Once the evaluation method has been chosen, how should I design the experiment?


title: Exclude highly correlated instruments

text: I have 3 financial instruments: A, B and C. And have 3 available correlation pairs:               (A, B) = 0.9     (B, C) = 0.9     (A, C) = 0.5      I need to exclude highly correlated instruments. i.e. only B, because A and C have low correlation (0.5). How I can do it algorithmically on every count of instruments?


title: Describing a best critical region of size $\alpha$ for testing $H_0 : \theta = 0$ against the alternative hypothesis $H_1 : \theta =1$

text: > Suppose a sample of size 10 is drawn from a distribution with probability > density function $f(x, \theta) = 2x^{\theta}(1-x)^{1-\theta}$ if $0<x<1$ and > $0$ otherwise, where $\theta \in \\{0,1\\}$. Describe a best critical region > of size $\alpha$ for testing $H_0 : \theta = 0$ against the alternative > hypothesis $H_1 : \theta =1$. We need $\frac{L(0)}{L(1)} \leq k$ for some $k < 1$. We find that $$\frac{L(0)}{L(1)} = \frac{2(1-x_1)2(1-x_2)...2(1-x_{10})}{2x_12x_2...2x_{10}} = \frac{(1-x_1)(1-x_2)...(1-x_{10})}{x_1x_2...x_{10}}.$$ Now I want to simplify $\frac{(1-x_1)(1-x_2)...(1-x_{10})}{x_1x_2...x_{10}} \leq k$ to see if I can turn this into a distribution that looks more familiar. But for some reason I wasn't able to do anything, so I was wondering if somebody would give me a hint so I can continue... Thanks in advance


title: Predicting one daily variable from another in SPSS

text: Note: There are similar questions to this one, but they don't seem to get at quite what I'm trying to figure out. I have a week's worth of daily data with a number of variables, including nighttime sleep and daytime mood. I want to predict mood during the day from sleep the night before. The questions are asked in the morning, so S1 is last night and M1 is yesterday, meaning I want to see if M2 is predicted by S1, M3 by S2, and so on. I'm having a hard time figuring out the best way to do this in SPSS. Any ideas?


title: Random Forest Regression Overfitting - Quantile Test on Test Data

text: I have fit a random forest regression model to training data (used 65% of data for training). The data has approximately 40,000 observations and 100 features. I fit a random forest regression in R with the following parameterization:               randomForest(formula = Response ~ .,                  data = crs$dataset[, c(crs$input, crs$target)],                  ntree = 500, mtry = 32, importance = TRUE, replace = FALSE, na.action = na.roughfix)      My understanding is that for Random Forest Regression problems, it is best to use approximately 1/3 of the candidate variables for each tree (rather than square root for classification problems) so that is why I have tried 32 variables per tree. After applying the model to my test holdout data set (approximately 35% of data) the model appears to be overfit which I am confused by as I thought Random Forests were supposed to be rather resistant to overfit (which has been my experience in prior usage of them). Here is a comparison of the average predicted vs. average actual value on Test data sorted ascending by predicted value (predictions grouped into deciles).                Prediction_Decile  Avg_Prediction  Avg_Actual  Ratio:Actual/Predicted      1  4,570   6,343   1.388       2  5,939   7,085   1.193       3  6,789   7,429   1.094       4  7,576   7,982   1.054       5  8,320   8,981   1.079       6  9,105   8,796   0.966       7  9,954   8,657   0.870         8  10,977  9,306   0.848       9  12,304  9,814   0.798        10 14,653 10,195 0.696 As you can see the ratio of Actual to predicted value is steadily decreasing as the predictions increase which is why I think I am overfitting. Any tips or advice on what may be causing this or how to tune the model to avoid this problem? The model appears to be doing a decent job of ordering the test observations,but a much poorer job of fitting them.


title: comparison of sets of binary sequences

text: I want to compare three sets of binary sequences and assess whether they differ statistically. Each sequence is 100 bits long (the length 0-100 represents the axis of an elongated tissue : if there is a 1 at position i, it means that a cell from that particular lineage was found at this position). The sets contain 26, 20 and 29 sequences respectively. Each set comes from a different pool of cells and I want to know whether the 3 pools differ significantly as far as their capacity to colonize along the tissue is concerned.


title: Logistic Regression Sample Size and Outcome Period

text: I am to run logistic regression on approximately 900,000 observations (each a unique customer) in order to determine propensity to purchase a given product or propensity to hold a given product. I would prefer building propensity to buy models but am not sure if I have enough data/how to best construct the outcomes. I must build models for 7 different products, and those models are to be treated independently (i.e. the outcome period for each model can be different). The following table shows the products that I'm going to model and the respective count of unique customers that either purchased the products from December 13-Jan 14, Nov 13-Jan 14, Oct 13-Jan 14 or simply hold the product irrespective of date. ![Product holding table](http://i.stack.imgur.com/h2lUu.png) I have almost no data for most products in the months Dec-Jan and Nov-Jan, but I start seeing more if I look at the months Oct-Jan. **Here - my questions, given the total of 900,000 observations:** 1) Is there any sense at all in building propensity to buy models for those products where I have data on product purchase between Oct-Jan (i.e. setting the dependent variable of those customers who bought a given product between Oct-Jan to 1 ("event")), or do I simply have too little data and must therefore rely on propensity to hold models for all products? 2) How do I determine the "non-event" (e.g. 0) segment of my dependent variable? Should I randomly sample from the population of customers who do not hold the product (or who did not buy the product if I can use propensity to buy models)? How many observations should I draw (i.e. do I want a 50-50 breakdown of events and non-events, or should I use more non-events, as I have more information on them? How many more non-events should I use?)?


title: When is the standard error of the mean impossibly large for a given data range, when we know the sample size?

text: While reading an article (published and peer-reviewed, referenced below) investigating proportions of different tissue types in the human brain I came upon a table presenting data on the different subgroups in the study. Data was presented as means +/- SE (which I assume to be standard error, although I could not see it specified in the text) and range. One example is the age of one of the subgroups: n= 23 mean age = 66 SE = 2.9 range of age 60-69 This puzzled me. If the SE is 2.9, then given the sample size of 23 the SD of the age should be: SD = sqrt(23) * 2.9 = 13.9 This would lead to a standard deviation that is quite a lot larger than the specified range for the age, which should not be possible. I'm a beginner in the world of statistics, so my question is simply if it from this is possible to conclude that one or several of the values given must be incorrect, or if I'm missing something that would make this data to make sense. Thank you! Reference: Guttmann CR, Jolesz FA, Kikinis R, et al. White matter changes with normal aging. Neurology 1998;50:972-978.


title: Poisson process thinning females and males arriving

text: Rock tickets are sold at a ticket counter. Females and males arrive at times of independent Poisson processes with rates 30 and 20. > What is the probability that the first three customers are female **My Work** Let $F(t), M(t)$ be the number of females, males respectively, that arrive up to time $t$. Then if I condition on $M(t) + F(t) = 3$, I get $$\Bbb P(F(t) = 3, M(t) = 0 | M(t) + F(t) = 3) = {\exp(-30t) (30t)^3/3! \cdot \exp(-20t)(20t)^0/0! \over \exp(-50t)(50t)^3 /3!}$$ Is this accurate? My answer, after canceling, works out to ${30^3 \over 50^3}$ which seems a bit too simple? It's exactly $\Bbb P (\text{female arrival})^3$.


title: How to set initial values for Holt-Winters backcasting?

text: I'm trying to figure out how to set initial seasonal values for backcasting of the multiplicative Holt Winters. As pointed out in this thread: Holt Winters initialization using backcasting like SPSS, you can set: $$ {\rm Seasonal\ Index}_i = \frac{m_i + b_i}{\sum(m_i + b_i)}. $$ This method produces the same results as SAS and SPSS when $m_i + b_i > 0$. But what should be done when $m_i + b_i < 0$?


title: Factoring in self selection bias

text: I am hoping to use some self selection survey data (it's from those awful things that pop up on the start page of a website asking if you have 5 minutes to spare). This is for business decision making so obviously I don't want to overstate the importance of any findings. I have found one public access paper that deals with "Heckman Corrections" here:- http://onlinelibrary.wiley.com/doi/10.1111/1467-6419.00104/pdf Am I right in thinking that there are not any statistical methods in common use for correcting for self selection bias? Are there any simple ones that would be useful for a quick and dirty decision making process, even if they're not perfect?


title: Notation of a resampled mean

text: For a work that I have to document, at some point I do N resampling of the initial sample and generate N means out of them (some kind of bootstrapping). Pretty easy, however I don't know how to write it mathematically, I was thinking of something like: $\tilde x = \sum_{k\in\Gamma} x_k$ and then describe $\Gamma$ as a set. Any idea?


title: An MGF exercise

text: I would appreciate some help in the following exercise: Let the Random Variable $Y$ have the moment generating function $$M(t)=\exp\\{q(\theta)-q(\theta+t) \\}$$ where $q(.)$ is a function. If $E[Y]=\theta$, prove that that $Y\sim N(\theta,1)$. **Hint: Consider $M\prime (0) =\theta$ and solve the differential equation.** * * * Starting from $M\prime (0)$ then, we see that $$M\prime (0)=-q \prime (\theta)=\theta $$ A solution of this differential equation is $ q (\theta)=C\exp\\{-\theta \\}$, with $C$ being a constant which I am having trouble definitizing. I fear that I have reached a dead end. I think I need to show that the MGF takes the form $\exp\\{\theta t +(1/2) t^2 \\}$, i.e. the MGF of a normal distribution with mean $\theta$ and variance $1$. In order to do that, however, I have to definitize the constant which due to the lack of boundary conditions, I don't seem to be able to. Any advice on how to proceed from here? Thank you.


title: R ggplot2/ezPlot: Plotting 3x3 RM ANOVA design per quantile of covariate with error bars

text: First of all, my data: http://www.pastebin.ca/2599202 (I hope this is not too inconvenient, because I fail creating good fitting example data) What I basically need, is a plot like: ![enter image description here](http://i.stack.imgur.com/6GO8L.png) I plotted the repeated measures factor `time`(x-axis, 3 levels) against `ias` (continuous dependent variable) for my 3 experimental groups. I did this 4 times (for each quantile of my trait-measure MIHT, `miht.binned`, .25 - 1.00). I have to admit I am not really an R professional and the `ggplot2` manual simply is an overkill for me. I created the plot with `ezPlot` (from `ezANOVA`) and only managed to do a bit layout tweaking with `ggplot2`:                  PlotIAS = ezPlot(                                          data = MyData           , dv = .(ias)           , wid = .(id)           , between = .(GROUP, miht.binned)           , within = .(time)           , x = .(time)           , split = .(GROUP)           , col   = .(miht.binned)           , x_lab = 'time of measurement'           , y_lab = 'IAS Score (Mean)'           #, do_bars = FALSE           , type = 3         )                   PlotIAS = PlotIAS +             theme(             panel.grid.major.y = element_line(colour = "gray80", size = NULL, linetype = NULL,                                                 lineend = NULL)             ,panel.grid.minor.y = element_line(colour = "gray90", size = NULL, linetype = NULL,                                                lineend = NULL)             ,panel.grid.major.x = element_blank()                        ,panel.grid.minor.x = element_blank()             ,legend.background = element_rect(fill = NULL, colour = "black")                   ,panel.background = element_rect(fill = "white", colour = "white", size = NULL,                                               linetype = NULL)           )          print(PlotIAS)      I did'nt find any information about these error bars ezPlot creates. They seem to be all the same. I just need to have error bars with SE or CI. I don't know if it is possible to add these in my ezPlot-based code (and how?) or if one has to create a complete new ggplot object for that (which is an overcharge for me...). Help is highly appreciated.


title: Feature selection methods comparison

text: I recently run a project that involves a feature selection step before further pattern recognition. The number of features for our data set is very large and instead of running greedy forward/backward feature selection. I simply run the Student T-test on the features and ranking them according to the p-values to find the most 'powerful' features for my classification task. However, I got some argument from my colleagues that this kind of Student T-test based feature selection method is not optimal. For example, the forward or backward feature selection can find the best features until the classification performance can't be improved (therefore, they can find the optimal feature subset). I am wondering if anyone knows some references (papers, thesis) that run some experiments to prove that Student T-test based feature selection method can perform similarly like the exhaustive methods? Or is there any mathematical background to underpin this? Thanks a lot. A. I found a relevant paper on PlosOne: A paper


title: Analysing significance in grouped variables

text: I have imputed a dataset which has 20 observations. Participants were asked to rate the importance of 13 variables under a certain circumstance, from 0 to 25. I imputed this using `mice(imp=25, max_it=25)`, and got a nice dataset after a few hours. I'm now analysing it using boxplots per circumstance (thus 13 boxplots at a time), but wonder how I can measure the significance per circumstance of boxplots. ![Groups of boxplots](http://i.stack.imgur.com/2Kojw.png) In the above example it is clear that the circumstance `Revenue: 0 to 2 Million` is different from the `Revenue: 10 to 50 Million` and possibly from `2 to 10 Million` but how can I measure this statistically? My data has a normal distribution. A friend recommended ANOVA, but that's univariate and I feel that 20 observations is too little... I saw: How to compare two groups on a measure of social skills that includes 5 subscales where each subscale is number correct out of 12? but not really sure if that applies since its non-parametric.


title: Calculating boxplots with imputed data

text: I have imputed a dataset which has over 200 variables, and 20 observations. In worst case, 80% of the data is missing, in the best case 100% is available. 5 out of the 20 participants provided data for all the variables. Participants were aksed to rate the importance of 13 variables under a certain circumstance, from 0 to 25. I imputed this using `mice(imp=25, max_it=25)`, and got a nice dataset after a few hours. I'm now analysing it using boxplots per circumstance (thus 13 boxplots at a time), but wonder how I should do this. From a Dutch textbook, I understood that you do the tests for each imputation, and then average the results. Does that mean that I should calculate the data for each variable/imputation, and then divide by 25? It seems a bit crude... I also wonder how much it differs for these descriptive statistics...


title: Model based on many variables

text: Let's say I have a matrix of values for many different variables Y1..Y1000 at X=1,2,3..,10. Some of these variables are directly correlated with X, some follow different shapes (e.g. a normal distribution) and some are just random. I want to build a model to predict X based on given values of Y1..Y1000. What would be the correct approach for this? I assume a simple linear regression would not be feasible because of the number of variables and the fact that not all variables are linearly dependent on X.


title: Metropolis Random Walk Algorithm - Covariance Matrix

text: I need to calculate an integral of a 3 variable density $f(x_1, x_2, x_3)$, defined in a domain $D$, in the region $$T=\\{x \in R³|f(x) \geq M \\}$$ using MCMC , where $M$ is a given positive constant. Defining the function $g(x)=1\;if\;f(x)\geq M$ and $g(x)=0\;if\;f(x)<M$,we have $$\int_T f(x)dx = \int_Dg(x)f(x)dx$$ So the MC estimative to the integral will be $$ \bar{I}= \frac{1}{N} \sum_{i=1}^N g(x_i),\; x_i \sim f$$ My value $M$ is $f(a,b,c)$, where $(a,b,c)$ is the maximum of the function $f$ with the constraint $$a = \rho \; c \; \Gamma\left(1 + \frac{1}{b}\right)$$ for $0<\rho<1$. My problem is when I try to use the Metropolis Random-Walk Algorithm with a trivariate normal (the mean is the last state) to find the sample distributed as $f$. If I change $\rho$, $(a,b,c)$ becomes closer or not from the point of maximum of the function (which I am using as initial state of my chain). If the point is close to my distribution maximum, I don't want the chain to move so distant when I run the chain because, in this case, $g$ will be zero almost all the time. But if I fix little variances in my covariance matrix, when the point is far from the distribution maximum (so I want the chain to walk a bit more in the distribution), it will give me wrong estimatives for the integral. So I think my covariance matrix will depend on $\rho$. But how can I find the right covariance matrix coefficients? I couldn't analyze the constraint so I could find the right relation. I hope my doubt is clear. Thanks in advance!


title: Principal Component Analysis: whether a variable is significantly loaded on a principal component or not?

text: Often, a variable is considered to be significantly loaded on a PC if its loading value in the loading table is above a cut off value (suppose 0.4 or 0.5 in some published cases). Is there any statistical/mathematical method to check whether a variable is significantly loaded on a particular PC or not?


title: Showing that a statistic is ancillary for a parameter

text: Working through a HW problem, and a hint is that for a decision rule $$T(X) = \frac{X_{(1)} + X_{(n)}}{2}$$ Then $$T - \bar{X} $$ is ancillary. Intuitively this makes complete sense, but I am failing to see how to show this. I thought about going to the pdf of an ordered statistic: $$g_i(x) = \frac{n!}{(i-1)!(n-i)!}[F(x)]^{i-1}[1-F(x)]^{n-i}f(x)$$, which for i = {1,n} reduces down to $$g_1(x) = n[1-F(x)]^{n-i}f(x), \:\:\:\: g_n(x) = n[F(x)]^{n-i}f(x)$$ But I would probably need the convolution from there to get the pdf of T. I feel like I'm way off and there is something simpler here to show that subtracting off $\bar{X}$ makes T independent of theta. Oh by the way $X \sim N(\theta, \sigma^2)$, $\sigma$ known Thanks! * * * In response to comments: Then $$T - \theta = \frac{1}{2} \left[ (X_{(1)} - \theta) + (X_{(n)} - \theta) \right]$$ and since $X_{(1)}$ is a location family $$X_{(1)} = Z_{(1)} + \theta$$ where $Z_{(1)}$ is the first ordered statistic of a standard normal vector, which implies $$X_{(1)} - \theta = Z_{(1)} \sim \left[ 1 - \int\limits_{-\infty}^z (2\pi)^2\exp(\frac{-z^2}{2})dz \right]^{n-1} (2\pi)^{-1/2}\exp(\frac{-z^2}{2})$$ which does not depend on theta. Similaryly $X_{(n)} - \theta$ also does not depend on $\theta \Longrightarrow T - \theta$ is ancillary for $\theta$. One big problem though, I am subtracting $\theta$, not $\bar{X}$. Unless were talking about a first order ancillary statistic (where $E\bar{X} = \theta$), I am once again confused, lol. Again, thanks for any additional help!


title: Unassigned classes in Multiclass SVM with One vs All Approach

text: How do you handle unassigned classes in multiclass support vector machines (multiclass SVM) with the One vs All approach? Lets say my training data has three classes A, B, and C. I use 3 SVM one-vs- all classifiers. For a particular test instance, all classifiers say 'Not A', 'Not B' and 'Not C' respectively. How do I assign this instance to one of the classes? I tried using an 'unassigned class' class, but the error rates are too high. Another question:   When two or more classifiers assign a test instance as positive, which one should I go with. I read somewhere that you should use "output function of the SVM". Now what's that? Matlab gives the following svmStruct:   1. SupportVectors   2. Alpha   3. Bias   4. KernelFunction    5. SupportVectorIndices   6. Others How do I use these to generate the "output value"? Also, I did not quite understand Platt's technique. How do I implement it using elements of svmStruct. This question is related but went unanswered.


title: Which test should I run?

text: I have RNA sequencing data (= expression of a certain protein) from the lung tumor of 14 patients, 5 non-smokers and 9 smokers. If I want to compare the expression level of that protein of smokers vs non- smokers, which test should I use? The dataset doesn't seem parametric, so I usually try a Wilcoxon test then, but these are not paired data, no? Thank you very much!


title: Using R and hidden Markov models to predict pseudorandom binary sequences better than P > .5?

text: I would to predict a pseudo-random binary sequence of 0's and 1's. I am thinking of using the HMM package in R. I have a binary sequence like ... 0 1 0 0 1 1 0 1 x(n+1)? with thousands of values. What I would like is to get automatic estimation for x(n+1), computing their prob.matrix automatically every time. The question is to know, if a HMM (or any other algorithm) can provide P > .5 for binary events that pass random tests, BUT are generated by deterministic algorithms, like LCG. Something like that; http://www.stanford.edu/class/stats366/hmmR2.html BUT using a data file with the binary sequence, so prob.matrix can't be estimated before. thank you


title: How to calculate probability of getting 1 heads in 2 coin flips or 1 heads in 12 flips, and so on

text: If I wondered about the probability of getting:   * Only one heads in two tosses - 2/4   * Only one head in three tosses = 3/8 or 37.5% But I just counted on my fingers, how do you do it for big numbers?


title: How to analyse a moderator effect in hierarchical multiple regression with covariates?

text: I am trying to find an appropriate method to examine the significant interaction effect in my hierarchial multiple regression analyses. I am looking at moderation analysis but most of of the methods (macros etc) do not take into account the covariates in my regression equations. I was advised to try out the Stata programme but at the moment, I have no access to it.   * Is there another statistical program that I can perform moderator regression with covariates?   * How can I perform simple slopes analysis with covarites?


title: Plotting sparklines in R

text: I would like to use R to plot out something like this: ![Something like this](http://i.stack.imgur.com/oY7UO.png) It would seem possible but highly complex to keep track of the coordinates, width, height, etc. Intuitively it would seem best to treat each cell as a new plot and transform the coordinates for each cell. Is there a way to do this in R? thanks!


title: How to get started with applying item response theory and what software to use?

text: ### Context I have been reading about item response theory, and I find it fascinating. I believe I understand the basics, but I am left wondering how to apply statistical techniques related to the area. Below are two articles that are similar to the area I would like to apply ITR in:   * http://www.jstor.org/stable/4640738?seq=7   * http://www.ncbi.nlm.nih.gov/pubmed/21744971 The second being the one I would actually like to extend at this point in time. I have downloaded a free program called jMetrik, and it seems to be working great. I think it may be too basic as far as IRT goes, but I am unsure. I know the "best" way would likely involve learning R; however, I don't know if I can spare the time to tackle that learning curve. Note that we have some funding to purchase software, but from what I see, there doesn't seem to be any great IRT programs out there. ### Questions   * What are your thoughts on the effectiveness of jMetrik?   * How would you suggest I go forward in applying IRT?   * What are the best programs for applying IRT?   * Do any of you use IRT regularly? If so, how?


title: Estimate the nearest of N random points in a box in E^d?

text: I have N uniform-random points $p_j$ in a box in $E^d$, $a_i \le x_i \le b_i$, and want to estimate the expected distance of the point nearest the origin in $L_q$: $\quad$ nearest( points $p_j$, box $a_i .. b_i$, $q$ ) $\;\equiv\;$ $\min_{j=1..N}$ $\sum_{i=1..d}$ |$p_{ji}|^q$ The box might straddle 0 in some dimensions, $a_i < 0 < b_i$.   Also I'd like the estimator to work for $0 < q < 1$ too, the "fractional metric". (16 March) Let's try to do a simpler case: $L_1$, unit cube 0 $\le x_i \le$ 1.   Geometrically, we want (correct me) the diagonal slice of the unit cube with volume $\frac{1}{N+1}$.   (Does anyone have a picture or applet of a 3d cube sliced into equal-volume slices ?) If d is big enough for a central limit theorem to hold,   $\quad \sum_{i=1..d} uniform_i \ \sim\ \mathcal{N}( \frac{d}{2}, \frac{d}{12} )$;   so $\mathcal{N}^{-1}( \frac{1}{N+1} )$ gives approximately the cut, and the expected nearest distance, that I want. In Python with scipy.stats, this is               def cutcube( dim, vol ):         """ cut a slice of the unit cube in E^dim with volume vol             normal approximation: cutcube( 3, 1/6 ) = .339 ~ 1/3             vol 1/(N+1) -> E nearest of N random points in the cube ?         """         return norm.ppf( vol, loc=dim/2, scale=np.sqrt( dim/12 )) / dim          cutcube( dim=2, vol=1/10 ) = 0.24     cutcube( dim=4, vol=1/10 ) = 0.32     cutcube( dim=8, vol=1/10 ) = 0.37     cutcube( dim=16, vol=1/10 ) = 0.41     cutcube( dim=32, vol=1/10 ) = 0.43


title: Putting a confidence interval on the mean of a very rare event

text: I'm simulating an extremely rare event (the detection of weighted photon packets in a highly absorbing material). For example, I may simulate the transmission of 1e9 photons and only detect 10 of them (and their weights vary from $0^+$ to 1). I'm calculating the mean value of this occurring (what is the mean of $\sum w_{received}/N_{transmitted}$, where $w_{received}$ is the weight of the received photons and $N_{transmitted}$ is the number of transmitted photons) and put a confidence interval on this. My understanding is that the population mean, $\mu$ lies within one standard Normal deviation of $\bar{x} \pm S_N/\sqrt{N}$, where $\bar{x}$ is the sample mean, $S_N$ is the sample variance, and $N$ is the number of samples. My question is this, how do I know when I've received enough "positive" samples, or received photons, to get an accurate estimate? Since I'm dealing with such a small probability, the total number of samples, $N$ is extremely high and $S_N/\sqrt{N}$ will always be small. I've seen notes saying the CLT only holds for samples sizes > 30, but would that apply to something like this? i.e. $N_{received}$ > 30? _Summary_ : I'm drawing samples from an unknown distribution and computing their mean. Sample values range from 0 to 1, with 0 being far more common. I'm computing the mean (which is very close to 0) and I want to put a confidence interval on this. **Edit** From my description above, there's two ways to think about the problem. The first is a binomial distribution of yes/no results (yes, received; no, not received). However, the more complex, and in practice, more useful, distribution is somewhat of a weighted binomial (not sure if that's the right term), where each received photon "packet" has a power associated with it, that ranges from 0 to 1. That is the case I care most about putting a confidence interval on. i.e. my sample values will predominately be 0's with a small number of samples with values 0 < x <= 1 **Edit 2** The PDF will look something like the crude drawing below. Note that it is not to scale - there will be _much_ more 0-weights than non-zero weights. This is an _example_ pdf of samples. I don't know what this will look like before I do an experiment - I just know my samples will range from 0 to 1, with a predominate amount being 0. ![Example PDF of sample values](http://i.stack.imgur.com/Vqs7F.png) **Edit 3** The term "weight" seems to be a bit misleading. "Sample value" or "power" might be a better term. My simulation tracks photon movement through water. Most photons never reach the receiver and therefore have a power of 0. As photons move through the water their power is reduced (to model a group of photons moving with a portion being absorbed as they move along). As these photons are received they have weights that range from a max of 1 (no loss from the travel) down towards a number approaching 0.


title: Using QDA for Non-Gaussian distributions

text: I am evaluating a Quadratic Discriminant Analysis (QDA) classifier on a high- dimensionality feature set. The features come from highly non-Gaussian distributions. However, when I transform the features to each have a Gaussian distribution in two different ways, the resulting classifiers performs worse than QDA applied on the raw features with the following three metrics:   * Accuracy (classes are balanced)   * Area under the ROC   * A probabilistic metric The first way of transforming each feature to a Gaussian distribution disregarded the class. For each feature, it found the parameters of the the corresponding distribution, used the distribution's CDF function or approximation for each data point, and then used the Gaussian's inverse CDF. The second way did the same thing, but included class labels and treated data from each class independently. Any idea why this occurs? I have confirmed that it is not due to ...   1. A bug in the code   2. Distributions changing from train set to test set


title: Generative model that penalizes clumping of data

text: I'm interested in modeling a generative process that encourages data to be "evenly distributed" over its support, i.e. clumping of data points is penalized. For example, if I have a mixture distribution and I set all components to have equal weight, it will tend to generate points equally between all the components, but this property isn't enforced. The problem here is that if we evaluate a data set under this likelihood, it gives the same probability to a data set that is evenly distributed as it would to a data set that is clumped under a single component. To illustrate a "better" model that does penalize clumping, consider the following generative process. Beginning with a mixture distribution with $N$ components, sample a point in the standard way: first sample an index $i$ from the component weights, and then sample a point from that component. Now decrease the weight of the $i$-th component to be a certain percentage of its original weight (call this percentage $p$), and distribute that mass to the other $N-1$ components. Repeat $M$ times. This process generates points from the mixture distribution in a way that encourages points to be evenly distributed between the components. The parameter $p$ adjusts how strongly this even-distribution of points is enforced--when $p=1$, this is a typical mixture distribution. When $p=0$ and $M = N$, this enforces a strict one-to-one correspondence between data points and mixture components. The problem with this process is that although it's straightforward to sample from, evaluating the data likelihood under this model appears to be intractable, due to the need to consider all possible orders in which the data could have been selected. Are there any well-known models that penalize clumping in this way?


title: How well does R scale to text classification tasks?

text: I am trying to get upto speed with R. I eventually want to use R libraries for doing text classification. I was just wondering what people's experiences are with regard to R's scalability when it comes to doing text classification. I am likely to run into high dimensional data (~300k dimensions). I am looking at using SVM and Random Forest in particular as classification algorithms. Would R libraries scale to my problem size? Thanks. EDIT 1: Just to clarify, my data set is likely to have 1000-3000 rows (perhaps a bit more) and 10 classes. EDIT 2: Since I am very new to R, I will request the posters to be more specific where possible. For example, if you are suggesting a workflow/pipeline, please be sure to mention the R libraries involved in each step if possible. Some additional pointers (to examples, sample code etc.) would be icing on the cake. EDIT 3: First off, thanks everyone for your comments. And secondly, I apologize, perhaps I should have given more context for the problem. I am new to R but not so much to text classification. I have already done pre- processing (stemming, stopword removal, tf-idf conversion etc.) on my some part of my data using tm package, just to get a feel for things. tm was so slow even on about 200docs that I got concerned about scalability. Then I started playing with FSelector and even that was really slow. And that's the point at which I made my OP. EDIT 4: It just occurred to me that I have 10 classes and about ~300 training documents per class, and I am in fact building the termXdoc matrix out of the entire training set resulting in very high dimensionality. But how about reducing every 1-out-of-k classification problem to a series of binary classification problems? That would drastically reduce the number of training documents (and hence dimensionality) at each of the k-1 steps considerably, wouldn't it? So is this approach a good one? How does it compare in terms of accuracy to the usual multi-class implementation?


title: p-value with multimodal PDF of a test statistic

text: I have opened a thread about p-value under the title "Understanding p-value" and gotten two answers and some comments. I think my questions in the thread is somewhat diverse and want to clarify my question more explicitly based on the discussion in the thread. Two different definitions of the p-value were suggested in the thread. **definition 1** > The p-value is $\int_{\\{x\,:\,f(x) \le f(x_o)\\}} f$. **definition 2** > The p-value is $\int_{\\{x\,:\,x_o \le x\\}} f$. In both of the definition, $f$ is the PDF of a chosen test statistic under the null hypothesis and $x_o$ is the observed value of the test statistic. I think the two definitions are clear and complete enough. (The p-value concerns data, a null hypothesis and a chosen statistic only. It does not concern the alternative hypothesis or other things.) The role of the p-value is to quantify how likely the observation is under the null hypothesis. Small p-value means the observed data is weird (ie. unlikely) under the null hypothesis and the assumed null hypothesis should be rejected. The definition 1 measures this weirdness in terms of $f(x_o)$, the probability density of the observed test statistic. So the definition integrates $f$ over the values of the test statistic that have smaller probability density (ie. more weird) than the observed one. The definition 2 measures the weirdness in terms of the distance of $x_o$ from the most likely value of the test statistic, if the most likely value is well defined. So the definition integrates $f$ over the values from the observed one to tail (ie. more weird region). If $f$ is unimodal, both of the two definitions seem reasonable. If $f$ is multimodal, however, I think the definition 2 is not reasonable. For an example, let's assume that $f$ is bimodal and $x_o$ is somewhere in the low probability density region between the two peaks. Then the most likely value is not well defined and the distance of $x_o$ from the most likely value cannot be reasonable measure of the weirdness. The p-value calculated along the definition 2 may be very large, whereas the observation $x_o$ is obviously weird because of its low probability density. The definition 1 still works in this case as it gives small p-value. I am not a statistician and I don't know which one of the definitions is "the right one" that statisticians usually use. Most of the materials I have seen before explain p-value in the sense of the definition 2. But, I encountered the definition 1 in Zag's answer of the old thread for the first time and was persuaded. What is the exact definition of the p-value? If it is not the definition 1, I'd like to know rationale for the right one and shortcomings of the definition 1.


title: Is there a site to post my survey to so I can get a sample representative of the population?

text: This is only for my high school senior project, so it doesn't have to be perfect. I'm doing a project on global warming, and I want to survey people for their opinions. I know that if I use a convenience sample of my classmates, I'll have a lot of bias. I was wondering if there was a site on the Internet where I could post my survey for random people to answer, so that I could get as close as possible to an SRS. If not, I'd be interested in some other suggestions for ways to carry out my survey.


title: Find the exact distribution of the MLE estimator and $n(\theta-\bar{\theta})$ exact and limiting distribution

text: I have another problem and the density is something like this: $$f(x,\theta)=\frac{1}{\theta} \text{if }0\le x\le\theta \text{ for some }\theta>0$$ $$f(x,\theta)=0 \text{ otherwise}$$ Given the sample $X_1,...,X_n$ be iid. I found the MLE to be $X_n$ which is the the maximum of the sample. I'm supposed to find the exact distribution but I'm not sure what do they want? Do I derive the density and CDF of a sample maximum? Lastly I'm supposed to find the exact distribution and asymptotic distribution of $$n(\theta-\hat{\theta})$$ I know about Lindeberg Levy-Central Limit Theorem but I'm not sure how to apply that into this case or deriving the exact distribution.


title: Interpreting multinomial logistic regression output in R

text: I am trying to perform multinomial logistic regression on my data which is as below(just the header). ![enter image description here](http://i.stack.imgur.com/LLZjK.jpg) "category" is my target variable and all other variables are independent variables. category has values like 1,2,3,4,5,6,7,8,9. My main motive is to predict category from independent variables. Here is summary of my data set. ![summary](http://i.stack.imgur.com/hXuNg.jpg) I used following command to perform multinomial logistic regression:               > mod=multinom(category~hlen+iplen+ipttl+iptype+tcpsport+tcpdport+tcpsec+tcpack+tcpwindow+tcpchksum+date_time, data=train)      ![command](http://i.stack.imgur.com/bZMWw.jpg) I got the following output but i don't know how to interpret it? ![enter image description here](http://i.stack.imgur.com/c2QQR.jpg) What should be the starting point or is there other way to do so? (For example I know how to interpret linear and logistic regression output to create regression equations) Please right click any image and select view to see it clearly. Thank you,


title: Complete interpretation of Cox regression output in R

text: I need help with complete interpretation of this output(below). Can someone explain all (almost all) estimators in this output: what they mean and what is the substantive interpretation. In addition: Which tools would you recommend to test for these estimators? Thank you! * * *               inter.coxph <- coxph( Surv(time,censor)~age+drug+age*drug, method="breslow")     summary(inter.coxph)           n= 100                   coef exp(coef) se(coef)      z       p       age  0.0942     1.099   0.0229  4.110 0.00004     drug  1.1859     3.274   1.2565  0.944 0.35000     age:drug -0.0067     0.993   0.0337 -0.199 0.84000               exp(coef) exp(-coef) lower .95 upper .95       age     1.099      0.910     1.051      1.15     drug     3.274      0.305     0.279     38.42     age:drug     0.993      1.007     0.930      1.06                Rsquare= 0.295   (max possible= 0.997 )      Likelihood ratio test= 35    on 3 df,   p=1.21e-007      Wald test            = 32.2  on 3 df,   p=4.83e-007      Score (logrank) test = 35.2  on 3 df,   p=1.13e-007


title: SPSS: Comparing Regression Coefficient for 2 Models

text: Hope you guys could help me with a question I've been stuck on for a while. I'm currently writing my thesis on how MRT (the railway system in Singapore) accessibility affects prices of public housings (HDBs). So far I've gotten quite far but have reached a wall. I have collected 11,000 HDB transactions, all which are within walking distance of MRT station (<750m to the nearest MRT Station or <10min walk to the nearest MRT Station). Apart from the usual structural variables:   * _Age_   * _Floor_   * _Size_ I also managed to obtain locational variables such as:   * Time taken to walk to nearest railway station ( _time_walk_ )   * Time taken to commute, via train, from the nearest station to the CBD station ( _time_train_ )   * Total Traveling Time, _time_walk + time_train_ ( _TTT_ ) With these variables, I ran a multiple regression with _Price_ as the DV and _Age, Floor, Size, time_walk_ , and _time_train_ as the IVs. This then produced regression coefficients for the DV. With these coefficient it allowed me to analyse the quantitative impact that, with everything else constant, (i) each additional walking minute, and (ii) each additional commuting minute on the MRT, had on the pricing of HDB. However, what I would like to investigate is, do residents living at different distances from the CBD value _time_walk_ differently? I understand that I can't create 3 models (shown below), each containing only the relevant details (eg. 0-9 mins train time, 10-19 mins train time ...etc) as the n number would be different, thus, comparing the coefficient estimates wouldn't be fair Model A: (0-9 mins time_train): How would Walking_Time affects house pricing? Model B: (10-19 mins time_train): How would Walking_Time affects house pricing? Model C: (20-29 mins time_train): How would Walking_Time affects house pricing? Any advice would be greatly appreciated. Regards, Slim


title: How to encode the time in an Markov Jump Process

text: Assume I have a time continuous Markov process with only a few (let's say 2, in reality a few more, like 5) states. I can measure these states in real data and I want to estimate the transfer probability matrices of the process. The states occur in random order and in random length. So for example state 1 hold on for 60 secs, then the system switches to state 2 and stays there for only 5 secs before returning to state 1 for 10 secs and so on. My question is, how do I incorporate the time when calculating the TPMs? My idea is to introduce new substates, e.g. state1_short, state2_mid and so on. The problem will be to determine what short, long etc. shall be. Is there another way, what are the keywords I could search for? Thank you very much!


title: Integration via importance sampling - expectation and variance

text: I need to calculate $$\int_A f(x) \; dx$$ via Monte Carlo Importance Sampling and, for simulations purposes, I need to calculate the variance and mean of my estimator. In importance sampling method, we have a random variable $X$ with pdf $g$ defined only in $A$. Then, the Monte Carlo estimative of the integral is $$\frac{1}{n}\sum_{i=1}^n \frac{f(x_i)}{g(x_i)}, x_i \sim X$$ which converges to the given integral, since $$\int_A f(x) \; dx = \int_A \frac{f(x)}{g(x)}g(x) \; dx$$ I need to calculate the expectation and variance (or standard deviation) of the Monte Carlo estimative for my integral. Of course that, since it is an simulation, I need to calculate it using sums and in function of my $n$. How can I do this? I thought that $$ \sigma = \left(\frac{1}{n} \sum_{i=0}^n \left(\frac{f(x_i)}{g(x_i)}-I\right)^2 \right)^{\frac{1}{2}}$$ but I think that is not correct. For real, I need to calculate the relative error of my simulation, which I think is given by the ration of the standard deviation and mean, but with this formula, it does not decreases as expected. Thanks in advance.


title: Is statistical significance calculated differently for a sub-set of responses than it would be for the total response set?

text: Please pardon my lack of terminology - not a stat person. Here's the background & 2 questions: If I have a survey sample set of 100 (within a total population of >10 million), then I believe that would give me a 9.8 confidence interval at 95% confidence level. That I understand. However, what I don't understand is how to figure out what the confidence interval/level (or in my layman's terms "statistical significance", correct?) would be of a specific answer within a specific question within that survey? To be more specific, 100 individuals answered the question "how many servers does your organization run?" and the 100 responses are broken down as follows: "1-10" - 23 responses "11-20" - 20 responses "21-30" - 16 responses "31-40" - 19 responses "41-50" - 13 responses "51+" - 9 responses Does the confidence level/interval change for EACH response option (ie "1-10" or "11-20") within that response set? And if so, do you figure it out the same way you would for the overall response total (100) - meaning, do you assume "23" is your sample set (rather than "100") for a population of >10 million? It just seems something is off here with my thinking. Help? And a follow-up question. If I needed to slice the survey data BASED on the responses to the above question, how many more responses would I need in order to get EACH response option (ie "1-10", "11-20") into something statistically significant given the total market population? Many thanks for your time.


title: Genetic algorithm with constraint in R

text: I am trying to estimate a set of parameters using a genetic algorithm in R with the 'GA' package. So far I am doing something very simple (which works):               library(GA)     df <- data.frame(DM=c(1000,1500),                  c1=c(50,75),                  c2=c(90,105))          func <- function(x1, x2, theta) theta[1]*x1 + theta[2]*x2     fitnessL2 <- function(theta, x1, x2, y) {                   f <- -sum((y - func(x1, x2, theta))^2)                  }     GA2 <- ga(type = "real-valued", fitness = fitnessL2,               x1 = df$c1, x2 = df$c2, y = df$DM, names = c("a", "b"),                min = c(0, 0), max = c(100, 100))      which produces values for a,b. However, I want to constrain a and b not individually, but as something like a + b == 20. Can someone please let me know if this is possible? Thanks in advance.


title: On distance between parameters in Ridge regression

text: In ridge regression we know that ![](http://i.stack.imgur.com/3YeSy.jpg) as an estimate of $\beta$ and this gives the minimum sum of squares of the residuals: ![](http://i.stack.imgur.com/1vkIm.jpg) And we know that ![](http://i.stack.imgur.com/6rmkv.jpg) The question is how to demonstrate that ![](http://i.stack.imgur.com/UwgBB.jpg)


title: Running fixed-effects model stata

text: I'm running a fixed effects model. My independent variable is store presence, and one of my dependent variables, `i.county`, measures county fixed effects. I use `xi: regress store_presence i.county other_var othervar2` where `county` is the US county code, a string variable. But the regression output reports that some variables "have been omitted due to collinearity". What should I do to fix this error and capture county-fixed effects in my model?


title: Is it possible to convert a Rayleigh distribution into a Gaussian distribution?

text: ...and how might we do this? If possible, I am curious if outliers in the Rayleigh distributed data would also remain outliers in the new Gaussian distributed data. Thanks.


title: Regression problem

text: Regression was estimated using OLS. We get y=a0 + a1x1 + a2x2 + error. We know covariance matrix ∑ of our estimator. **1\. How to get confidence interval for a1/a2 ratio?** **2\. In what case would the confidence interval be bigger - if there is positive correlation between a1 and a2 or if there is no correlation?** I have some feeling it has to do with Cholesky's decomposition of ∑ . First one is only concerned with description of the procedure. I got it at my exams today and had no idea how to suss it out ;) any ideas?


title: What is the difference between the 'Pivot table' and 'Contingency table'?

text: Is there any difference? Or it they are basically the same? Just excel users call it 'pivot table' and Statisticians call it 'contingency table'? Thanks,


title: How to Calculate silhouette coefficient in SPSS for clustered data set?

text: I am having a pre clustered dataset with data and the action cluster identified for it using a custom clustering method. I am looking to calculate silhouette coefficient on this clustered dataset using SPSS to determine the quality of clusters created; any idea how i can do that?


title: How to get summary statistics from "resistant regression" - lqs - in R?

text: I've been using the `lm` function in R to do demand modeling (tons of steel to be predicted by various economic indicators). I used $R^2$ and $F$ to report on the strength of the model. However, when I use the R function `lqs` ("resistant regression") and then type in `summary(model_name)` I do not get any statistics that I can use to report on the strength of the regression model. Any suggestions? EDIT: Thanks for your quick response. I don't have a problem with lqs(). The problem is that when I type in summary(Model) I do not get any goodness of fit information (e.g., adjusted R squared) as I do when I enter summary(x) where X is a model created using the lm function. I'd like to have something to show the strength of the model. I"m using MASS. See below. Regards, Bill Yarberry library(MASS)               M10 = lqs(agri ~ p12+p1+p11+p5+p8+p6+p25+p50+p35, data = agri_data2) summary(M10) Length Class Mode     crit 1 -none- numeric     sing 1 -none- character coefficients 10 -none- numeric     bestone 10 -none- numeric     fitted.values 103 -none- numeric     residuals 103 -none- numeric     scale 2 -none- numeric     terms 3 terms call     call 3 -none- call     xlevels 0 -none- list     model 10 data.frame list


title: Estimating occupancy rates from arrival rates

text: Let's say the number of hotel rooms in a city is X. We know the arrival rate of visitors every day. We don't know the current occupancy. Is it possible to   1. estimate the occupancy rate over all hotel rooms (assume 1 visitor / room)   2. estimate the distribution of stay duration (assume every visitor stays at least 1 day) On the flip side, if we knew just the occupancy rate, can we estimate the total number of hotel rooms from the arrival data?


title: How to search for a statistical procedure in R?

text: Is there an R package, website, or command that will allow one to search for a specific statistical procedure they desire? For instance, if I wanted to find a package that had the Box-Cox Transformation, the website/package/command might return "MASS" and refer me to the `boxcox()` function. It is fairly straightforward with something like Box-Cox, but I was hoping it would allow me to find more difficult procedures or search by what the function does ("Concatenating columns to a data frame" might turn up `cbind()`). Does something like this exist?


title: Comparing the rate of admitting to not quitting smoking between two groups?

text: One group of patients was interviewed by doctors and 2% patients admitted to smoking. Another group of patients was interviewed by nurses and their admittance rate to smoking was 10%. How would I show that the difference in rates is statistically significant?


title: Logrank test for trend (proportional hazards)

text: How one can create a logrank test for trend and does it differ from normal logrank test? Any suggestions or literature? Maybe some R examples and functions?


title: precision recall breakeven point

text: There are two popular measures that aggregate $Precision$ and $Recall$, there are $F1$ and $\text{Precision Recall Breakeven point}$. $F1$ can be calculated easily by formula, but how to calculate $\text{breakeven point}$? I have some experiment and get all four values: true positive, false positive, true negative, false negative, so now I can calculate Precision and Recall but how to find breakeven point in this case.


title: Interpreting binary logit model?

text: I have some output from a binary logit regression and I have some difficulties interpreting it properly. Particularly I would like to know whether the coefficient of “male” tells me about the probability of outcomes 0 and 1 and what role the intecept plays in that example. Dependent Variable is political preference.                           "Estimate" "Std. err." "z-score"     (Intercept) −5.183     1.094       −4.738     Male         0.652     0.318        2.053


title: What the dimension of an exponential family tell us about that family?

text: In Wikipedia it is stated that: > A vector exponential family is said to be curved if the dimension of $$ > {\boldsymbol \theta} = \left (\theta_1, \theta_2, \ldots, \theta_d \right > )^T$$ is less than the dimension of the vector $$ {\boldsymbol > \eta}(\boldsymbol \theta) = \left (\eta_1(\boldsymbol \theta), > \eta_2(\boldsymbol \theta), \ldots, \eta_s(\boldsymbol \theta) \right )^T.$$ with > $$ f_X(x|\boldsymbol \theta) = h(x) g(\boldsymbol \theta) > \exp\Big(\boldsymbol\eta({\boldsymbol \theta}) \cdot \mathbf{T}(x)\Big) $$ but what happens when the dimension of $\boldsymbol \eta$ is smaller than the one of $\boldsymbol \theta$ ?


title: Response-surface methodology vs. 3-level full-factorial

text: What is the difference between a Box-Behnken design, a central composite design and 3-level full-factorial design? What are the advantages or disadvantages of one design over the others?


title: bounds of wavelet coefficients of function

text: Let $A_k^p:=\\{f:\, \|(ix)^k\hat{f}(x)\|_p\leq 1\\}, k \in Z_+, p \in (1, \infty)$. I am wondering what the application of the lower and upper bounds of wavelet coefficients of the function on the function space $A_k^p$? I.e., if one would know upper and lower bounds for this coefficients, would it be helpful for some questions? If yes, what kind of questions? Thank you.


title: LDA vs. perceptron

text: I am trying to get a feel for how LDA 'fits' within other supervised learning techniques. I have already read some of the LDA-esque posts on here about LDA. I am already familiar with the perceptron, but just learning LDA now. How does LDA 'fit' into the family of supervised learning algorithms? What might be its drawbacks vs those other methods, and what might it be used better for? Why use LDA, when one could just use, say, the perceptron for instance?


title: A test of whether a sample has a specified distribution

text: Suppose we have an iid sample $X_1, \dots, X_n$. We want to test if it has a standard normal distribution. One test statistic I learned from Brockwell and Davis's Introduction to Time Series and Forecasting is to   * first find an interval $(-b, b)$ which has a probability $95\%$ under the standard normal distribution,   * then compute the ratio $r$ of the sample points falling into $(-b,b)$. The book continues to say if $r$ is not equal to $95\%$, then reject the null. But I think it is not right. What I think should be done instead is to find the distribution of $r$ under the null (and then find the rejection region of $r$). but I am not sure how to do that? This test should have very low power,isn't it? Because the test statistic only captures a little information about the sample's distribution? Thanks and regards! PS: in case my understanding is incorrect, the original text from the book says: ![enter image description here](http://i.stack.imgur.com/Xo7rS.png) Note that the book tells how to test if a sample is iid, but its underlying idea is to test if the sample autocorrelations is iid with distribution $N(0,1/n)$, which is the same question as in my post.


title: Group categorical variables in glmnet

text: Consider the following fit: `fit3a=glmnet(x,g4,family="multinomial",type.multinomial="grouped")` How do I indicate which columns in `x` are categorical/multinomial? Is there an option to specify the index of the grouped variables? The documentation describes the option `type.multinomial` as follows: > If "grouped" then a grouped lasso penalty is used on the multinomial > coefﬁcients for a variable. This ensures they are all in our out together. > The default is "ungrouped".


title: Why is GARCH(1,1) = ARCH(infinity)?

text: How do you prove that GARCH(1,1) = ARCH(infinity)? Need some guidance on how to start on this.


title: Mathematical definition of causality

text: Let $Y$ and $X$ be random variables. $E(Y|X)$ is the conditional mean of $Y$ given $X$. We say $Y$ is not causally related to $X$ if $E(Y|X)$ does not depend on $X$, which implies it is equal to $E(Y)$. Now, let's go along with this definition of causality for a second. By the law of iterated expectations, $E(XE(Y|X)) = E(E(XY|X)) = E(XY)$. This means that if $E(Y|X)$ does not depend on $X$, if it is equal to $E(Y)$, then $E(X)E(Y) = E(XY)$. In other words: If $X$ and $Y$ are not causally related, then $X$ and $Y$ are uncorrelated! - This makes no sense and I know this must be wrong. Have I defined causality incorrectly? What have I done wrong? In econometrics we generally assume $E(Y|X) = b_0 + b_1X$. So $E(Y|X) = E(Y)$ is equivalent to $b_1 = 0$. The logic applies in this specific scenario too.


title: Generalized Moran's I test for multinomial (polychotomous) logit / probit model

text: Generalized Moran's I test is suggested by Kelejian and Prucha (2001). But, as far as I know, there is no empirical work using the test for multinomial discrete choice model. I am looking for a scripting that handles generalized Moran's I test in multinomial logit / probit models. Moreover, I would be interesting in understanding (4.16), (4.17) in Kelejian and Prucha (2001). If someone understands generalized Moran's I test, I would appreciate it someone could point that out to me. Reference: Kelejian, H. H., Prucha, I. R.: On the asymptotic distribution of the Moran I test statistic with applications, Journal of Econometrics, 104 (2), 219-257, 2001. (pdf)


title: How big are regularization parameters values?

text: I wanted to know how big are the regularization parameter values for ridge or lasso. I have seen most of the places generally using values like 0.1 or 0.01 but in some of my experiments the cross validation gives regularization parameters with values something like 8.6 and 172. Are these values fine or am I doing something wrong?


title: Reverse PCA with nonsquare matrix

text: I am working with PCA to do dimensionality reduction to a set of data. I have 1600 data points with 36 variables and I want to have a matrix with a new data set with 6 principal components. I managed to do it in matlab, so I have:               xtrain (normalized) <1600x36>      and I write:               coeff=pca(xtrain,'NumComponents',6)       which gives me a <36x6> matrix So far so good. Now I want to know how to reverse the process. If I'm given a data set of how do I use the coeff matrix to pass it to a "n x 36" representation? I would like to do this in matlab but a simple explanation can help me.


title: Using Random Number Table

text: In a random number table, the random numbers are given as first thousand, second thousand, third thousand, and so on. Every thousand has 25 rows and 1 to 4, 5 to 8 up to 37–40 columns. Suppose I have to draw a sample of size 40 from a population of 10000, and I am using the 1 to 5 column for selecting the random numbers. When the 1 to 5 column of first thousand has been finished, but I have not arrived at my desired sample size yet, what is the process to select the random numbers? Will I need to go to the `column 6 to 10` of the `first thousand`? Or go to the `column 1 to 5` of the `second thousand`? Another table here is the 1st page. If I decide to use the 1 to 5 number columns, when the 1 to 5 number columns of the 1st page are finished, will I move to the 1 to 5 number columns of the 2nd page or stay on the 1st page using 6 to 10 columns?


title: Drug affects individuals but ANOVA not significant: references?

text: Let's suppose that I measure the effect of a drug on blood pressure in several people **many times**. For some people, the drug increases blood pressure significantly and for others decreases it significantly. So the drug is affecting the individuals. Because the direction of the effect is different for each person, it is possible that if I perform a global statistical analysis such as ANOVA, the drug shows not significance. So, it is incorrect to conclude that the drug does not affect people by using the ANOVA. I would like to know references that talk about this problem for being able to write about it using more correct terminology. More specifically, this is what I have. I collected a dose-response curve for each subject with and without drug, fit a cumulative normal distribution to each condition and estimate the mean and confidence intervals by bootstrapping. So I collect many data for each subject. For each subject the mean with and without drug is significantly different(no overlap of confidence intervals) but, for some observers the drug increases blood pressure while for others decreases blood pressure. So if I use the means to perform an ANOVA, the effect of the drug is not significant. So, it seems to me that report only the ANOVA without caring about the individual effects is not correct, but this is sometimes done in my field. So I was looking for some references that pointed this out. EDIT: I think I chose the wrong example using drugs and blood pressure. So I asked another question using an example closer to my real problem fitting a model of each subject and perform a t-test over the estimated parameters


title: Getting probabilities over 1 in positive and unlabeled learning

text: I have a question regarding PU-Learning, which deals with learning from positive-labeled (no labeled negative examples) and positive/negative- unlabeled data. Particularly, my question is about the paper Learning classifiers from only Positive and Unlabeled data. This paper converts a non-traditional classifier $g(x)$ which is learned from labeled/unlabeled datasets and outputs the probability of an example of being labeled, into a traditional classifier $f(x)$ which outputs the probability of an example to be positive. I have implemented the first proposed approach (section 2, "learning a traditional classifier from nontraditional input"). However, the probability $f(x)=g(x)/p(s=1|y=1)$ becomes greater than 1 for several examples (this should not be greater than 1, since it is a probability) For learning the non-traditional classifier $g(x)$, I am using a non- traditional dataset composed of   1. 100,000 labeled examples randomly chosen from the whole labeled data, and    2. 100,000 unlabeled examples randomly chosen from the whole unlabeled data.  And for estimating $p(s=1|y=1)$--with the first proposed estimator--I am using a validation dataset composed of 30,000 labeled instances randomly chosen from a separate labeled dataset. According to the paper, the probability $f(x)$ is guaranteed to be well-formed if (1) and (2) overlap in the example space. In my case, (1) and (2) have a 20% of overlapped examples, but I am still getting non well-formed probabilities (>1). How could I achieve a well-defined probability for $f(x)$? **EDIT 8/21**   1. According to the paper, $g(x)$ must be a classifier that produces correct probabilities as its output, like Logistic Regression, or a calibrated classifier such as Naive Bayes/SVM, in order to get the approach to work. I am simply using Maximum Entropy classifier (from the NLTK package), also known as Logistic Regression. Therefore I guess there should not be any problem with this.   2. Although the paper states that the labeled (1) and unlabeled (2) datasets for training $g(x)$ are "samples from overlapping regions in feature space", such datasets seem to be disjoint (see experiment from Section 5, P and U are disjoint). I have tried with disjoint labeled and unlabeled datasets as well, however I am still getting probabilities over one for $f(x)$.   3. Balancing datasets does not make any difference, neither.    4. I have tried the second estimator proposed(e2). Still getting probabilities over 1 for $f(x)$.  **EDIT 8/29**   1. In section 5, which describes an example with real-world data, we have a set P of labeled+positive examples from the database TCDB, and a set U of unlabeled examples randomly sampled from the database SwissProt, being P and U disjoint. Then they use P,U to learn a non traditional classifier. I think that this approach can be applied to my problem, since I also have a set of positive examples, and a set of unlabeled examples randomly sampled. What do you think ?   2. I have tried the third estimator proposed(e3). Since it takes the maximum probability of the examples from the positive-labeled set (in my case around 0.98), I am no longer getting probabilities over one. However, since this estimate is only based in one example and not the average between all examples, therefore this does not look like a good estimator to me. Any thoughts regarding the validity of this estimator? Thanks


title: How to compare median survival between groups?

text: I'm looking into median survival using Kaplan-Meier in different states for a type of cancer. There are quite big differences between the states. How can i compare the median survival between all the states and determine which ones are significantly different from the mean median survival all across the country?


title: Is it o.k. to bootstrap sample of a table from its proportions - and how to do so (in R)?

text: I wish to create a confidence interval for a statistic calculated on a table (let's say the p.value of chisqr.test). For this, I might sample a bootstrap sample with the same size, from a table with the proportions of the table I have. And my questions are:   1. Is this a correct procedure to use, or is there some catch here that I might be missing?   2. How can a table be "sampled" in R? (my first instinct is to "open up" the table using melt, then expend the table to rows, sample them, and fold the results back to a table. Is there a better way?) p.s: I am aware of the "simulate.p.value" switch in chisqr.test - my question is more general - does it make sense to use bootstrap in such a way. Sample table if someone wishes to show something:               (x <- matrix(c(12, 5, 7, 7), ncol = 2))     prop.table(x)          # undoing a table     as.data.frame(as.table(x)) # but how do we open up the rows now? - **update**: this was answered in the comments thanks to chl.      Thanks.


title: Interrater reliability for events in a time series with uncertainty about event time

text: I have multiple independent coders who are trying to identify events in a time series -- in this case, watching video of face-to-face conversation and looking for particular nonverbal behaviors (e.g., head nods) and coding the time and category of each event. This data could reasonable be treated as a discrete-time series with a high sampling rate (30 frames/second) or as a continuous-time series, whichever is easier to work with. I'd like to compute some measure of inter-rater reliability, but I expect there to be some uncertainty in _when_ events occurred; that is, I expect that one coder might, for example, code that a particular movement began a quarter second later than other coders thought it started. These are rare events, if that helps; typically at least several seconds (hundreds of video frames) between events. Is there a good way of assessing inter-rater reliability that looks at both of these kinds of agreement and disagreement: (1) do raters agree on _what_ event occurred (if any), and (2) do they agree on _when_ it occurred? The second is important to me because I'm interested in looking at the timing of these events relative to other things happening in the conversation, like what people are saying. Standard practice in my field seems to be to divide things up into time slices, say 1/4 of a second or so, aggregate the events each coder reported per time slice, then compute Cohen's kappa or some similar measure. But the choice of slice duration is ad-hoc, and I don't get a good idea of uncertainty in time of events. Best thought I have so far is that I could compute some kind of reliability curve; something like kappa as a function of the size of the window within which I consider two events as being coded at the same time. I'm not really sure where to go from there, though...


title: classic non-monotone missing data sets

text: I am developing a method to handle non-monotone missing data. What are some classic non-monotonic missing data sets, especially in the area of public health?


title: Why my tree model only took 1 variable?

text: I have the following data frame:               $ Status           : Factor w/ 2 levels "Active","Withdrawn": 1 1 1 1 1 1 1 1 1 1 ...     $ Age              : num  41 45 43 24 34 29 36 33 28 24 ...     $ Tenure           : num  6 8 20 8 8 9 2 16 47 49 ...     $ Current.Hire.Date: Date, format: "2002-01-22" "2000-03-31" "1989-06-12" ...     $ Gender           : Factor w/ 2 levels "Female","Male": 2 2 1 2 2 1 2 1 2 1 ...     $ Function         : Factor w/ 10 levels "","Finance","HR",..: 10 2 10 8 9 10 6 8 10 9 ...     $ Grade            : Factor w/ 21 levels "11","12","13",..: 15 20 10 13 19 17 14 17 17 14 ...     $ Retiree          : Factor w/ 2 levels "No","Yes": 2 2 2 1 1 1 1 2 1 1 ...     $ Race             : Factor w/ 5 levels "Asian","Black",..: 5 5 5 5 1 5 5 5 5 5 ...      I want to use decision tree to predict the binary outcome which is `Status`. I used the following script:               > tree.turnover_nonprod = tree(Status~Age+Tenure+Race+Gender+Function+Grade+Retiree, data = turnover_nonprod)     > summary(tree.turnover_nonprod)     Classification tree:     tree(formula = Status ~ Age + Tenure + Race + Gender + Function +      Grade + Retiree, data = turnover_nonprod)     Variables actually used in tree construction:     [1] "Function"     Number of terminal nodes:  2      Residual mean deviance:  1.352 = 6774 / 5010      Misclassification error rate: 0.4419 = 2215 / 5012       why the model only took `Function` but not the rest of the variables? I am sure I am doing something wrong here but don't know what. Thanks,


title: Alternative Specific Variables in R

text: I am building a discrete choice model (rail and auto). I have the cost of a trip for each mode (rail versus auto). The utility equations I am building are denoted below: \begin{align*} V_\text{auto} &= \beta_1\cdot\mbox{IVTT}_\text{auto}+ \beta_2\cdot\mbox{OVTT}_\text{auto} \\\ V_\text{rail} &= \beta_0+ \beta_1 \cdot\mbox{IVTT}_\text{rail} + \beta_2 \cdot\mbox{OVTT}_\text{rail}+ \beta_3×\mbox{Cost}_\text{rail} \end{align*} The code I have is:               data_mlogit<-mlogit.data(newdata,shape="wide",choice="choice",                              alt.levels=c("auto","rail"),varying=c(6:8,10:12))     model <- mlogit(choice ~ ivtt+ovtt+cost|0|0, data_mlogit)       How can I make the code understand to use the cost rail variable only without using the cost auto variable? I don't want the cost variable (auto) to be part of the utility equation for the auto mode.


title: Can the mutual information of a "cell" be negative?

text: Please forgive me if this is not the right Stack Exchange (and for inventing terms). For discrete random variables X and Y, the mutual information of X and Y can be defined as follows: $I(X;Y) = \sum_{y \in Y} \sum_{x \in X} p(x,y) \log{ \left( \frac{p(x,y)}{p_1(x)\,p_2(y)} \right) }, \,\\!$ I will define the mutual information of a "cell" $x_0$ to be: $CI(x_0,Y) = \sum_{y \in Y} p(x_0,y) \log{ \left( \frac{p(x_0,y)}{p_1(x_0)\,p_2(y)} \right) }, \,\\!$ I'm not sure if this quantity goes by another name. Essentially I'm restricting focus to a single state of variable X (and then the full MI can be calculated by summing all the cell MIs). My question: is it guaranteed that $CI(x_0,Y) \ge 0$? We know $I(X;Y)\ge0$ and we know that the pointwise mutual information can be negative. I feel like CI should be nonnegative and that I might be missing some obvious proof.


title: Long term weather data - how to get there most painlessly

text: I'd like to get daily or weekly weather data (temperature and precipitation) from specific measurement stations around the USA. I don't need anything extrapolated to my specific coordinate, a reasonably close measurement station is enough. Preferably I would like to have that data as far as possible back in time. At least the seventies. I would like to work with the data in matlab (and/or R). How do I get there the easiest? I see Wunderground/NOAA has some sort of things. But nothing seems very easy, wunderground requires you to download thousands of text files and figure out how to make one timeseries.


title: Cost function of neural network is NON-CONVEX?

text: The cost function of neural network is $J(W,b)$, and it is claimed to be **non-convex**. I don't quite understand why it's that way, since as I see that it's quite similar to the cost function of Logistic Regression, right? If it is non-convex, so the 2nd order derivative $\frac{\partial J}{\partial W} < 0$, right?


title: What are the major differences between the parameter estimation of a simple linear regression vs multiple?

text: This is the standard technique that I have used to calculate the parameter of a single linear regression problem with one independent variable. How does the parameter estimation problem differ in multiple linear regression with 'k' independent variables. I assume I'd treat each independent variable in a multiple linear regression problem separately, which would then give me 'k' different sub simple linear regression problems, and then I'd add them parameters to give me my final model. I am requesting to be corrected or improved upon. ![enter image description here](http://i.stack.imgur.com/lYevl.gif)


title: Warning message in auto.arima

text: I am using `auto.arima()` for prediction, and getting the following warning message. I want to know if I can ignore this warning message or if I should be worried. `Warning message:`   `In auto.arima(forecast_data_ts) :`   `Unable to fit final model using maximum likelihood. AIC value approximated`


title: Cannot replicate odds ratio

text: I am trying to replicate some odds ratio calculations which are reported in this paper. For example in the paper on page 3 in the Strome et al. data for tonsil: it states there are 52 cases and 48 controls. HPV % cases = 40% and controls = 6%. Converting the percentages to raw figures, I get the following 2×2 table:                              Case  Control     HPV positive     21        3        HPV negative     31       45      Therefore the odds ratio should $=\frac{21\cdot45}{31\cdot3} = 10.16$. However the paper reports this particular odds ratio as 18.2 on the following page. Does anyone know why this is? Are they using a different formula?


title: Can I use the Williams type contrast from `multcomp` package to perform a Williams test(Williams, 1971)?

text: Using Williams type contrast in the multiple comparison test with `multcomp` package will give different test statistics and test p-values for the same dataset and same hypothesis. I would like to know the difference between the two procedures. Can I use the Williams type contrast from `multcomp` package to perform a Williams test (Williams, 1971)? I see lots people do it, but I am not sure if this is correct. An example is given below. The dataset is `warpbreaks`. The result from `glht` is:               Simultaneous Tests for General Linear Hypotheses     Multiple Comparisons of Means: Williams Contrasts          Fit: aov(formula = breaks ~ tension, data = warpbreaks)          Linear Hypotheses:              Estimate Std. Error t value Pr(>|t|)         C 1 == 0   -14.72       3.96  -3.718 0.000827 ***     C 2 == 0   -12.36       3.43  -3.604 0.001166 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Adjusted p values reported -- single-step method)      The contrast matrix for the test is:                        (Intercept) tensionM tensionH     C 1           0      0.0      1.0     C 2           0      0.5      0.5     attr(,"type")     [1] "Williams"      And the result from applying the original Williams' test (1971) is:                 Dose     Mean       MLE     Tbar pvaltab WillTab.95%     1    L 36.38889 -36.38889 0.000000      NA          NA     2    M 26.38889 -26.38889 2.349153     <5%    1.673529     3    H 21.66667 -21.66667 3.458476     <1%    1.753529       WillTab.99% critical_tab.95 critical_tab.99     1          NA              NA              NA     2    2.400588       -29.26492       -26.16994     3    2.444118       -28.92437       -25.98464      The code I used is here, and `EcotoxTests` is a package I wrote myself for performing a Williams test.               {r}     library(multcomp)     amod <- aov(breaks ~ tension, data = warpbreaks)     g <- glht(amod, linfct = mcp(tension = "Williams"))     summary(g)     require(EcotoxTests)     dat1 <- data.frame(dose=warpbreaks$tension,breaks=warpbreaks$breaks,nums=warpbreaks$breaks)     calcWilliams(dat1)


title: Joint vs. marginal prediction intervals for path forecasts (with k-family wise error rate)

text: I am trying to become comfortable with the bootstrapping of joint prediction regions described in this paper: http://www.nccr-finrisk.uzh.ch/media/pdf/wp/WP748_A3.pdf This calculates the prediction intervals, $I_j$ for all time points $j$, such that all datapoints of a forecast path lie within this region with probability $1-\alpha$, i.e.: $P(y_j \in I_j \forall j) = 1-\alpha$ The authors of the paper extent the idea to the "k-family wise error rate" which narrows the intervals by requiring that at least k datapoints of the forecast path lie outside the region: $P(\sum_1^N1[y_j \notin I_j] >= k) = 1-\alpha$ where $1[.]$ is the indicator function giving 1 when the condition matches and 0 when it is false. N is the number of forecast data points. By playing around, I discovered that for $k=N/2$ the k-family error rate joint prediction intervals seem to be identical to the marginal prediction intervals, i.e.: The following seem to give identical $I_j$ $P(\sum_1^N1[y_j \notin I_j] >= N/2) = 1-\alpha$ $P(y_j \in I_j) = 1-\alpha$ Is this pure coincidence? How can this be shown formally?


title: Possible classification techniques to use when each feature is a probability distribution

text: I am working with some data where the features have a _temporal_ aspect (e.g. how often does a feature occur between $t_{begin}$ and $t_{end}$). I am trying to build a binary classifier for this data. The problem, however, is that each feature _is_ a probability distribution. Normally, something like a SVM classifier would work: each object represented by a vector, where the $i^{th}$ index of the vector is a number (e.g. tf-idf, etc.). I have no idea how to proceed where each object is represented by a set of probability distributions. The simple vector representation will clearly not be suitable. I've searched far and wide, but haven't found anything that suits this kind of data. Any pointers/ideas will be greatly appreciated. Also, I wanted to add that I have looked into distances between distribution and density functions (e.g. Wasserstein metric, Kolmogorov-Smirnov statistic ), but I don't see how computing these metrics will help in the classification.


title: R: One period our cross validation with time series

text: I have quarterly data with one causal variable (X) and one dependent variable (Y). 30 such observations. I have the X variable for a quarter, and I'm seeking to predict that quarter's Y. The structure of my time series requires one seasonal dummy every 4th quarter - holiday shopping spiking in Q4. I want to know what the error rate is for a one period out forecast (using X), for the last 5-10 periods. How could I implement this in R?


title: How should I analyze repeated-measures individual differences experiments?

text: I do quasi-experimental individual differences psychology research. I examine how people who differ in a cognitive ability (that I measure) perform on another task that always at least involves within-subject manipulations (and sometimes between-subject)–DVs are usually response time and/or accuracy. For this question I'd like to focus on response times (let's assume they are normally distributed). I then infer from the ability-task relations what it means theoretically for the cognitive ability. The nature of this work is correlational and involves repeated measures where each subject completes many task trials(mostly I’m not interested in changes over time, just the overall difference). Researchers in my field often create categorical variables from the cognitive ability scores and compare the upper and lower quartiles of the distribution with a repeated-measures ANOVA. Because the cognitive ability is measured continuously, I am looking for an analytic strategy that treats the cognitive ability in this way. I’ve recently been investigating mixed-models, thinking that I can treat the people as a random effect grouping variable and the cognitive ability score as a random effect nested within people. I would like to examine interactions between this nested random effect (cognitive ability) and the fixed effects of the experimental treatments by doing model comparisons. Does this seem like a reasonable analytic strategy? Am I thinking about this right? What are some other ways (the simpler, the better) that I can take advantage of repeated measures -remove experimental subject variance- while also maintaining the cognitive ability measure as a continuous measure? Any citations or examples in R are appreciated. In a typical experiment, I would expect there to be anywhere from 1-3 categorical IVs with 2-4 levels(measured by multiple trials), and 1 continuous IV(cognitive ability). The exact nature of the categorical variables changes from study to study. The designs are fully crossed.


title: For count data from a survey, do variance corrections for survey design imply that the Poisson distribution will not accurately model the counts?

text: I have categorical count data that comes from a complex survey. Each unit of analysis in the survey (`household`, `individual`, etc.) is put into one and only category per dimension, ranging from 2 to 20 categories per dimension. The number of dimensions exceeds 50. The survey reports a variance estimate corrected based on the clustering and stratification in the survey design for each count. I have three closely related questions:   1. Does the fact that the variances for the cell counts depend on something other than and additional to the value of the count or the underlying cell probability imply that these cells are not distributed according to the Poisson distribution?   2. If not, why?   3. If so, can we say what distribution is probably the correct one, or what are good candidates?    I was thinking negative binomial, but that just gives you one more coefficient, while the survey results include a distinct variance correction for each cell.


title: Shapiro test and transformation

text: I am playing with `shapiro.test` from R and checking for non-normality of error variance.               > shapiro.test(residuals(fit))                  Shapiro-Wilk normality test          data:  residuals(fit)      W = 0.9525, p-value = 0.0003303          > shapiro.test(residuals( lm( sqrt(V1)~V2 ,data=market) ) )                  Shapiro-Wilk normality test          data:  residuals(lm(sqrt(V1) ~ V2, data = market))      W = 0.8895, p-value = 5.89e-08          > shapiro.test(residuals( lm( log(V1)~V2 ,data=market) ) )                  Shapiro-Wilk normality test          data:  residuals(lm(log(V1) ~ V2, data = market))      W = 0.7698, p-value = 1.95e-12          > shapiro.test(residuals( lm( 1/(V1) ~ V2 ,data=market) ) )                  Shapiro-Wilk normality test          data:  residuals(lm(1/(V1) ~ V2, data = market))      W = 0.3954, p-value < 2.2e-16      The p-value for normality test is <0.001. So I did a transformation on `V1`, `log(Y)`, `inverse(Y)` and `sqrt(Y)` but their p-values gets even smaller. Does this mean that these transformations don't work? I also did a `boxcox` transformation, and I get p-value of 0.3. So in that case, `boxcox` is the best remedy only in this case? Or should I just use `boxcox` in the future for remedy of non-normality? The rest of the usual transformation methods are useless?


title: Should an SVM grid search show a high-accuracy region with low accuracies around?

text: I have 12 positive training sets (cancer cells treated with drugs with each of 12 different mechanisms of action). For each of these positive training sets, I would like to train a support-vector machine to distinguish it from a negative set of equal size sampled from the experiment. Each set has between 1000 and 6000 cells, and there are 476 features (image features) of each cell, each scaled linearly to [0, 1]. I use LIBSVM and the Gaussian RGB kernel. Using five-fold crossvalidation, I have done a grid search for log₂ C ∈ [-5, 15] and log₂ ɣ ∈ [-15, 3]. The results are as follows: ![Results of grid search](http://i.stack.imgur.com/8GiiB.png) I was disappointed that there is not a single set of parameters that give high accuracies for all 12 classification problems. I was also surprised that the grids do not generally show a high-accuracy region surrounded by lower accuracies. Does this just mean that I need to expand the search parameter space, or is the grid search an indication that something else is wrong?


title: Comment on this method of sampling from any distribution?

text: Say I have a distribution, either described by a probability density function that is integrable and continuous, or by a set of discrete probabilities over a finite set of symbols. I want to efficiently generate samples from this distribution. I don't care about references to other algorithms, I want to know if this idea is essentially correct, or how to improve it. Here is how I currently do it for the discrete case :   1. form the cumulative probability function ( defined at the ith symbol as the sum of all probabilities from the first symbol to the ith symbol ).    2. Then, there is an interval, [0,max(CDF)] and it is trivial to output a IID random number anywhere in that interval using a regular PRNG.    3. Since every symbol is now associated with an interval, we output the symbol in whose interval the IID number landed. It seems completely obvious to me that the output sequence will have the characteristics of the given discrete distribution, and is optimal up to the original PRNG. I want now to extend to the continuous case (this may be the Ziggurat algorithm, and if so, boy am I happy I jumped straight to that). But can we not simply integrate the PDF to get a CDF, again form the interval [0,max(CDF)] and again shoot IID randoms at this interval, and again, and this is the tricky part for me, _associate_ **_each_** _sample symbol_ (of which there are infinitely many) _with an_ **_interval_**. Then output the sample symbol in whose interval the IID number landed. So, my question is : how can I associate the sample symbols, with intervals? Presumably by discretizing (putting into buckets) the sample symbol interval, and giving anything in a bucket the same probability. The size of buckets could be worked out so as to, over the intended output stream length, not cause any detectable deviation from the underlying continuous distribution. Also, I am doing this in JavaScript. It's a DIY so I am not interested in trusting or using someone else's code. But nor am I interested in spending 2 days on it. I'm asking here as I am not an expert in distributions, and I want to have the ideas corrected before I try to implement this.


title: Which type of regression analysis should I use for consumption pattern data?

text: I conducted a survey asking people of various income groups about their annual expenditure under various consumption brackets such as food, clothing, housing etc. I then calculated the amount spent by people of different income brackets in these categories as a percentage of their total income. Now, to plot this data and find the best fit line which type of regression analysis should I use. Also are there some other good statistical tools that I can use for exploiting this data?


title: Pitfalls in doing GLMM or interpreting GLMM results

text: I know this is a general question, but it would benefit CVers to learn about common pitfalls, things to consider and caveats on how to do rigorous GLMM modeling. Could those with experience on GLMMs share any insights they might have?


title: Can you use heteroskedastic time series variables within a regression model?

text: We are working on a multivariate linear regression model. Our objective is to forecast the quarterly % growth in mortgage loans outstanding. The independent variables are: 1) Dow Jones level. 2) % change in Dow Jones over past quarter. 3) Case Shiller housing price index. 4) % change in Case Shiller housing price index over past quarter. In a stepwise regression process, all above variables were selected. And, variables 1) and 3) were surprisingly significant. But, that is when used in combination. Somehow, there is something about the difference in those two indeces that does partly explain the % change in mortgage loans outstanding. For my part, I find variables 1) and 3) problematic. This is because I believe they are heteroskedastic. And, their respective coefficients confidence interval are therefore unreliable. I also think they may have multicollinearity issues with their related variables based on % change. They may also cause some autocorrelation issues. However, at first it seems some of my concerns may be overstated. After graphing the residuals of the whole model they seem OK. They don't trend upward. So, it appears the heteroskedasticity is not an issue for the whole model. Multicollinearity is not too bad. The variable with the highest VIF is around 5 much lower than the usual threshold of 10. Nevertheless, I am still somewhat concerned that even though the whole model seems OK; the specific regression coefficients of the mentioned variables may not be (or more specifically the related confidence intervals).


title: What test to use for overlapping categories

text: i am analyzing a surveillance tool where people could be collected by one or both methods. So it kinda looks like a Venn diagram. What test would I use to test if there's significant differences by the methods for different characteristics? The only statistical package I have access to is SAS so please recommend something that can be done in that. Thanks.


title: How to interpret the frailty terms in a coxph.penal model?

text: I have the following `coxph` model:               Surv(Y) ~ genotype + treatment + age + frailty(batch) + genotype:frailty(batch)      `batch` is a random effect (two replicates of the experiment were done with different subjects). Events can only happen once per subject. The summary output looks like this:                                         coef    se(coef) se2     Chisq DF   p           genotypeTG                -0.2504 0.44224  0.43550  0.32 1.00 5.7e-01     treatmentX                -0.6695 0.16379  0.16370 16.71 1.00 4.4e-05     age                        0.0135 0.00839  0.00829  2.59 1.00 1.1e-01     frailty(batch)                                      1.39 0.94 2.2e-01     genotypeTG:frailty(batch) -1.3479 0.40020  0.39353 11.34 1.00 7.6e-04     genotypeTG:frailty(batch)         0.00000  0.00000       1.00              My question is: given that there is a significant interaction between genotype and batch, is the genotype coefficient a summary of the two batches, or is it just for the batch the model is using as the reference group? If the former, how do I extract the batch-specific estimates? If the latter, how do I construct the overall coefficient and its SE? Thanks. * * * Update: by playing around with `predict()` and `coef()` I found:               > coef(myfit)                     genotypeTG              treatmentX                        age               gamma:batch1               gamma:batch2 genotypeTG:frailty(batch)1     genotypeTG:frailty(batch)2                    -0.25044857             -0.66945067                 0.01350297                -0.17817716                 0.15117985                -1.34789198                             NA           > cbind(mynewdata,predict(myfit,mynewdata))       genotype treatment  age  batch   predicted     1       WT         Z    0 batch1  0.64416111     2       TG         Z    0 batch1 -0.95417945     3       WT         X    0 batch1 -0.02528956     4       TG         X    0 batch1 -1.62363012     5       WT         Z    0 batch2  0.97351812     6       TG         Z    0 batch2  0.72306955     7       WT         X    0 batch2  0.30406745     8       TG         X    0 batch2  0.05361887      ...assuming no bugs in `predict.coxph.penal()`, this means that the first three coefficients are ONLY for `batch1`. For groups that have a level which interacts with the frailty term, I add the interaction coefficient for that batch in addition to the coefficient/s for those variables. To get `batch2` estimates, I do the above but also SUBTRACT the `batch1` frailty intercept, add the `batch2` frailty intercept, and then add other coefficients as appropriate. The interaction term for `batch2` is NA but apparently it doesn't get added. My question now reduces to, how do I get the overall coefficients and SEs (i.e. the equivalent of what `summary.lme()` would report if this were an `lme` fit)? Also, should I be concerned that one of the frailty interaction coefficients is NA and that the summary omits some of the frailty coefficients? * * * Or, am I creating unnecessary problems for myself by allowing an interaction between a frailty term and a fixed effect? My reasoning was that bidirectional `stepAIC()` identified this interaction as improving the fit of the model; this implies that the response of transgenic organisms may have been systematically different from that of normal ones between the two replicate experiments. Therefore, ignoring the interaction would amount to ignoring this source of error by treating the experiments as though only the baseline hazard differs between them. Is there a flaw in this reasoning somewhere?


title: Comparison between MDL and BIC

text: I'm currently studying Hidden Markov Models. There's a set of observations from which I need to determine the optimal number of states. After having found the maximum likelihood using Baum-Welch, I considered two model selection criteria for determining the optimal states. These are Minimum Description length (MDL) and Bayesian Inference criterion (BIC). However, with MDL, the number of states=2 whereas with BIC it's 4. Does this mean that MDL performs better than BIC?


title: Hierachical Predictors in a Regression

text: Note: Mainly this question pertains to predictions from a model. If the unit of analysis of a regression (or any predictive model really) is the individual retail store and these stores are organized into geographic regions, when does it make sense to include variables that describe the geographic region (i.e. a coarser grain than the unit of the analysis) versus just adding a factor (nominal variable) for each region? For example, adding a factor for each region would literally be a variable that takes on values 'A', 'B'......etc. Other descriptions of the region might be total number of households, per capita income. etc Given this set-up, there are two cases I can think of and I would like to know it it makes any sense to add the additional descriptions of the region (versus just the factor): 1) We don't expect the descriptions to change over the time we will use the model for prediction (e.g. We have 2013 total household information per region and we wont get an update until 2016). 2) We do expect the descriptions to change over the time we will use the model for prediction (e.g. We have 2013 total household information per region and we will get an update next month - or whenever we will next use the prediction).