title: Asymptotic normality and normalization wrt variance

text: Let $X_n, n \in \mathbb N$ be a sequence of random variables with finite variances. As $n \to \infty$, are the following two equivalent:   * $X_n \to N(0, \sigma^2)$ for some $\sigma^2 \in [0, \infty)$,   * $\frac{X_n}{\sqrt{Var(X_n)}} \to N(0,1)$? **Motivation of my question:** The asymptotic normality of MLE is usually given with its asymptotic variance being inverse of Fisher information under some regularity conditions: $$ \sqrt{n}\big(\hat\theta_\mathrm{mle} - \theta_0\big)\ \xrightarrow{d}\ \mathcal{N}(0,\,I^{-1}(\theta_0)). $$ All of Statistics by Wasserman, however states ![enter image description here](http://i.stack.imgur.com/vmuAA.png) I was wondering if the two results (or conclusion parts only) about asymptotic normality of MLE are equivalent? Thanks and regards!


title: If ϵ is uniformly distributed, then a linear probability model is appropriate? Can I find any Literature?

text: A latent variable model involving a binomial observed variable $Y$ can be constructed such that $Y$ is related to the latent variable $Y^*$ via $ Y = \begin{cases} 0, & \mbox{if }Y^*>0 \\\ 1, & \mbox{if }Y^*<0\. \end{cases} $ The latent variable $Y^*$ is then related to a set of regression variables $X$ by the model $Y^* = X\beta + \varepsilon$. This results in a binomial regression model. The variance of $\varepsilon$ can not be identified and when it is not of interest is often assumed to be equal to one. If $\varepsilon$ is normally distributed, then a probit is the appropriate model and if $\varepsilon$ is log-Weibull distributed, then a logit is appropriate. If $\varepsilon$ is uniformly distributed, then a linear probability model is appropriate.


title: Monte Carlo Integration Interval Probability

text: Use MC integration to estimate the probability that X * exp(X) < 2.5, assuming that X ~ Gamma(1.2,3.7)               data mcprobdata;     /* generate the sample */     call streaminit(23891);     count = 0;     do i=1 to 2000;         p = rand('uniform');         x = quantile('Gamma', p, 1.2, 3.7);         count = count + (x*exp(x) < 2.5);     end;          /* calculate the estimate */     t = count/2000;     se = sqrt(t*(1-t)/2000);          run;      Is my t variable correctly calculating this interval probability?


title: dealing with exponentials in python - infinities and overflows

text: In a machine learning algorithm that I'm using, I need to get the exponential values of something in one of the steps. This is the step that I'm dealing with right now: ![enter image description here](http://i.stack.imgur.com/f65Na.png) I've already got all the 1+g_j(X_i) etc etc calculated, there's no problem in that. Let's call it `calculated_value`. That's fine. The problem is, I am getting infinities and overflows when implementing this function in code. What's the solution for this? What would the logarithmic version of the same function, that I could use in place of the function above? In case it matters, β is going to be used as a coefficient to some variables, including partially to calculate the weights for a weighted regression solution.


title: Cummulative Mixed Model in R variable names

text: I'm trying to fit a cumulative link mixed model clmm() in Rstudio. I'm currently having issues with the diagnosing what is wrong with my model from the output I am getting. The output I got from my model changes the variable names. This only happens to dummy variables but not the scale variables. Does anyone know what this indicates, or if I should be worried about it? The names also changed the first time fitting the model. I fixed this by changing the contrast of some parameters. These variables now seem to have the right names. However the cohort and gender dummies do not. Regards, Oscar Edit: Fixed import issues. Remaining warnings can be ignored according to responses from other posts. Added second output link. I edited the education categories slightly. This seems to have effected the iteration process a lot, now I have a lot more warnings from clmm. I will look at the distribution of the education variables again to see if I can get more equal distribution across categories. link to output: https://www.dropbox.com/s/a2qrf5yvxpid9ct/Rquestions.docx link to output 2: https://www.dropbox.com/s/8tah94zr4dw0xs1/Rquestions.rtf


title: Boxplot with respect to two factors using ggplot2 in R

text: I am very new to R and to any packages in R. I looked at the ggplot2 documentation but could not find this. I want a box plot of variable `boxthis` with respect to two factors `f1` and `f2`. That is suppose both `f1` and `f2` are factor variables and each of them takes two values and `boxthis` is a continuous variable. I want to get 4 boxplots on a graph, each corresponding to one combination from the possible combinations that `f1` and `f2` can take. I think using the basic functionality in R, this can be done by               > boxplot(boxthis ~ f1 * f2 , data = datasetname)       Thanks in advance for any help.


title: Tobit model with R

text: Does anybody know where to find good application and examples (besides the manual and the book applied econometrics with R) using the tobit model with the packages AER? ### Edit I'm searching for a command to compute the marginal effects for y (not for the latent variable y*). It seems to be $\phi(x\beta/\sigma)\beta$, where $\phi$ is the std.normal cumulative distribution function. But how can I compute those effects with R?


title: Why are ties so difficult in nonparametric statistics?

text: My nonparametric text, Practical Nonparametric Statistics, often gives clean formulas for expectations, variances, test statistics, and the like, but includes the caveat that this only works if we ignore ties. When calculating the Mann-Whitney U Statistic, it is encouraged that you throw out tied pairs when comparing which is bigger. I get that ties don't really tell us much about which population is bigger (if that's what we're interested in) since neither group is bigger than the other, but it doesn't seem like that would matter when developing asymptotic distributions. Why then is it such a quandary dealing with ties in some nonparametric procedures? Is there a way of extracting any useful information from ties, rather than simply throwing them away? EDIT: In regards to @whuber's comment, I checked my sources again, and some procedures use an average of ranks instead of dropping the tied values completely. While this seems more sensible in reference to retaining information, it also seems to me that it lacks rigor. The spirit of the question still stands, however.


title: Analyzing discrete choice panel data with mlogit in R

text: I searched around and saw some high level discussion on mlogit and discrete choice panel models (here and here) but I need a more concrete answer than those. I am hoping somebody can point out what I am doing wrong with mlogit and point me in the right direction. I have a dataset where I ask participants to respond to 8 choice sets (CS) in random order. ResID is the participant identifier and CS is the choice set number. Each choice set is built of three products identified over three attributes (Price, Source and ext) at two levels (high low). The variables Compatibility and PriceBtw are individual level variables. I am trying to fit a mixed multinomial logit model to this data, with both product and individual level attributes. I got the data in a format supported by mlogit by               OSSdata_V0<-mlogit.data(data1, shape="long", choice="Choice", alt.levels=c(1,2,3), id="Strata", varying=4:6)      then I tried to fit the fixed effects part of the model by:               ml.OSS<-mlogit(Choice~Price+Source+Ext, data=OSSdata_V0)     Error in if (abs(x - oldx) < ftol) { :     missing value where TRUE/FALSE needed      I went through the data. I see no missingness anywhere. I was also able to run the fixed effects model in clogit and mclogit, the problem only occurs in mlogit. When I try to integrate the random effects part I get a different error message.               ml.OSS<-mlogit(Choice~Price+Source+Ext | PriceBtw, data=OSSdata_V0)     Error in solve.default(H, g[!fixed]) :     Lapack routine dgesv: system is exactly singular: U[10,10] = 0      My data is structured as follows and is in a data.frame called data1:                   Strata    ResID    CS    Price    Source    Ext    Product    Choice    Compatibility    PriceBtw time     1     1|1    1    1    1    0    1    A    1    1    $150    212     2     1|1    1    1    1    1    1    B    0    1    $150    212     3     1|1    1    1    0    1    0    C    0    1    $150    212     4     1|2    1    2    1    1    1    B    0    1    $150    212     5     1|2    1    2    0    1    0    C    0    1    $150    212     6     1|2    1    2    1    0    0    D    1    1    $150    212     7     1|3    1    3    1    1    0    E    1    1    $150    212     8     1|3    1    3    0    1    1    F    0    1    $150    212     9     1|3    1    3    0    0    1    G    0    1    $150    212     10    1|4    1    4    1    1    1    B    0    1    $150    212     11    1|4    1    4    1    0    0    D    1    1    $150    212     12    1|4    1    4    0    0    1    G    0    1    $150    212     13    1|5    1    5    0    0    0    H    0    1    $150    212     14    1|5    1    5    0    1    0    C    1    1    $150    212     15    1|5    1    5    1    0    1    A    0    1    $150    212     16    1|6    1    6    0    0    0    H    1    1    $150    212     17    1|6    1    6    1    0    1    A    0    1    $150    212     18    1|6    1    6    1    1    0    E    0    1    $150    212     19    1|7    1    7    0    1    1    F    0    1    $150    212     20    1|7    1    7    0    0    1    G    1    1    $150    212     21    1|7    1    7    1    0    0    D    0    1    $150    212     22    1|8    1    8    0    0    0    H    0    1    $150    212     23    1|8    1    8    0    1    1    F    1    1    $150    212     24    1|8    1    8    1    1    0    E    0    1    $150    212     25    2|1    2    1    1    0    1    A    1    2    $150    865     26    2|1    2    1    1    1    1    B    0    2    $150    865     27    2|1    2    1    0    1    0    C    0    2    $150    865     28    2|2    2    2    1    1    1    B    1    2    $150    865     29    2|2    2    2    0    1    0    C    0    2    $150    865     30    2|2    2    2    1    0    0    D    0    2    $150    865


title: Chi Square for this case?

text: I am not sure how to perform the statistical analysis on the following table. I did an experiment in which 12 participants had to choose between 3 conditions when provided with 3 stimuli.                 Stimulus  Condition1  Condition2 Condition 3         A            9          1          2         B           10          2          0         C            8          2          2      I want to prove that it is not by chance that Condition 1 is preferred rather than the other two conditions. How can I do this analysis? Maybe with a Chi Square test? If yes should I group Condition 2 and 3 against Condition 1? I use R, in case could you please provide an R example in order to analyze these data? Thanks in advance


title: Help with covariance

text: > A deck of 52 cards is shuffled and a bridge hand of 13 cards is dealt out. > Let X and Y denote, respectively, the number of aces and the number of > spades in the hand. > > (a) Show that X and Y are uncorrelated. > > (b) Are they independent? I know that if ${\rm Cov}(X,Y)=0$ then $X$ and $Y$ are uncorrelated. I also know that $X$ and $Y$ are independent since the probability of choosing aces does not affect the probability of drawing spades. I tried using ${\rm Cov}(X,Y)=E[XY]-E[X]E[Y]$ but got stuck on trying to find the expected value of $X$ and $Y$. Is this a valid approach to the problem or is there something I am missing?


title: Smart application of transformations in digit recognition

text: So I'm working on some machine learning code to recognize the digits 0-9 using singular value decomposition, followed by a least squares comparison with the first 15 or so singular vectors as my basis. Just with that alone the accuracy is quite good, but I'm hoping to improve it by incorporating direct 2-norm digit comparison against the entirety of my training examples when there is a close call between the 1st and 2nd min in my least squares SVD algorithm. So for example if I have a test digit which is classified based on my SVD algorithm as a $2$ but where it was a really close call whether it was maybe a $3$, then I would compute the first order approximation of the test digit's x,y-shifts, x,y,diagonal-stretches and rotations, and compare it to these same transformations precomputed for all my training digits. What I'm trying to decide is if whether or not this is really a smart way to use transformations to improve my accuracy. The amount of comparisons (vector subtraction followed by 2-norm) needing to be done when one of these close calls occurs is nearly a million (I have 60,000 training digits), which seems prohibitively expensive. Furthermore, its not clear to me that this type of transformation comparison algorithm is even a good tie breaker when my SVD algorithm fails to pick a clear winner. Moreover, since all I'm using is the first order approximation to these transformations, large transformations cause significant divergence from the non-linear transformation manifold, which makes me wonder whether just direct comparison to the training digits alone, not bothering at all with comparing transformations, would yield nearly as good results. Does anyone have any thoughts on how I might better incorporate transformations into my digit recognition code? Do I need a second order approximation? Should I compare only to some random subset of my training digits to cut down on computation costs? etc.


title: The differences between models via their resampling distributions.

text: The `caret` package offers the ability to make statistical statements about the performance of different models used for classification. According to the description, it uses asymptotic simultaneous conﬁdence intervals for many-to- one comparisons of proportions. I have read the papers by Horthon et al (2005) and by Eugster et al. 2008 but yet still I cannot summarize what the results mean nor how it was arrived at. I am trying to interpret this graph (taken from the `caret` page. ![enter image description here](http://i.stack.imgur.com/UvFrO.png) Can anyone offer a brief description on an interpretation of the graph, especially the x-axis scale (-0.4 to 0.6). Any comments will be welcomed. Kurt


title: What is Precision-Recall Curve?

text: I have a data mining assignment where I make a content-based image retrieval system. I have 20 images of 5 animals. So in total 100 images. My system returns the 10 most relevant images to an input image. Now I need to evaluate the performance of my system with a Precision-Recall curve. However, I do not understand the concept of a Precision-Recall curve. Let's say my system returns 10 images for a gorilla image, but only 4 of them are gorillas. The other 6 images returned are other animals'. Now my precision is `4/10 = 0.4` ( relevants returned )/ (all returned) and my recall is `4/20 = 0.2` (relevants returned) / (all relevants) So I only have a point `<0.2,0.4>` not a curve. How do I have a curve, i.e. a set of points? By changing the number of images returned (which is fixed = 10 in my case) ? Any help would be appreciated, Thanks !


title: Statistics for conditions with different number of levels

text: I have measured the brain activity in a group of 10 subjects who had to perform a task in `real` and `imagined` conditions, the task is about stopping a hand movement when a stop signal appears. In the `real` condition I have 2 sub-levels (succesfully stopped `SS` and unsuccesfully stopped `US`), and I don't have these sub-levels in the imagined condition, there aren't any overt response to classify as SS or US in imagination. My main objective is to compare the brain activity (using event related potentials) between the `real` and the `imagined`conditions, a simple t-test can do it. But this activity is modulated if the task is `SS` or `US`, and taking both only as the `real` condition dilutes some effects. Is there any way to compare between: `imagined-real`, `imagined-SS`, `imagined-US` and `SS-US`? I think that a one-way ANOVA with four levels is not a good way to compare these data...


Are samples from sliding windows independent samples? E.g. if I have window size of 90 seconds counting the number of cars on a street and I output the average within the window every second for 30 seconds, do I have 30 independent samples or not? I'd say yes, as it looks (to me) like sampling with (partial?) replacement. But I'm not sure. I'm asking, because I thought if the samples were independent and I collected n>30, the central limit theorem could be used for further calculations.


title: How to calculate 'weighted' pooled/composite statistics for two given groups given group statistics

text: I have two sets of independent samples from the same distribution. For each, I have calculated   1. sample weighted mean (u1, u2)   2. sample weighted std deviation (d1, d2)   3. sample weighted std error (e1, e2) The question below gives a way to calculate the pooled statistics in the 'unweighted' case How to calculate pooled variance of two groups given known group variances, means, and sample sizes? So, how can I calculate this for the weighted case? Edit 1: To clarify, EACH sample has a 'weight' associated with it which represents a 'confidence' associated with it 0 < weight < 1 Edit 2: To clarify further, each sample represents return on money invested in a particular investment. The weight of a sample is the fraction of the total money invested in a given interval. The two groups of samples are from two different time intervals. The total money invested in the two intervals is not the same


title: Pre-post treatment design: accounting for reduced effect of treatment in baseline high scorers

text: I'm planning a study in which I want to test the effect of a treatment on a dependent psychometric variable. I expect subjects who score lower at baseline to benefit more from the treatment (larger score difference pre-post), and I expect the treatment to have less of an effect on subjects who score high at baseline (smaller score difference pre-post). This is not due to a shortcoming of the psychometric measure but due to the underlying psychometric construct. I have empirical data of a small ( _n_ = 10) pre-study with identical design that I could use to quantify this effect. Is there a more elegant solution to test this than to perform a median split?


title: ANOVA normality assumption for which variables?

text: When I want to conduct a repeated measures ANOVA pre and post test scores, which variables need to be tested for normality? Each the pre and post scores or the difference between them? I suppose the difference as with the t-test. Correct?


title: ARIMAX for modelling daily sales

text: I am trying to model daily sales for a take out restaurant. They are only open on business days - no holidays or weekends - as their primary clients are office workers on their lunch breaks. Below is what two years of the daily sales time series looks like. ![enter image description here](http://i.stack.imgur.com/UhHDR.png) The days with zero sales, as you can see above, are the days that the restaurant was closed due to a public holiday (Easter Monday etc.). There is definitely a weekly pattern: sales tend to peak on Thursdays. Furthermore, the presence of a holiday changes the sales pattern in the surrounding weeks (the week before and the week after). You might notice that there are sales spikes before or after certain holidays. An example of this: if Monday is a holiday, sales tend to be much lower on the Friday before that long weekend - presumably office workers leaving work early. There are also yearly seasonal patterns. Sales are lower in the summer, for example, presumably, in part, because many office workers are taking their vacations. My approach has been to use a ARIMAX model to fit the data (using R). I've followed the approach suggested by Rob Hyndman here. The difference is that I'm using only Mon-Fri, so my frequency is 5, and I've added dummy variables for all of the days that the restaurant is closed (holidays). The model fit is not very good so far, of course. I haven't done anything to take into account the effect of a holiday on the surrounding days. Further, I'm including the holidays as days with sales equal to zero, so this must throw off the model. Here is what R returns:               ARIMA(0,1,1)(1,0,1)[5]                              Coefficients:          ma1     sar1    sma1      Mon       Tue       Wed      Thu     Day    Newyears   FamilyDay  GoodFriday      Easter  VictoriaDay   CanadaDay    CivicDay       -0.804  -0.2608  0.3255  54.4530  113.8052  152.0052  -6.3025  0.0388  -1545.1973  -1604.5038  -1581.6740  -1586.8710    -1628.253  -1437.6075  -1181.0054     s.e.   0.028   0.4641  0.4529  23.2788   23.3748   23.4900  23.4367  1.5546    117.6128    113.6446    113.8825    114.5609      114.786    112.7561    114.0031        LabourDay  Thanksgiving       -1310.3416    -1332.8028     s.e.    113.5081      113.5179          sigma^2 estimated as 28269:  log likelihood=-3305.08     AIC=6646.16   AICc=6647.56   BIC=6722.2      I think I should include the month of the year as a dummy variable to capture yearly seasonal effects as well. My questions:   1. What can I do to capture "long weekend effects"? Should I included a dummy variable for every Friday that proceeds a long weekend etc?   2. How should I deal with the holidays that the restaurant was closed for? If I remove them, then the week lengths will not be the same. If I include them, then they are outliers that throw everything off.   3. What else can I do to improve my model? Thanks very much for any input.


title: Why is Shapiro–Wilk test considered to be the best normality test?

text: I have read somewhere in the literature that the Shapiro–Wilk test is considered to be the best normality test because for a given significance level, $\alpha$, the probability of rejecting the null hypothesis if it's false is higher than in the case of the other normality tests. Could you please explain to me, using mathematical arguments if possible, how exactly it works compared to some of the other normality tests (say the Anderson–Darling test)?


title: How to solve for the parameters in a logistic function?

text: I want to find the parameters of a logistic function. I read the guide here. It has a very clear explanation, but it did not have the final solution that I need. Now, we will consider a basis logistic function: $$h_\theta(x^{i})=\frac{1}{1+e^{-\theta_0-\theta_1x}}$$ We want to find $\theta_0$ and $\theta_1$ subject to minimum cost function: $$J(\theta)=-\frac{1}{m}\sum_{i=1}^{m}y^{i}\log(h_\theta(x^{i}))+(1-y^{i})\log(1-h_\theta(x^{i}))$$ We use the notation: $$\theta x^i:=\theta_0+\theta_1 x^i_1 $$ Then $$\log h_\theta(x^i)=\log\frac{1}{1+e^{-\theta x^i} }=-\log ( 1+e^{-\theta x^i} ),$$ $$\log(1- h_\theta(x^i))=\log(1-\frac{1}{1+e^{-\theta x^i} })=\log (e^{-\theta x^i} )-\log ( 1+e^{-\theta x^i} )=-\theta x^i-\log ( 1+e^{-\theta x^i} ),$$ and $$J(\theta)=-\frac{1}{m}\sum_{i=1}^m \left[y_i\theta x^i-\theta x^i-\log(1+e^{-\theta x^i})\right]=-\frac{1}{m}\sum_{i=1}^m \left[y_i\theta x^i-\log(1+e^{\theta x^i})\right],~~(*)$$ where the second equality follows from: $$-\theta x^i-\log(1+e^{-\theta x^i})= -\left[ \log e^{\theta x^i}+ \log(1+e^{-\theta x^i} ) \right]=-\log(1+e^{\theta x^i}). $$ All I need now is to compute the partial derivatives of $(*)$ w.r.t. $\theta_j$. As $$\frac{\partial}{\partial \theta_J}y_i\theta x^i=y_ix^i_j$$ $$\frac{\partial}{\partial \theta_j}\log(1+e^{\theta x^i})=\frac{x^i_je^{\theta x^i}}{1+e^{\theta x^i}}=x^i_jh_\theta(x^i), $$ The above steps are correct, but it did not have the solution for $\theta_0$ and $\theta_1$. We can take the derivative and put it equal $0$. Namely, $$\frac{\partial J(\theta)}{\partial \theta_0}=y-\frac{e^{\theta_0+\theta_1 x}}{1+e^{\theta_0+\theta_1 x}}=0~~(*)$$ $$\frac{\partial J(\theta)}{\partial \theta_1}=yx-\frac{xe^{\theta_0+\theta_1 x}}{1+e^{\theta_0+\theta_1 x}}=0 ~~(**) $$   * Are $(*)$ and $(**)$ correct?    * Assume that $y$ and $x$ are known. How do you represent $\theta_0$ and $\theta_1$ by $y$ and $x$?


> **Possible Duplicate:**   >  When is logistic regression solved in closed form? Why is there no analytic solution for logistic regression? I was trying to derive a solution similar to the normal equations of linear regression, and I am not sure where it goes wrong. Please provide a mathematical proof.


title: Column Means Significance Tests in R

text: New to the site. I am just getting started with R, and want to replicate a feature that is available in SPSS. Simply, I build a "Custom Table" in SPSS with a single categorical variable in the column and many continous/scale variables in the rows (no interactions, just stacked on top of each other). The table reports the means and valid N's for each column (summary statistics are in the rows), and select the option to generate significance tests for column means (each column against the others) using alpha .05 and adjust for unequal variances. Here is my question. How can I replicate this in R? What is my best option to build this table and what tests are available that will get me to the same spot? Since I am getting used to R, I am still trying to navigate around what is available. Many thanks in advance!


title: How can I get R to stop autocompleting my quotes/parens?

text: When I type a left paren or any quote in the R console, it automatically creates a matching one to the right of my cursor. I guess the idea is that I can just type the expression I want inside without having to worry about matching, but I find it annoying, and would rather just type it myself. How can I disable this feature? I am using R 2.8.0 on OSX 10.5.8.


title: Application of wavelets to time-series-based anomaly detection algorithms

text: I've been beginning to work my way through Statistical Data Mining Tutorials by Andrew Moore (highly recommended for anyone else first venturing into this field). I started by reading this extremely interesting PDF entitled "Introductory overview of time-series-based anomaly detection algorithms" in which Moore traces through many of the techniques used in the creation of an algorithm to detect disease outbreaks. Halfway through the slides, on page 27, he lists a number of other "state of the art methods" used to detect outbreaks. The first one listed is **wavelets**. Wikipeida describes a wavelet as > a wave-like oscillation with an amplitude that starts out at zero, > increases, and then decreases back to zero. It can typically be visualized > as a "brief oscillation" but does not describe their application to statistics and my Google searches yield highly academic papers that assume a knowledge of how wavelets relate to statistics or full books on the subject. I would like a basic understanding of how wavelets are applied to time-series anomaly detection, much in the way Moore illustrates the other techniques in his tutorial. Can someone provide an explanation of how detection methods using wavelets work or a link to an understandable article on the matter?


title: Interaction in generalized linear model

text: I have 2 continuous variables as my predictors and the interaction between them, so 3 effects all together (when I center my predictors only the interaction is significant). I am using a binary logit model except where I have fixed the value of the number of trials at 20, for my dependent variable. My problem has to do with further understanding the nature of the interaction. I know that in generalized linear models the interaction is more complex compared to the linear model because of the link function. I have read that interpreting the sign of the parameter estimates is very limited, so I am trying to find a way to further understand the interaction, and graphing it would be nice as well. Unfortunately, I only have SPSS at my disposal. Jeff Gill (http://www.artsci.wustl.edu/~jgill/papers/interactions3.pdf) wrote about a method called first differences which seems very similar to the inteff command in stata. If anybody could tell me how to how to further understand and test the interaction in a generalized linear model using spss or hand calculations I would be eternally grateful.


title: Question on re-indexing data for graphing

text: A file contains data for 1975-2012 of housing prices across a panel of countries. However, it's indexed to 2005, so each country's housing price series converges to 100 in 2005, then they spread apart again to 2012. This makes graphing all data series on one chart for comparison's sake unintelligible. Without the original data, could one properly reindex the data from 1975 in a statistically correct manner, and how could someone do this in either R or excel?


title: Alternative to multiple Fisher's Exact tests

text: I have count data for a pair of test and control populations of samples for a variable that takes a value of either A or B for a sample:                           A   |   B              -------+-------     Control     7   |   1     Test        3   |   5      The hypothesis is that the test condition stimulates the conversion of the observed parameter from a state with a value of A to a state with a value of B, and I can use the Fisher Exact test on the data to accept or reject this hypothesis. Now suppose we have three such pairs of test and control populations, with varying numbers of observations carried out for each population for the test or control condition.                               1                2                 3                      A   |   B        A   |   B         A   |   B              -------+------   -------+-------   -------+-------     Control     7   |   1       70   |  10         7   |   1     Test        3   |   5        3   |   5       300   | 500      We can conduct the Exact test thrice, once for each pair, and obtain three P values, which we somehow assess together (majority rule?) to come to a conclusion. But, is there some other statistical test that can be used, perhaps one that yields a single P value, or is more robust? [Adding information about the actual experiment] RNA molecules produced from a certain gene in the cells of the human body get edited because of exposure to the test condition, resulting in the conversion of a specific nucleotide base of the RNA from being an adenine (A; equivalent to the hypothetical state A above) to a guanine (G; equivalent to the hypothetical state B above). The editing is assayed by sequencing the RNA molecules, with varying number of molecules being sequenced for the control and test conditions. The three tables allude to data from three humans.


title: Logistic regression on 3 covariates each with unequal sample size

text: I have 3 covariates $(x_1,x_2,x_3)$, each covariate takes only 2 values $\\{0,1\\}$ but each covariates have unequal sample size. $n(x_1)=79,\;n(x_2)=80,\; n(x_3)=77$ $x_1$:(60 zeros and 19 ones) $x_2$: (56 zeros and 24 ones) $x_3$: (34 zeros and 43 ones) My goal is to estimates: $$p = \frac{1}{1+e^{-(b_0+b_1x_1+b_2x_2+b_3x_3)}}$$ and fit a logistic regression model using all these three $x_i$'s , how can I do this?


title: Compare two time series

text: I have two time series that come from the same system. One is taken from the whole system, and the other is taken from, say, 10% of the system. The two time series have the same frequency. Is there a measure to know how similar these two time series are?


title: VAR model with time series of different frequencies

text: In case I want to see the effect of two or more endogenous time series on each other, I use a VAR model. But how do I proceed if one data set is monthly, and the other one daily?


title: Choice of sampling method before using bootstrap validation method for computer models

text: Bootstrap method sounds to usually work for observational data. In my research, I want to use the bootstrap validation method for constructing a surrogate model for a time-consuming engineering model. The bootstrap validation method is to use the Efron-Gong optimism bootstrap to estimate the prediction error of surrogate models for model selection. Then firstly I need to create an original data set, i.e. the combinations of several inputs and outputs from this engineering models by using a sampling method. My question is whether it matters for choosing a sampling method, including Latin hypercube sampling (LHS), quasi-random method (Sobol or Helton), or simple random method. If I use LHS to construct an original data set, then I will bootstrap this LHS sample with replacement to create a new bootstrap sample, which may lose the advantages of using lhs method. This is because this new bootstrap sample is non-lhs sample. In this case, it would be better if I use quasi-random sampling in the first place? Any thoughts on this would be greatly appreciated.


title: Best predictive Cox regression model using c-index and cross-validation

text: I have a set of 40 genes expression data as binary variables (overexpressed yes or not) for predict recurrence in a kind of cancer in a sample of 77 patients, where 22 has recurred. This set is called discovery set. The objective is to find a subset of these genes to predict the recurrence in cancer. Times to recurrence and follow-up times for censored observations were saved. In a first step, I have fitted univariate Cox regression model, and I have discarded 13 genes of the original set of 40 genes, removing those that have P-values > 0.25 in the LRT. Then, I work only with 27 genes. In the second step, I am analyzing all the possible multivariate models of 3 variables, 4 variables, 5 variables, etc. I decide what is the more predictive model of $X$ variables fitting multivariate Cox regression and calculating the average c-index of 100 times 5-folder CV. For calculate the c-index I use the “index.corrected” reported by the `validate` function of the Harrell’s `rms` package, and I calculate the average of the 100 times I run the CV. I will use the AIC for compare the best model with $X+1$ variables and the best model with $X$ variable, to decide what it is the best predictive model. Finally, I have an independent validation set of patients ($n=105$) where I will check the model. I will use the linear predictor of the best model in the discovery set as prognostic index, dividing in tertiles of risk, and calculating the values of the prognostic index of these new patients. I will analyze statistical differences of the 3 groups with the log-rank test, HR with 95CI% and KM Plots.


title: 70% certain with a 70% success rate

text: We have some classifiers that classify text documents. Each classifier reports how likely their classification is correct based on what they know of the document. Classifiers do not know how accurate they themselves are. If a classifier reports that it is 70% likely that document belongs to class A and that classifier is correct in its classifications 70% of the time, what is the likelihood that the document actually is a member of class A. Surely not 70% of 70%, so what would the actual calculation be?


title: Dealing with ties, weights and voting in kNN

text: I am programming a kNN algorithm and would like to know the following: Tie-breaks:   1. What happens if there is no clear winner in the majority voting? E.g. all k nearest neighbors are from different classes, or for k=4 there are 2 neighbors from class A and 2 neighbors from class B?   2. What happens if it is not possible to determine exactly k nearest neighbors because there are more neighbors which have the same distance? E.g. for the list of distances `(x1;2), (x2;3.5), (x3;4.8), (x4;4.8), (x5;4.8), (x6;9.2)` it would not be possible to determine the k=3 or k=4 nearest neighbors, because the 3rd to 5th neighbors all have same distance. Weights:   1. I read it is good to weight the k-nearest neighbors before selecting the winning class. How does that work? I.e. how are the neighbors weighted and how is then the class determined? Majority vote alternatives:   1. Are there other rules/strategies to determine the winning class other than majority vote?


title: Is there a method to know if we can apply regression to certain dataset and obtain satisfactory results?

text: I am working on 84 images, where each image is assigned a target value Y (could be any measure such as blur or noise level). Y is already known for all images. Features X (predictors) are extracted from each image. There are many methods for feature extraction, and most of them return large number of features. The problem turns to be a regression problem, where we need to predict Y given X (the features). Having many methods to construct the features, I need a way to assess the ability of certain features to predict Y. I am thinking of using Support Victor Regression with Radial Basis kernel for prediction, I will also investigate other methods, such as decision trees and linear regression. **My question is:** given a dataset, can we analyze the data and know if we can apply regression analysis on it and obtain good results (e.g. > 70% accuracy)? (Either linear or non-linear regression) Please keep in mind that there are too many features, so analyzing the correlation between X and Y, and the multi-colinearity among X is infeasible, unless there is an automated way.


title: Repeated loess smoothing for time series data

text: I have 3 time series and I would like to determine if there are any common patterns in the data. Zuur mentioned in one of his books that a reasonably simple approach to test this is based on repeated loess smoothing. Zuur states that to visualize the general patterns of the time series, we fit a loess smoother through each series and then plot them all on the same graph. This generally captures the long term variations, to capture the short term variations (which I am interested in) we have to apply loess smoothing again, but now on the residuals. Below, is an example data set where we have hourly values of ozone concentrations for 3 countries, I would like to perform what Zuur has explained, but to this example data set. Could anyone offer some suggestions? Please note that my knowledge of R and statistics in general is low, I apologize if this is a simplistic question.               require(plyr)     require(lattice)          TopFolder <- list("http://www.nilu.no/projects/ccc/onlinedata/ozone/CZ03_2009.dat"                       ,"http://www.nilu.no/projects/ccc/onlinedata/ozone/CY02_2009.dat"                       ,"http://www.nilu.no/projects/ccc/onlinedata/ozone/BE35_2009.dat"     )     ## create variable for the data      data = ldply(TopFolder, header = TRUE, read.table, sep = "", skip = 3)     ## define Ozone levels     Ozone <- data$Value         Ozone[Ozone==-999] <- NA         Ozone <- data.frame(Ozone)         ## define Datetime -  need to concatenate arrays         DateTime <- paste(data$Date,data$Hour, sep = " ")         Date <- as.POSIXct(DateTime, format = "%d.%m.%Y %H:%M")         ## define countries          Countries <- c("Country1","Country2","Country3")         Country <- data.frame(Country = rep(Countries, each = 8760))         ## bind together         Dat <- cbind(Ozone, Country = Country)         Dat <- transform(Dat, Date=Date, Doy = as.numeric(format(Date,format = "%j")),                          Tod = as.numeric(format(Date,format = "%H")),                          DecTime = rep(seq(1,365, length = 8760),by = 3))         ## Select data for a month during the summer          NewDat <- subset(Dat, as.Date(Dat$Date) >= '2009-07-01 00:00:00' & as.Date(Dat$Date) <= '2009-08-01 00:00:00')           ## plot the data     xyplot(Ozone~DecTime | Country, data = NewDat, type  = "l", col = 1,            strip = function(bg = 'white',...)strip.default(bg = 'white',...))      So, here I have imported the data and taken a subset (for one month), and then plotted the values. I would now like to apply loess smoothing to represent the short term variation i.e. hourly.


title: What determines to which variable the explained variance is granted?

text: What happens when there is a linear regression of the kind `y ~ x1 + x2 + x3 + x4` What determines which variables will be granted the explained variance rather than the others given a correlation matrix (not assuming independence)? In case this really is unpredictable, why do we insist on using regression rather than simply looking at the correlation matrix: at least then you know you are not accounting for effects of other variables when looking at the correlation coefficient, while using regression you'd be oblivious to which variable gets the explained variance (and how correctly or wrong this is!).


title: Problems optimizing a likelihood function in R

text: I'm trying to maximize the log-likelihood $y_1 \log(p_1) + (1-y_1)\log(p_1) + y_2\log(p_2) + (1-y_2)\log(p_2)$. I have data on success of an experiment. Y identifies if person died or not, and X identifies control or treatment group. What's wrong with my code?               # MLE for the likelihood     y <- c(rep(1,39), rep(0,674-39), rep(1,22), rep(0,680-22))     x <- c(rep(0, 674), rep(1, 680))          binom.ll <- function(theta, y, x) {       y[x==0]*log(theta[1]) + (1-y[x==0])*log(1-theta[1])) + y[x==1]*log(theta[2]) + (1-y[x==1])*log(1-theta[2]))     }          theta.start <- c(0, 0)     ml.res <- nlm(binom.ll, theta.start, print.level=1, y=y, x=x, hessian=T)     ml.res          iteration = 0     Step:     [1] 0 0     Parameter:     [1] 0 0     Function Value     [1] 1.797693e+308     Gradient:     [1] -Inf    0          Error in nlm(binom.ll, theta.start, print.level = 1, y = y, x = x, hessian = T): non-finite value supplied by 'nlm'     In addition: Warning messages:     1: In nlm(binom.ll, theta.start, print.level = 1, y = y, x = x, hessian = T): NA/Inf replaced by maximum positive value     2: In nlm(binom.ll, theta.start, print.level = 1, y = y, x = x, hessian = T): NA/Inf replaced by maximum positive value


title: Treating non-stationarity of time series in seasonal adjusted data with R

text: I'm currently trying to use a variable x (and others) to explain a dependent variable y in a distributed lag model (with the long term goal of predicting variable y). The plot of variable x shows an evident seasonality at the end of the year: See http://i.imgur.com/8gGtaVS.png After having deseasonalized the data with decompose (with multiplicative components), adf and kpss tests indicate, that the seasonal adjusted data is still not stationary. Because there are more independent variables and I don't want to look deeper into investigation a cointegration relationship between those series, I thought it would be the most usual way to take the difference of both series (with the diff() function). Now there are 2 alternatives:   1. Take the diff of the already seasonal adjusted data. The problem with this approach is, that I'm not sure if this a good idea because I don't see how you can reseasonalize the time series for a forecast result in an easy manner.   2. Take the diff of the raw series. This leads to the following graph. See http://i.imgur.com/loDU2IE.png The time series is now stationary regarding adf and kpss tests, however there is still the seasonal pattern visible. Now I'm not sure if it is recommended to use decompose (with the multiplicative) method, to deseasonalize the diff of the time series, especially because there are zero values which have no effect when calculating the seasonal adjusted time series in the following way:               decompose(x_diff, "mult")$x / decompose(x_diff, "mult")$season      So, how should I proceed when I want to include (the diff of) x as a independent (and lagged) variable in a distributed lag model?


title: Is there any theory governing factors vs flags

text: I'm doing some work in `R` using the `gbm` package. I'm curious about the repercussions of treating categorical variables as factors or using `1/0 flags` for each individually. Is there any literature on how the 2 are handled either in `gbm` or in `R` in general?


title: R regression by groups and data visualization in one graph~

text: Hi currently I am conducting simple linear regression on two variables for data of different regions. I know that I can use the lmList function to get the coefficients at once for all the regions. But can I generate the Q-Q residual plot for all the regions in one graph with different panels at once? And for the output for lmList function, only the coefficients are displays, without R-square for each regression. How can I see that? Thanks a lot for help in advance!


title: Metrics for cluster evaluation

text: I make a set of clusters using some clustering algorithm. Precision, Recall, F Measure, Fallout and RI of individual clusters are calculated for testing the performance. How do I calculate the average precision, recall, f measure, etc.? Is it the average of the different clusters' precisions? How do I get a single precision number for the _N_ total clusters?


title: Selecting KDE bandwidth by n-fold cross-validation

text: I know that various questions already deal with this topic, but I don't think anyone has answered this before! I'm reading a paper in which the authors claim to use 20-fold cross validation to estimate the bandwidth of a Gaussian KDE. What does this mean? I've read up on cross validation for KDE bandwidth estimation, but every example I have found (including the 'standards' least squares and maximum likelihood) don't have an 'n-fold' mentioned, because they use the leave-one-out approach (hence n is always equal to the number of data points). Thanks for any light you can shed.


title: Performing Mutiple Correspondence Analysis (MCA) using R?

text: I want to apply MCA (multiple correspondence analysis) to a data table containing regions as individuals, and INFESTED AREA as variable and degree infested also and the type of trees infested , and I want to apply MCA to it using XLSTAT or any program, but I have problem in organizing this data, here a sample of it: A sample for the year 2013 ![enter image description here](http://i.stack.imgur.com/wOwdq.png) And other for 2012: ![enter image description here](http://i.stack.imgur.com/mV7BP.png) I want to find a way to organize these data, make them as categories to better fit for MCA analysis? And I want to know if I can apply MCA for different years on the same table?


title: Modelling the difference between two Hubbert curves

text: I have a large set of time series data, consisting of series from two different conditions, the averages of which are shown below. I would like to fit a model to this data, to test that a) the peak value is greater in the 'conflict' condition, and b) the peak occurs earlier in this condition.               ggplot(data, aes(x=Time, y=Variable, colour=Condition))          + stat_summary(fun.data=mean_se, geom="pointrange")      ![enter image description here](http://i.stack.imgur.com/3NQ8T.png) From my own research on this, I know that:   * Growth curve analysis is the usual way of modelling time series data like this, but I can't figure out how I would fit a polynomial for this shape of curve.      * I have some experience fitting GCA models using `lme4` in R, mostly following Dan Mirman's tutorials, but I'm still learning.   * Curves of this kind are referred to as Hubbert curves, and typically used to model oil production, as a symmetric logistic curve up and down.   * R package `grofit` is maybe useful for analyses of this kind, although I would rather know I'm barking up the right tree before investing time in learning how to use this. Can anyone point me in the right direction here?


title: getting from Edgeworth expansions to Cornish Fisher Expansions

text: I am currently trying to work out how to get from the Edgeworth expansion to the Cornish-Fisher expansion. I use van-der-Vaarts "Asymptotics Statistics" and Hall's book on Edgeworth expansions and the bootstrap. Unfortunately neither help much to get the details right (vdVaart p 338). tldr: Edgeworth expansions can be inverted to get Cornish Fisher expansions - but how? The paragraph im having trouble with is the following Let $\Phi, \phi$ be the normal cdf, pdf, $p_1$ a polynomial, $z_\alpha, \hat{\xi}_{n,\alpha}$ the normal and bootstrap upper $\alpha$ Quantiles. The book states that $$ 1-\alpha = \Phi(\hat{\xi}_{n,\alpha})+\frac{p_1(\hat{\xi}_{n,\alpha}\mid\hat{P}_n)\phi(\hat{\xi}_{n,\alpha})}{\sqrt n} + O_P(n^{-1}) $$ can be somehow inverted to obtain $$ \hat{\xi}_{n,\alpha}=z_\alpha-\frac{p_1(z\alpha \mid P)}{\sqrt n}+ O_P(n^{-1}) $$ apparently by Taylor expanding $\Phi, \phi$ and $p_1$ around $z_\alpha$ I guess it all makes sense and I can accept that it is possible, but I'd like to see formal justification.


title: 10 minutes wind speed data to 1 second wind speed data

text: I have a wind speed data series .txt file (1 year long), in which in each register I have the following info: date; hour; 10 minute wind speed average; 10 minute max value; 10 minute sigma An example of the data would be like this: 050206 0130 8.05 10.28 0.84 050206 0150 7.29 11.06 1.13 .... For each 10 minutes, the available information is the wind speed average, the maximum value in such 10 minutes and the standard deviation in such 10 minutes. Such statistics for each 10 minutes are calculated with 1 s sampling data (600 values for each 10 minutes). It can be assumed that the 1 s sampled values are Normally distributed. The problem is that such 1 s sampled values are not available and I need them for my study. So, the questions are: -Is there a mathematical/statistical solid well known method to generate the 1 s values for each 10 minute register assuming that they are Normally distributed?? -Is there a mathematical/statistical solid well known method to link each 1 s values with coherence?? I would like to implement a Python script in order to do so.


title: ANOVA Best test for correlation with large data set and many independent variables?

text: I'm trying to analyze some factors contributing to win rates of a game, there are several hundred factors but each game will only have a small subset of them (10-20). Some of the factors may be correlated (picking one ability will tend to blend well with other abilities, and some may be antagonistic and so will not often be picked together). Ultimately I want to know the correlation these factors have on winning, which is obviously a binary value. I have a database of several hundred thousand games so scale is not an issue, but I'm stuck at what sort of test to use. I've used ANOVA tests in university for data before so my thinking is maybe that would work, but I've never had so many variables to fit, and I'm not sure my data is normally distributed, I just wanted to check to make sure that a one-way ANOVA would make sense in this situation, or if I should be looking for another way to figure out the correlation on factors and winning. Thanks!


title: What exactly does 'representative sample' refer to?

text: When reading passages like the following: > > Based on a representative sample of 88 recent raids, we show that the > Turkana sustain costly cooperation in combat at a remarkably large scale, at > least in part, through punishment of free-riders. I wonder what 'representative sample' might refer to. Is it related to power calculations (for example) in statistical inference or is there some way to assess the number of samples required from the total population for it to be considered representative?


title: Seeding problem for a 3-way bracket

text: I apologize in advance if this is not the correct place for this question. I'm going to be running a 36-person tournament in a certain activity, with the first two rounds being 3-way competitions (36->12->4). My question is how to fairly seed these 3-way tournaments. Before the round of 36, I'll have a good idea of player strength, and the seeds will reflect that strength. I can either decide matchups by seed or by a rough estimation of player strength. My idea was to give the highest strength of schedule to the 36th seed, the 2nd-highest to the 35th seed, etc. Is there any good method out there for figuring out matchups so that this can happen, or is there another method which would be appropriate for seeding such matchups?


title: Log squared error accuracy justification: data mining competition

text: I was wondering what is the justification/rational for using the following error measure for judging accuracy instead of simply squared error? ![enter image description here](http://i.stack.imgur.com/LHMTG.png) Link to competition: Heritage Health


title: Variance of a marginal order distribution less than the variance of the full distribution?

text: Is there a proof/theorem that states that $$\mathrm{Var}(X_{kn}) \leq \mathrm{Var}(X)$$ where $X$ is a continuous random variable with distribution $F$, and $X_{kn}$ is the $k^{th}$ of $n$ order statistics from a sample of size $n$ from $F$?


title: Cross-Correlation for non stationary data

text: I am using cross correlation to find the correlation between two time series data say X and Y. I have read this somewhere that : " If X or Y contains auto-correlation or is not stationary about mean , the cross- correlation will not reflect the true relationship between X and Y. What to do in this case? Thanks in advance


title: Is the Kalman filter actually forecasting?

text: The state space equation is: $$Y_t = F_tθ_t + v_t\hspace{4em} \textrm{eq. 1}$$   $$θ_t = G_tθ_{t-1} + w_t\hspace{2.8em} \textrm{eq. 2}$$ $F_t$ in eq.1 are the independent variables and we can predict $Y_t$ provided we know $F_t$. Say that I am trying to forecast the number of hospitalizations due to flu based on humidity. In this case, based on eq.1 it looks like I need to know the humidity of next week to forecast the hospitalizations for next. In reality, I won't have any data about my independent variables for future too. I am actually not clear as to why this is called forecasting. Am I missing something here?


title: Is Time Series decomposition right approach?

text: I have sales data for last four year. I want to divide my data in two parts so that I can have base volume and promotional volume. I am thinking to use " Decomposition of Time Series " approach such that the base volume would be trend component and seasonal and error component would be promotional volume. Is this method would be good? is there any other method/way to break down the data in two parts. P.S : promotional volume is the volume that comes from promotions during festival seasons. Thanks in advance... Please reply......


title: Can a model for a resampled time series be useful for modeling the original time series?

text: Say we have a time series $x_i, i=1, \dots$. If it is hard to model, then is it possibly acceptable to model the resampled one $y_i := x_{ik}, i = 1, \dots$ for some fixed natural number $k$? (Assume no special requirement or restriction, just a general case.) If I have a model for $y_i$'s, can that be helpful in modeling the original $x_i$'s, and how? Thanks!


title: Logistic vs. OLS

text: So the title is a bit misleading. I have a categorical dependent variable (no, maybe, yes) in multiple observations over a span of several years that I have turned into "percentage that said 'yes'" So technically it is now "yes" vs. "no and maybe" and bounded by 0 (0%) and 1 (100%). I am using OLS because I am interested in the one unit increase association between my independent (another percent) variable and my dependent variable. None of the coefficients of my independent variable would result in my dependent variable crossing 0 (0%) and 1 (100%). Additionally, besides the loss of interpret-ability, I don't want to use logistic regression because the "prevalence" of my independent variable is very high (>65% in most cases), so the resulting OR from a logistic regression would be inflated. What do you all say? How would you explain the use of OLS in this case if someone was very adamant that I have a dichotomous outcome variable?


title: Linear regression terminology question -- Beta (β)

text: I was a bit confused with the meaning of $\beta$, and thought its usage was rather loose. In fact, it seems that $\beta$ is used to express two distinct concepts:   1. The generalisation of the sample "b coefficient" to the population concerned.   2. The standardized regression coefficients (regression coefficients obtained when all variables are standardized with a sd of 1). Would there be an alternative symbol to any one of the two significations above, to avoid this confusion?


title: How do I calculate a probability to encounter a pattern of numbers in a numbers list?

text: Imagine you have an array of integer numbers. You need to calculate the probability of encountering some integer sub-array in that list. For example: **the original list:** '1284726437594563495834905834095845634853' which is N integers from 0 to 9. **the question is:** what is the probability to encounter '456', order is significant, in that big data set? **my considerations:** The are N-2 possible sets of 3 numbers in big list. Each has (1/10)^3 variants. With order significant we need one of them, but any place of the pattern would match. So we have               (1/10)^3 * (N-2)      In case order is not significant, there are n! possible permutations for n items. Thus the formula is               (1/10)^3 * 3! * (N-2)      Am I right? UPDATE: The digits are distributed randomly, i.e. a possibilities of them are equal. The list is some sequence of these random digits. By "encounter" I mean to have at least one '456' sequence in the given data set. Looks like my considerations were wrong. Here is a corrected version: A probability to meet at least one sub-array '456': A probability not to meet any is (1 - (1/10)^3) ^ (N-2), so probability to meet one is               1 - (1 - (1/10)^3) ^ (N-2)      A most probable number of patterns to meet is               (N-2) * 1/10^3


title: Why are Gaussian process models called non-parametric?

text: I am a bit confused. Why are Gaussian processes called non parametric models? They do assume that the functional values, or a subset of them, have a Gaussian prior with mean 0 and covariance function given as the kernel function. These kernel functions themselves have some parameters (i.e., hyperparameters). So why are they called non parametric models?


title: Determining predictive ability of models for readmission risks for patients

text: I am studying risk prediction models for readmissions to hospitals and during the course of my research, I came across this survey paper that evaluates various risk prediction models published in various medical journals: https://www.dropbox.com/s/t1inuv8112tnnip/Survey%20of%20Predictive%20Models%20for%20Readmission%20Risk.pdf Towards the end of the article (page 1696), the paper states that: > Even though the overall predictive ability of the clinical models was poor, > we did find that high- and low-risk scores were associated with a clinically > meaningful gradient of readmission rates. I am a bit confused by this statement. If in practice, a high risk score from a given model is associated with high rate of readmission, doesn't it imply that the model has good predictive ability?


title: Using a particle filter for robot localization

text: I have a robot that has a GPS and velocity sensors. The GPS updates roughly every 1-2 seconds. I've been playing around with a Kalman filter that has been working pretty well. I just learned and finally think I understand particle filters so I'm wondering if a particle filter be useful to keep track of the robot's location in between GPS updates instead of the Kalman filter. My plan would look something like this:   1. Starting GPS coordinates.    2. Create N random particles distributed around the starting coordinate (2 meters is the typical accuracy of most GPS sensors)   3. Robot moves and records velocity data from the sensors   4. Move all of the particles based on a linear model,velocity data, and noise   5. With the next GPS sensor update weight the particles based on a Gaussian from the updated coordinates. Am I on the right track or is this worth trying to code or should I just stick with the Kalman filter since it's a linear system?


title: How do i normalize residuals?

text: I am trying to adjust a hierarchical multiple regression model and no matter which transformations I use (z-transformation, sqrt, cuberoot, inv, inv sqrt ...), I do not manage to get the residuals normally distributed. What can I do? Does anyone have any suggestions? I also checked out a related thread, the transformations suggested there did not help me. My residuals are normal according to D'Agostino Normality Test, but not according to Shapiro-Wilk (which is the crucial one according to my supervisor). I cannot use a non- parametric model. I would appreciate your help a lot! Thanks, brobdingnag! qq plot of residual


title: Looking for a Good Text on Statistical Analysis of Satisfaction Data

text: I am looking for a good textbook (or other resource) that covers the analysis of satisfaction data. Most of my data uses likert-type scales. Can anyone recommend something with examples?


title: How to find the most important factors or combination of levels in a finite data set

text: I’d like to discover which “calendar factors” (e.g. day of the week, month) have the strongest relationship with whether or not a particular product is sold at least once. For days when the product was sold, I have data like this:               Date        Day of Week  Day Type  Month     01-01-2014  Wednesday    Weekday   January     01-13-2014  Monday       Weekday   January     02-09-2014  Sunday       Weekend   February     02-15-2014  Saturday     Weekend   February     02-16-2014  Sunday       Weekend   February     02-23-2014  Sunday       Weekend   February      If we pretend that this is the entire dataset, we see that the product was sold mostly on Sundays in February. I’ve come up with a simple way to rank each individual calendar factor (e.g. month) as well as a few combinations (e.g. day type + month) by creating a ratio of “sale days” vs. the entire population: using the example data, Wednesday + January would have a rank of 1 / 5. However, I’m struggling to determine which factors and/or combination of factors are most significant without simply picking the highest-ranking one. I feel like there must be a stronger mathematical approach to solving this kind of problem. Any ideas? Also, how would you classify this problem?


title: Can i use a mixture model for when I have an omitted variable?

text: I plan to fit a GAM or GAMM. There is one categorical variable which I think is important for explaining Y (or Y*), but it is not in my dataset - it is measurable but has not been measured. Can I use a mixture model (GAMM) to compensate for the omission of this variable? If I had this variable, I would have just used a GAM. Can you recommend some software for me to use? Do you have examples? Is the gamm4 R package "the way to go"? Other details: Y is count data, so I plan to use a Negative Binomial model or "family". The purpose of this taks is to investigate the relationships between the variables. Although the dataset, with N rows, fits into memory of my machine, a matrix with NxN rows will not fit into memory.


title: Given a series of Posterior probability scores, how can FDR be estimated?

text: I have a series of events which I manage to compute their significance using an specific software. Given that the software it outputs as a result the significant events followed by a posterior probability (PP) per event, is there any way to process those events using R to compute the False Discovery Rate (FDR) associated, and filter them to a specific FDR combined along with specific PP? I read some papers where they use Benjamini-Hochberg FDR calculation, but I found so many alternatives which have confused me a bit. As an example, this is my list of the first 20 PP events found by the software. Each event is determined by the columns 2 to 4, and the 6th one corresponds to the PP value.               space start end width PP     1 chr1 9950 10100 151 0.999999999996379     2 chr1 10400 11500 1101 1     3 chr1 13250 13650 401 0.999999999997592     4 chr1 13900 14400 501 0.999999999999933     5 chr1 14550 14650 101 0.991426544211918     6 chr1 15150 15250 101 0.99992525793782     7 chr1 15500 15750 251 1     8 chr1 17100 17250 151 0.999907293781992     9 chr1 17850 18000 151 0.999825155927197     10 chr1 28800 29900 1101 1     11 chr1 34800 34950 151 0.999999928381097     12 chr1 46900 47050 151 0.999984423295619     13 chr1 91350 91500 151 0.999957532522137     14 chr1 134550 134700 151 0.909017332110283     15 chr1 136850 136950 101 0.916775268418563     16 chr1 151200 151300 101 0.999895368100688     17 chr1 235850 235950 101 0.999185144466072     18 chr1 236450 236600 151 0.999984423299768     19 chr1 237700 237900 201 1     20 chr1 324900 325050 151 0.999065471728184      Thanks for your help!


title: Logistic regression diagnostics

text: In linear regression, we need to check which assumption is violated. But for logistic regression, what assumptions do we have for the model? I was attempting to wrap my head around this thing. But the diagnostics I can come up with right are DFBETAS plot and Pearson residual plot.


title: Error while calculating metrics like AUC, ACC in R

text: I am trying to make predictions on my validation data using a decision tree model created with training data set. I am able to do that with success but, I am not able to calculate various metrics like Area under curve (AUC), and overall classification rate (ACC) using mmetric() command in R. Here is the screen shot with errors and commands that I used. Note: Right click the image and select view image to see it clearly. "t" is my training data set and "v" is validation data set. ![mmetric](http://i.stack.imgur.com/20uib.jpg) Note: You can download my training data set "t" and validation data set "v" as a .csv files from here: Click here I used the same commands for a different data set and it worked perfectly, so could be this due to some problem with my data set. Thanks, P


title: What is the best way of presenting a correlation? (in R)

text: I have two variables which show significant correlation (Spearman). I would like to graphically show the strength of the relationship, much like showing a linear regression fit with confidence bands. What would be the best correct way of doing that?


title: Nested test to remedy pseudoreplication

text: **Data structure:** I have two datasets from two protected areas that differ in protection status. Both areas contain 43 and 37 sites each. **Question:** I would like to know which test would be the best for testing whether the PA status has had an effect on:   1. the first axis of a PCoA (principal coordinates analysis) - i.e. species composition turnover (derived by constructing a bray curtis dissimilarity matrix) and    2. species richness per site (a continuous variable).  **Problem:** I understand that there is pseudoreplication present in this as I only have two areas. From what I have read, it seems that I either have to use an ANCOVA / GLM / mixed-effect model, where I define PA status as both a random effect and a fixed effect. I intended to nest sites within PA, but it seems that as there is only one datapoint per site it will not work as a nested object. For those familiar with R, here are some codes I have tried:               pcoaPAanovadata1 <- read.csv("PCoA\\data\\                                   combined data PCoA axis 1 with distance variables.csv",                                    header=T)          str(pcoaPAanovadata1)     'data.frame': 80 obs. of 7 variables:     PCOA:    num -0.2215 -0.3521 -0.0611 0.3434 -0.3624 ...     PA.stat: Factor w/ 2 levels "N","P": 1 1 1 1 1 1 1 1 1 1 ...     village: num 33.6 33.7 39.9 37.9 34 ...     road:    num 4.18 3.8 0.89 0.1 3.43 5.49 1.86 5.04 0.79 0.88 ...     track:   num 8.11 6.48 3.11 2.71 4.49 5.35 1.25 4.03 7.62 6.77 ...     site:    Factor w/ 80 levels "M1_11","M1_17",..: 1 2 3 4 5 6 7 8 9 10 ...     rich:    num 3.27 1.79 7.31 0.82 1.79 1.82 2.45 0.82 5.47 2.79 ...      compare community composition turnover at different PAs: below specifies a null model where the slope deviates as a result of the random effect               z0 <- lmer(rich ~ 1, random = ~ 1 | pastat/site, data = pcoaPAanovadata1)     summary(z0)     z1 <- lme(rich ~ pastat, random = ~ 1 | pastat/site, data = pcoaPAanovadata1)     summary(z1)     anova(z0,z1)      impacts of distance variables:               zz <- lme(pcoa ~ road, random = ~ 1 | pastat/site, data = pcoaPAanovadata1)     summary(zz)      The errors I get from the lme(linear mixed effect model):               Warning message:     In pt(-abs(tTable[, "t-value"]), tTable[, "DF"]) : NaNs produced      The error I get from the ANOVA:               Warning message:     In anova.lme(z0, z1) :     fitted objects with different fixed effects. REML comparisons are not meaningful.      Firstly, I was hoping to just clarify whether the test I am running is correct. Secondly, it'd be great if someone could tell me what the errors mean. I apologise if my question is poorly phrased, I am relatively new to R and the statistics I am using.


title: calculation of correlation time?

text: I have two time series of a two realizations of a random process, the process fluctuate around an average so i consider that is a ergodic process. I compute the time correlation function, but i need the characteristic time scale, better known as the correlation time, but i don't know how, This time characterizes the time for time correlation function decay to zero. I tried to fit a curve, but the resulting time its too big, so i think its not the correct way.


title: Interpreting the effects of a covariate in a 2 x 2 x continuous x continuous ANOVA?

text: As part of my social psych experiment, I conducted a 2(Alcohol norm: 'Low' vs. Control) x 2(Social life norm: 'Positive' vs Control) x Group identification (continuous, mean-centred) x Attitude (continuous, mean-centred) ANOVA on the dependent variable (alcohol consumption), and included a covariate (social pressure). I'm getting pretty confused at how to interpret my results, and which statistical tests to run next: 1) The effect of my covariate (social pressure) was significant, so I ran a simple linear regression to see the general relationship between social pressure and alcohol consumption. Is that ok? 2) Including the covariate in the ANOVA did not significantly alter any of the main effects from the same analysis without the covariate, so I concluded that social pressure did not predict alcohol consumption above and beyond the effects of the interaction between group identification, attitude, alcohol norm, and social life norm. That seems right in stats term, can anyone help me understand what that means to the layman? 3) Whilst trying to interpret my conclusion from Q2, I decided that maybe social pressure and alcohol norm weren't independent. When I included them in that ANOVA, would I have had to assume that they were independent? I conducted an independent samples t-test on social pressure with alcohol norm as the grouping variable, which confirmed that social pressure was significantly higher in the 'control' condition for alcohol norm than the 'low' condition. So when interpreting the main results, I concluded that it was possible social pressure was not significant because it was already captured in alcohol norm, and thus did not have any additional effect. Is this a valid conclusion? If not, could someone explain why? Also, would I have to make post-hoc adjustments for the t-test? Sorry for all the questions, but as you can tell I'm getting pretty confused here. Help to any or all of these questions would be massively appreciated


title: Posterior probability using forward backward algorithm in R

text: For one of my projects I need to find posterior probability of visiting a state S and emitting a symbol. I have built a HMM in R and later I get one observation sequence. But, I am not able to interpret the output.               > # Sequence of observations     > observations = c("L","L","R","R")     > # Calculate posterior probablities of the states     > posterior = posterior(hmm, observations) ## ---> hmm is hidden markov model     > print(posterior)             index     states         1       2       3         4          A 0.6037344 0.56639 0.43361 0.3962656          B 0.3962656 0.43361 0.56639 0.6037344      How do I interpret this? This matrix gives output w.r.t time. But, I need posterior probability of visiting a state S and emitting a symbol say A. How can I get it? * * * Thanks for the response.The description says it combines forward and backward probabilities to get posterior details...But,as i have read in papers, the forward-backward algorithm gives the probability of being in a state say A and emitting a symbol say L.. which i do not see here. Here is all the details of HMM - # Initialise HMM # init HMM format: A and B are hidden states,L and R are observations,transprob: initial # transition probabilities and emissionprobs IS EMISSION probability hmm = initHMM(c("A","B"), c("L","R"), transProbs=matrix(c(.8,.2,.2,.8),2), emissionProbs=matrix(c(.6,.4,.4,.6),2)) print(hmm) # Sequence of observations observations = c("L","L","R","R") # Calculate posterior probablities of the states posterior = posterior(hmm,observations) print(posterior) Description of POSTERIOR: This function computes the posterior probabilities of being in state X at time k for a given sequence of observations and a given Hidden Markov Model. Details The posterior probability of being in a state X at time k can be computed from the forward and backward probabilities: Ws(X_k = X | E_1 = e_1, ... , E_n = e_n) = f[X,k] * b[X,k] / Prob(E_1 = e_1, ... , E_n = e_n) Where E_1...E_n = e_1...e_n is the sequence of observed emissions and X_k is a random variable that represents the state at time k. But,in general how can we calculate posterior probability without using this algo


title: Does autocorrelation cause bias in the regression parameters in piecewise regression?

text: In simple linear regression problems, autocorrelated residuals are supposed not to result in biased estimates for the regression parameters. Can the same be said for piecewise regression? Suppose I want to fit a continuous, piecewise linear function of a single variable. Let's say for example we have data on shipping cost and weight of shipment. The function is piecewise because as the weight increases, at some point an additional rail car is required. We want to find the breakpoints and the slopes of the individual pieces. The model is fit, and for whatever reason, the residuals are found to be serially correlated in time. Could the regression parameters be biased? I have posted some data in a Google spreadsheet at this link: http://goo.gl/LrTv3 Suppose it is known that there are two breakpoints at (unknown) points x1 and x2. We want to fit the data to a model f(x) given by:               x < x1:      f(x) = a + m1*x     x1 < x < x2: f(x) = a + m1*x1 + m2*(x - x1)     x > x2:      f(x) = a + m1*x1 + m2*(x2 - x1) + m3*(x - x2)      I use the nlm function in R to find the unknown parameters x1, x2, m1, m2 and m3:               sqerr <- function(prm,y,x) {       a <- prm[1]       x1 <- prm[2]       x2 <- prm[3]       m1 <- prm[4]       m2 <- prm[5]       m3 <- prm[6]       sqerr <- sum((y-(a+ifelse(x<x1,m1*x,                     m1*x1+ifelse(x<x2,m2*(x-x1),                     m2*(x2-x1)+m3*(x-x2)))))^2)     }     data <- read.table("data.txt",header=T)     ai <- 0.4; x1i <- 0.4; x2i <- 0.7; m1i <- 0.0; m2i <- 0.8; m3i <- 3     prm <- c(ai,x1i,x2i,m1i,m2i,m3i)     uu <- nlm(sqerr,prm,data$Y,data$X)      Then I plot the residuals vs. the lag-1 residuals:               y <- data$Y         x <- data$X     a <- uu$est[1]         x1 <- uu$est[2]     x2 <- uu$est[3]         m1 <- uu$est[4]     m2 <- uu$est[5]         m3 <- uu$est[6]     resid <- (y-(a+ifelse(x<x1,m1*x,m1*x1+ifelse(x<x2,m2*(x-x1),m2*(x2-x1)+m3*(x-x2)))))     plot(resid[1:149]~resid[2:150])      There is clearly some sequential correlation. So my question is, are the regression parameters biased because of this? I have an old paper by Kadiyala (A Transformation Used to Circumvent the Problem of Autocorrelation, Econometrica Vol. 36, No. 1, Jan. 1968) that states: > "It is well known (see Watson [7] and Watson and Hannan [8]) that simple > least squares estimators, though unbiased (when the independent variables > are "fixed variates"),are, in general inefficient in the presence of > autocorrelation among the disturbances." It seems that by "simple least squares" he means linear equations of the form y = a + bx (that is the example used in the paper). But I have seen other papers that seem to imply that the estimators (i.e., regression parameters) are unbiased no matter what type of model you have. I don't think it's true in general.


title: Justification for use of $\chi^2(1)$ in Wald and score test

text: In a recent exam, we were asked to justify the use of the $\chi^2(1)$ distribution in performing the Wald or Rao's score test. There was only 1 mark for this (approx. 2.5 mins worth of time). My answer was > The Wald and score test statistics are based on various approximations to > the log-likelihood ratio, which are valid and equivalent in large samples > when $H_0$ is true. For example, the approximation for the 2nd Wald > statistic is > > $$ 2\log(LR)\simeq > (\hat{\theta}_n-\theta_0)^2E\\{-\ell''(\theta)\\}|_{\theta_0} > =(\hat{\theta}_n-\theta_0)^2I(\theta_0) $$ > > where $I(\theta_0)$ is the Fisher information at $\theta_0$. Then, using > asymptotic normality, > > $$\hat{\theta}_n\approx N \left ( \theta_0,\frac{1}{ni(\theta_0)}\right )=N > \left ( \theta_0,\frac{1}{I(\theta_0)}\right )$$ > > which yields $(\hat{\theta}_n-\theta_0)\sqrt{I(\theta_0)} \approx N(0,1)$ > > Since $2\log(LR)$ is approximately equal to the square of the LHS of this > when $n$ is large, it is approximately distributed as the square of a > standard normal random variable, that is, as $\chi^2(1)$ The marker wrote "insufficient" and I got zero for this. As it is a summative exam, they won't give any feedback, nor engage in any discussion about it. I was wondering if anyone here can explain what I have missed or where I went wrong. I'm not very good with latex so I hope I didn't make any mistakes in the typing ! This is for an elective module in statistical theory in the final year of an undergraduate maths degree. Thanks ! Edit: There is a formal procedure to have my script remarked, but for the sake of 1 mark, and since I passed quite comfortably, I don't really want to rock the boat.


title: Bayesian inference with the change of variables

text: I have a following question which I haven't been able to solve on my own. Imagine that I have measured some property X from data and for simplicity the posterior is Gaussian: $P(X|D) \sim Normal(X_0, \Sigma) $ Now I assume that I have some particular model M with parameters $\phi$ where $X=F(\phi)$ and $F()$ is an arbitrary function. And now I'm trying to express the inference on $\phi$, e.g. $P(\phi|M D)$ I'm particuarly interested in the case where the $F$ mapping isn't simple $R^n$ -> $R^n$ and doesn't cover the whole space of X. E.g. in 2-D case $x_1=cos(\phi)$ and $x_2=sin(\phi)$ I'm interested in a derivation from basic rules a formula for $P(\phi| M, D) $. From naive intuition My guess the formula will be $P(\phi | M, D) \sim P(\phi |M) * Gaussian(F(\phi)-X_0,\Sigma)$ But I can't really derive it from Cox rules. I would appreciate any help with this, it should not be very complicated, I think.


title: Optimal non-linear IV estimator

text: I know that the feasible optimal GMM estimator is consistent and asymptotically efficient. I also know that the fully-parametric MLE estimator is more efficient than GMM provided that we know the distribution of the data. Does this mean that in an ideal world, to obtain the "optimal" estimates for a nonlinear model with instrumental variables, one should use the MLE approach (provided that we know the relevant data distributions)? If so, if we have more instruments than endogenous variables, can we still use the MLE method? Thanks!


title: How to check out an appropriate survival model?

text: Maybe this question is too general, but I think it's worthy to post to clarify me some standard procedures in survival analysis. As we know, there are a list of method and models like parametric (exp, weibull, AFT), nonparametric (KW, NA), and semi-parametric (cox prop, prop) methods in survival analysis. Suppose I have a data frame like this:               set.seed(123)     require(survival)     df<-data.frame(time=as.integer(rnorm(100,50,5)),                     status=rbinom(100,1,0.7),                     age=rnorm(100,60, 5),                     gender=rbinom(100,1,0.5))     df$time<-ifelse(df$time>=50,50,df$time)     df         time status   age gender     1     47      1 63.94      0     2     48      0 63.85      0     3     50      1 61.66      0     ...     98    50      1 53.20      0     99    48      1 56.68      1     100   44      0 62.43      0      I'd like to know what's an appropriate way to check which model I should use for this survival analysis.


title: Is it good practice to standardize your data in a regression with panel/longitudinal data?

text: In general, I standardize my independent variables in regressions, in order to properly compare the coefficients (this way they have the same units: standard deviations). However, with panel/longitudinal data, I'm not sure how I should standardize my data, especially if I estimate a hierarchical model. To see why it can be a potential problem, assume you have $i = 1, \ldots, n$ individuals measured along $t=1,\ldots, T$ periods and you measured a dependent variable, $y_{i,t}$ and one independent variable $x_{i,t}$. If you run a complete pooling regression, then it's ok to standardize your data in this way: $x.z = (x- \text{mean}(x))/\text{sd}(x)$, since it will not change t-statistic. On the other hand, if you fit an unpooled regression, i.e., one regression for each individual, then you should standardize your data by individual only, not the whole dataset (in R code):               for (i in 1:n) {       for ( t in 1:T) x.z[i] =  (x[i,t] - mean(x[i,]))/sd(x[i,])      }      However, if you fit a simple hierarchical model with a varying intercept by individuals, then you are using a shrinkage estimator, i.e, you are estimating a model between pooled and unpooled regression. How should I standardize my data? Using the whole data like a pooled regression? Using only individuals, like in the unpooled case?


title: How to calculate the variance of a partition of variables

text: I'm running an experiment where I'm gathering (independent) samples in parallel, I compute the variance of each group of samples and now I want to combine then all to find the total variance of all the samples. I'm having a hard time finding a derivation for this as I'm not sure of terminology. I think of it as a partition of one RV. So I want to find $Var(X)$ from $Var(X_1)$, $Var(X_2)$, ..., and $Var(X_n)$, where $X$ = $[X_1, X_2, \dots, X_n]$. EDIT: The partitions are not the same size/cardinality, but the sum of the partition sizes equal the number of samples in the overall sample set. EDIT 2: There is a formula for a parallel computation here, but it only covers the case of a partition into two sets, not $n$ sets.


title: ratio, MSE of basic demand to MSE of aggregated demand for MA(1) process

text: I have a basic demand series that follow MA(1) process,I've applied non- overlapping aggregation approach and then SES on both basic and aggregated series to obtain forecasts, and then I disaggregated forecasts resulted from aggregated series and finally I calculate ratio of MSEbasic/MSEdagg, where MSEbasic=Var(dt-ft) and MSEdagg=Var(dt-(FT/m), dt is basic series at time t, ft is forecast of basic series, FT is aggregated forecast and m is aggregation level and alpha is smoothing constant. this ratio shows that for highly negative value of Theta, basic approach works better and for the rest of Theta parameter aggregated approach works better, I'm trying to understand why aggregation does not work for high negative Theta,I was wondering if you guide me. bests Roji![enter image description here](http://i.stack.imgur.com/1jq9d.jpg)


title: What form of analysis should be used for loss given default estimates?

text: I am working on BASEL II IRB models and we have to estimate loss based on historic defaults. There are different outcomes/scenarios we have identified that a default can encounter that will affect the LGD like CURE, FORCED SALE, LOSS, PARTIAL CURE. We have the data needed to be able to allow the classification to happen, but I am unsure of what form of analysis is best suited to process this data and produce the necessary segmentations and subsequent estimates of loss for each one. We have data currently in a format with one row per month per default for every month including and since the account first defaulted. A number of scenarios involve activity over a number of months for a defaulted account. I am unsure what form of analysis I shoud be performing. I considered classification trees but as a criteria can span multiple months (in rows) I don't think this would work and I also considered time series analysis but this is a new area for me and I'm not sure whether it can pull out events. What form of analysis should I be using and are there any specific resources for it in R that you think would be helpful?          **UPDATE** I have reduced the dataset to a series of 6 month snapshots and classifed each account at each position such that I have a dataset like `ID,Age,LGDstatus,LTV,BTV,[OTHER CHARACTERISTICS],Bal,BalChangeSinceInception` where a predictive model needs to be built such that we can predict how much BalChangeSinceInception is expected for each LGDstatus at each point, and the likelihood that an account will go to the various LGDstatus. Simplisticly then I could build a prediction which says... An account has gone into default, based on these x factors, by month Y we expect a loss of z% based on an a% probability of experiencing LOSS and a b% probability of experiencing a CURE However, that may be the wrong approach and direction would be invaluable.


title: P-value adjustment in correlational analysis

text: I'm currently doing some large correlational analysis: Let's just say I wanted every kendall-correlation and it's associated p-value from 20x30 Items in a neat matrix - should I be doing p-value adjustment correcting for multiple comparisons, eg. holm (nothing's significant) or should I just use the passthrough option adjustment="none" (lots of significance) and be done with it? Sample Size is 145-153 depending on pairwise deletion, Korrelations are quite low(-.14/+28), using R with corr.test from psych-package Goal is exploring possible associative relationships between items Thanks in advance! JS


title: Gauss-Markov assumptions

text: I was trying to determine which of the Gauss-Markov assumptions allow us to see that $b_1$ is an unbiased estimator of $\beta_1$. I have a feeling it's that $X_{i}$ is not random, but is there anything else that I'm missing?


title: How to perform model selection in GEE in R

text: I would like to do model selection for generalized estimating equations (GEE). Pan (2001) is most frequently cited for developing a method using QIC. I am wondering if anyone knows of a way to do this in R? I am currently using the package 'geepack' for my GEE analysis. I heard about 'yags' but it is not available for R version 2.13 or 2.14 and because of other things I'm doing, I'd rather not try to revert to an older version of R. Any Help would be appreciated. Thanks, Dan


title: Multiple regression with both positive and negative data

text: I have a problem regarding interpretation of a multiple regression. For $Y=a+bX_1+cX_2$, the influence of independent variable $X_1$ on $Y$, the dependent variable, is positive (coefficient $b$ is positive). This means that an increase of one unit of $X_1$ leads to an increase of one unit of $Y$, which is consistent with the economic interpretation. But in my dataset, $X_1$ has also negative values. For those values, the interpretation is not correct anymore. An increase with 1 unit of $X_1$ will lead to a decrease of Y, meaning an inverse relation. What am I doing wrong? Is it better to use the data as absolute values ABS(X)?


title: Is the Gaussian copula (for d=2) with normal margins identical to the bivariate normal?

text: I am not sure about this: In the 2 dimensional case, if I consider the Gaussian copula, is this identical to the bivariate normal distribution, in the case I choose the normal distribution for the margins? Does this hold for multidemsionality? So is a d-dimensional Gaussian copula with normal margins identical to the multivariate normal distribution?


title: Why is this representing the left tail?

text: In this source about the Clayton copula on page 18 they write: > It has been used to study correlated risks because it exhibits strong left > tail dependence and relatively weak right tail dependence. Anecdotal and > empirical evidence suggests that loan defaults are highly correlated during > recessionary times. Similarly, researchers have studied the “broken heart > syndrome” in which spouses’ ages at death tend to be correlated. When > correlation between two events, such as performance of two funds or spouses’ > ages at death, is strongest in the left tail of the joint distribution, > Clayton is an appro- priate modeling choice. Now, I am not understanding, why these examples mean, that the left tail is correlated? Why is the spouses' ages at death the left tail? And why is recessionary times the left tail? Is it, because in resessionary times negative returns occur and this is representing the left tail?


title: Term frequency/inverse document frequency (TF/IDF): weighting

text: I've got a dataset which represents 1000 documents and all the words that appear in it. So the rows represent the documents and the columns represent the words. So for example, the value in cell $(i,j)$ stands for the times word $j$ occurs in document $i$. Now, I have to find 'weights' of the words, using tf/idf method, but I actually don't know how to do this. Can someone please help me out?


title: Cluster-robust estimation versus Random Effects estimation

text: My dependent variable is the decision of individuals for a value between 1 and 8 in an economic experiment. Every individual repeated that decision 20 times. I want to regress individual characteristics on that choice. I expect that the decisions of one individual are correlated with each other. So I need some kind of panel model. Should I go for random effects or cluster-robust estimation? RE seems to be prevalent in the literature, but cluster-robust was a suggestion by my supervisor, but it was simply her intuition. What are the reasons for either? What should I pick?


title: Supervised classification with different numbers of predictors per observation

text: I would like to build a supervised classifier on data that contains a different number of predictors depending on the observation, and I'm looking for the best approach to use. For example (which isn't quite my real dataset - mine is biology), imagine I want to build a learner for people's road accident risk based on the features of their cars. Clearly people can have more than one car (and the number of cars can by itself be a predictor), and I'm wondering what approaches exist for using the features of each of their cars as predictors. I can of course just use some kind of pooling strategy for the features (for example, counting the total number of air fresheners and having a binary variable describing whether or not at least one car has a spoiler), but what if, say, the least accident-prone people have one car that's a convertible and one that's a family van, with a different set of features in each? Pooling strategies won't describe this situation precisely, leading to a 'dilution' of the best set of predictors. Another strategy that comes to mind is models with missing data - in which case I could have a variable for each feature of each type of car, which could be set to NA for people who don't own a certain type. But what if somebody has two cars of the same type and this, in turn, is different from having just one with respect to their accident risk? In addition, random forests, for example, while allowing for missing data, lose the ability to assess variable importance when it is present - which isn't great in my situation. Surely people across disciplines - including actuarial science I suppose - have come across this type of problems? I would greatly appreciate if you could point me to some relevant papers/books.


title: Is the explanatory variable in the generalized linear model viewed as a random variable or parameter?

text: In the generalized linear model, the distribution of a random variable $Y$ is assumed to be an exponential family distribution, and it is written in terms of an explanatory variable $X$ and parameter $\theta$. I was wondering if $X$ is viewed as a parameter of the distribution of $Y$, or a random variable and the distribution of $Y$ mentioned earlier is in fact a conditional one given $X$? If $X$ can be viewed as a random variable, there seems to me a problem here. In the generalized linear model, the distribution of a random variable $Y$ is assumed to be an exponential family distribution. But since the distribution of $Y$ depends on $X$, it is actually that $P(Y|X=x)$ is an exponential family distribution. $X$ can have any distribution, and the unconditional distribution of $Y$ may not be an exponential family distribution. So how shall we understand the problem? Thanks!


title: Paired t-test with subsamples in R?

text: I've been trying to find a similar "problem", but I didn't, so I'm going to try to explain what am I looking for. I have for different sites A,B,C and D. In each site I have a forest and a fallow=> AFo, AFa, BFo,BFa... In each forest and fallow I have 4 quadrates: AFo1, AFo2, AFo3, AFo4, AFa1,... I measured different variables such as pH, electric conductivity, total phosphorus, etc in all quadrates. My aim is to compare the forests to the fallows, but I have to take into account the variability between the sites and the quadrates. I have made a two way crossed ANOVA, but my teacher told me that I should do a t-test which takes into account the subsamples (the quadrates). My fixed factor would be forest/fallow, and my random factor the sites. Can anyone help me, please? :) Thank you very much in advance! Lili * * * More or less, but normally she wants me to perform a t-test on all the AFo-s , while taking into account the subsamples 1,2,3 and 4. Because what I would like to do is not to compare the quadrate1 from the forest with quadrate 1 from the fallow, but the mean of the 4 quadrates of the forest with the mean of the 4 quadrates of the fallow. Hm. don't know whether I'm expressing myself clearly..


title: Sensitivity Analysis for Mulitple Dependent and Independent Variables

text: I have Mulitple Vectors of Independent Variables and Multiple Vectors of Dependent Variables from a Monte Carlo Simulation. For each Dependent Variable, I want to know the contribution of each Independent Variable to the uncertainty of the Dependent Variable. I have the MATLAB Statistics Toolbox. I am looking for the general approach. The specific approach with MATLAB Stats Toolbox would be a PLUS.


title: Is there a name for this type of bootstrapping?

text: Consider an experiment with multiple human participants, each measured multiple times in two conditions. A mixed effects model can be formulated (using lme4 syntax) as:               fit = lmer(         formula = measure ~ (1|participant) + condition     )      Now, say I want to generate bootstrapped confidence intervals for the predictions of this model. I think I've come up with a simple and computationally efficient method, and I'm sure I'm not the first to think of it, but I'm having trouble finding any prior publications describing this approach. Here it is:   1. Fit the model (as above), call this the "original model"   2. Obtain predictions from the original model, call these the "original predictions"   3. Obtain residuals from the original model associated with each response from each participant   4. Resample the residuals, sampling _participants_ with replacement   5. Fit a linear mixed effects model with gaussian error to the _residuals_ , call this the "interim model"   6. Compute predictions from the interim model for each condition (these predictions will be very close to zero), call these the "interim predictions"   7. Add the interim predictions to the original predictions, call the result the "resample predictions"   8. Repeat steps 4 through 7 many times, generating a distribution of resample predictions for each condition from which once can compute CIs. I've seen "residual bootstrapping" procedures in the context of simple regression (i.e. not a mixed model) where residuals are sampled as the unit of resampling and then added to the predictions of the original model before fitting a new model on each iteration of the bootstrap, but this seems rather different from the approach I describe where residuals are never resampled, people are, and only _after_ the interim model is obtained do the original model predictions come into play. This last feature has a really nice side- benefit in that no matter the complexity of the original model, the interim model can always be fit as a gaussian linear mixed model, which can be substantially faster in some cases. For example, I recently had binomial data and 3 predictor variables, one of which I suspected would cause strongly non- linear effects, so I had to employ Generalized Additive Mixed Modelling using a binomial link function. Fitting the original model in this case took over an hour, whereas fitting the gaussian LMM on each iteration took mere seconds. I really don't want to claim priority on this if it's already a known procedure, so I'd be very grateful if anyone can provide information on where this might have been described before. (Also, if there are any glaring problems with this approach, do let me know!)